--- title: "DSolve" language: "en" type: "Symbol" summary: "DSolve[eqn] solves a differential equation eqn. DSolve[eqn, u, x] solves a differential equation for the function u, with independent variable x. DSolve[eqn, u, {x, xmin, xmax}] solves a differential equation for x between xmin and xmax. DSolve[{eqn1, eqn2, ...}, {u1, u2, ...}, ...] solves a list of differential equations. DSolve[eqn, u, {x1, x2, ...}] solves a partial differential equation. DSolve[eqn, u, {x1, x2, ...} \\[Element] \\[CapitalOmega]] solves the partial differential equation eqn over the region \\[CapitalOmega]." keywords: - ODE - DAE - DDE - PDE - Abel equation - Abramov algorithm - Airy's differential equation - analytical solution - Bernoulli equation - Bocharov techniques - Bronstein algorithm - closed-form solution - diff - diffeq - differential-algebraic equations - differential equations - Hamilton-Jacobi equations - homogeneous differential equations - inhomogeneous differential equations - Kamke - Kovacic algorithm - linear differential equations - Mellin transform methods - multivariate differential equations - nonlinear differential equations - ordinary differential equations - partial differential equations - piecewise differential equations - Riccati equations - solving differential equations - symmetry reduction - hyperbolic PDE - parabolic PDE - elliptic PDE - heat equation - diffusion equation - wave equation - Laplace equation - potential equation - Schrodinger equation - Burgers' equation - Sturm-Liouville problem - eigenvalue problem - Poisson equation - Helmholtz equation - Black-Scholes equation - integral equation - integro-differential equation - Volterra integral equation - Fredholm integral equation - Abel integral equation - singular integral equation - weakly singular integral equation - dsolve - pdsolve canonical_url: "https://reference.wolfram.com/language/ref/DSolve.html" source: "Wolfram Language Documentation" related_guides: - title: "Equation Solving" link: "https://reference.wolfram.com/language/guide/EquationSolving.en.md" - title: "Calculus" link: "https://reference.wolfram.com/language/guide/Calculus.en.md" - title: "Differential Equations" link: "https://reference.wolfram.com/language/guide/DifferentialEquations.en.md" - title: "Differential Operators" link: "https://reference.wolfram.com/language/guide/DifferentialOperators.en.md" - title: "Partial Differential Equations" link: "https://reference.wolfram.com/language/guide/PDEModelingAndAnalysis.en.md" - title: "Mathematical Data" link: "https://reference.wolfram.com/language/guide/MathematicalData.en.md" - title: "Fractional Calculus" link: "https://reference.wolfram.com/language/guide/FractionalCalculus.en.md" - title: "Scientific Models" link: "https://reference.wolfram.com/language/guide/ScientificModels.en.md" - title: "Integral Transforms" link: "https://reference.wolfram.com/language/guide/IntegralTransforms.en.md" - title: "Symbolic Vectors, Matrices and Arrays" link: "https://reference.wolfram.com/language/guide/SymbolicArrays.en.md" - title: "Solvers over Regions" link: "https://reference.wolfram.com/language/guide/GeometricSolvers.en.md" related_functions: - title: "DSolveValue" link: "https://reference.wolfram.com/language/ref/DSolveValue.en.md" - title: "NDSolve" link: "https://reference.wolfram.com/language/ref/NDSolve.en.md" - title: "AsymptoticDSolveValue" link: "https://reference.wolfram.com/language/ref/AsymptoticDSolveValue.en.md" - title: "WhenEvent" link: "https://reference.wolfram.com/language/ref/WhenEvent.en.md" - title: "DEigensystem" link: "https://reference.wolfram.com/language/ref/DEigensystem.en.md" - title: "DEigenvalues" link: "https://reference.wolfram.com/language/ref/DEigenvalues.en.md" - title: "DFixedPoints" link: "https://reference.wolfram.com/language/ref/DFixedPoints.en.md" - title: "DStabilityConditions" link: "https://reference.wolfram.com/language/ref/DStabilityConditions.en.md" - title: "NDEigensystem" link: "https://reference.wolfram.com/language/ref/NDEigensystem.en.md" - title: "NDEigenvalues" link: "https://reference.wolfram.com/language/ref/NDEigenvalues.en.md" - title: "GreenFunction" link: "https://reference.wolfram.com/language/ref/GreenFunction.en.md" - title: "CompleteIntegral" link: "https://reference.wolfram.com/language/ref/CompleteIntegral.en.md" - title: "Solve" link: "https://reference.wolfram.com/language/ref/Solve.en.md" - title: "RSolve" link: "https://reference.wolfram.com/language/ref/RSolve.en.md" - title: "Integrate" link: "https://reference.wolfram.com/language/ref/Integrate.en.md" - title: "DifferentialRoot" link: "https://reference.wolfram.com/language/ref/DifferentialRoot.en.md" - title: "StreamPlot" link: "https://reference.wolfram.com/language/ref/StreamPlot.en.md" - title: "ItoProcess" link: "https://reference.wolfram.com/language/ref/ItoProcess.en.md" - title: "SystemModelSimulate" link: "https://reference.wolfram.com/language/ref/SystemModelSimulate.en.md" - title: "FractionalD" link: "https://reference.wolfram.com/language/ref/FractionalD.en.md" - title: "CaputoD" link: "https://reference.wolfram.com/language/ref/CaputoD.en.md" - title: "TruncateSum" link: "https://reference.wolfram.com/language/ref/TruncateSum.en.md" related_tutorials: - title: "Differential Equations" link: "https://reference.wolfram.com/language/tutorial/Calculus.en.md#9252" - title: "Differential Equation Solving with DSolve" link: "https://reference.wolfram.com/language/tutorial/DSolveOverview.en.md" - title: "Symbolic Solutions of PDEs" link: "https://reference.wolfram.com/language/tutorial/SymbolicSolutionsOfPDEs.en.md" - title: "Introduction to Fractional Calculus" link: "https://reference.wolfram.com/language/tutorial/FractionalCalculus.en.md" - title: "Numerical Differential Equation Solving with NDSolve" link: "https://reference.wolfram.com/language/tutorial/NDSolveOverview.en.md" - title: "Implementation notes: Algebra and Calculus" link: "https://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.en.md#8313" --- # DSolve DSolve[eqn] solves a differential equation eqn. DSolve[eqn, u, x] solves a differential equation for the function u, with independent variable x. DSolve[eqn, u, {x, xmin, xmax}] solves a differential equation for x between xmin and xmax. DSolve[{eqn1, eqn2, …}, {u1, u2, …}, …] solves a list of differential equations. DSolve[eqn, u, {x1, x2, …}] solves a partial differential equation. DSolve[eqn, u, {x1, x2, …}∈Ω] solves the partial differential equation eqn over the region Ω. ## Details and Options * ``DSolve`` can solve ordinary differential equations (ODEs), partial differential equations (PDEs), differential algebraic equations (DAEs), delay differential equations (DDEs), integral equations, integro-differential equations, and hybrid differential equations. * Different classes of equations solvable by ``DSolve`` include: | | | | ------------------------------ | ------------------------------- | | u'[x] == f[x, u[x]] | ordinary differential equation | | a ∂xu[x, y] + b ∂yu[x, y] == f | partial differential equation | | f[u'[x], u[x], x] == 0 | differential algebraic equation | | u'[x] == f[x, u[x - x1]] | delay differential equation | | u'[x] + ∫abk[x, t]u[t]\[DifferentialD]t == f | integro-differential equation | | {…, WhenEvent[cond, u[x] -> g]} | hybrid differential equation | * The output from ``DSolve`` is controlled by the form of the dependent function ``u`` or ``u[x]`` : | | | | | -------------------- | ------------------ | -------------------------------- | | DSolve[eqn, u, x] | {{u -> f}, …} | where f is a pure function | | DSolve[eqn, u[x], x] | {{u[x] -> f[x]}, …} | where f[x] is an expression in x | * With a pure function output, ``eqn /. {{u -> f}, …}`` can be used to verify the solution. » * The dependent variable in ``DSolve[eqn]`` is inferred from ``eqn`` and may be specified either explicitly as ``u[x]`` or as a pure function ``u`` for autonomous equations. » * ``DSolve`` can give implicit solutions in terms of ``Solve``. » * ``DSolve`` can give solutions that include ``Inactive`` sums and integrals that cannot be carried out explicitly. Variables ``K[1]``, ``K[2]``, ``…`` are used in such cases. * Boundary conditions for ODEs and DAEs can be specified by giving equations at specific points such as ``u[x1] == a``, ``u'[x2] == b``, etc. * Boundary conditions for PDEs can be given as equations ``u[x, y1] == a``, ``Derivative[1, 0][u][x, y1] == b``, etc. or as ``DirichletCondition[u[x, y] == g[x, y], cond]``. * Initial conditions for DDEs can be given as a history function ``g[x]`` in the form ``u[x /; x < x0] == g[x]``. * ``WhenEvent[event, action]`` may be included in the equations ``eqn`` to specify an ``action`` that occurs when ``event`` becomes ``True``. * The specification ``u∈Vectors[n]`` or ``u∈Matrices[{m, n}]`` can be used to indicate that the dependent variable ``u`` is a vector-valued or a matrix-valued variable, respectively. Alternatively, ``u`` can be specified as a ``VectorSymbol`` or ``MatrixSymbol``. » » * The region ``Ω`` can be anything for which ``RegionQ[Ω]`` is ``True``. * ``N[DSolve[…]]`` calls ``NDSolve`` or ``ParametricNDSolve`` for differential equations that cannot be solved symbolically. * The following options can be given: | | | | | ------------------------ | ------------- | --------------------------------------- | | Assumptions | \$Assumptions | assumptions on parameters | | DiscreteVariables | {} | discrete variables for hybrid equations | | GeneratedParameters | C | how to name generated parameters | | Method | Automatic | what method to use | | IncludeSingularSolutions | False | whether to include singular solutions | * ``GeneratedParameters`` controls the form of generated parameters; for ODEs and DAEs these are by default constants ``C[n]`` and for PDEs they are arbitrary functions ``C[n][…]``. » --- ## Examples (185) ### Basic Examples (3) Solve a differential equation with dependent variable y[x]: ```wl In[1]:= DSolve[y'[x] + y[x] == a Sin[x], y[x], x] Out[1]= {{y[x] -> E^-x C[1] + (1/2) a (-Cos[x] + Sin[x])}} ``` Include a boundary condition: ```wl In[2]:= DSolve[{y'[x] + y[x] == a Sin[x], y[0] == 0}, y[x], x] Out[2]= {{y[x] -> -(1/2) a E^-x (-1 + E^x Cos[x] - E^x Sin[x])}} ``` --- Get a "pure function" solution for ``y`` : ```wl In[1]:= DSolve[{y'[x] + y[x] == a Sin[x], y[0] == 0}, y, x] Out[1]= {{y -> Function[{x}, -(1/2) a E^-x (-1 + E^x Cos[x] - E^x Sin[x])]}} ``` Substitute the solution into an expression: ```wl In[2]:= FullSimplify[y''[x] + y[x] ^ 2 /. %] Out[2]= {(1/4) a E^-2 x (a - 2 E^x (-1 + a Cos[x] - a Sin[x]) + E^2 x (a + 2 Cos[x] - 2 Sin[x] - a Sin[2 x]))} ``` --- Specify only the first argument of ``DSolve`` : ```wl In[1]:= DSolve[y'[x] == y[x]] Out[1]= {{y[x] -> E^x C[1]}} ``` Solve the same equation with an implicit independent variable $\unicode{f817}$ : ```wl In[2]:= DSolve[y' == y] Out[2]= {{y -> Function[{\[FormalX]}, E^\[FormalX] C[1]]}} ``` ### Scope (117) #### Basic Uses (17) Compute the general solution of a first-order differential equation: ```wl In[1]:= DSolve[y'[x] == 3y[x], y[x], x] Out[1]= {{y[x] -> E^3 x C[1]}} ``` Obtain a particular solution by adding an initial condition: ```wl In[2]:= DSolve[{y'[x] == 3y[x], y[0] == 5}, y[x], x] Out[2]= {{y[x] -> 5 E^3 x}} ``` --- Plot the general solution of a differential equation: ```wl In[1]:= sol = DSolve[{y'[x] == -y[x]}, y[x], x] Out[1]= {{y[x] -> E^-x C[1]}} ``` Plot the solution curves for two different values of the arbitrary constant ``C[1]`` : ```wl In[2]:= Plot[{y[x] /. sol /. C[1] -> 2, y[x] /. sol /. C[1] -> 3}, {x, 1, 4}] Out[2]= [image] ``` --- Solve a first-order differential equation using a one-argument specification: ```wl In[1]:= DSolve[y' == 5y] Out[1]= {{y -> Function[{\[FormalX]}, E^5 \[FormalX] C[1]]}} In[2]:= DSolve[y'[x] == 5y[x]] Out[2]= {{y[x] -> E^5 x C[1]}} ``` --- Plot the solution of a boundary value problem: ```wl In[1]:= sol = DSolve[{y'[x] == -3y[x] ^ 2, y[1] == 2}, y[x], x] Out[1]= {{y[x] -> (2/-5 + 6 x)}} In[2]:= Plot[y[x] /. sol, {x, 1, 4}] Out[2]= [image] ``` --- Solve a boundary value problem using a one-argument specification: ```wl In[1]:= DSolve[{y'[x] == x y[x], y[0] == 3}] Out[1]= {{y[x] -> 3 E^(x^2/2)}} ``` --- Verify the solution of a first-order differential equation by using ``y`` in the second argument: ```wl In[1]:= eqns = {y'[x] == y[x] + a Cos[x], y[0] == 3}; In[2]:= sol = DSolve[eqns, y, x] Out[2]= {{y -> Function[{x}, (1/2) (6 E^x + a E^x - a Cos[x] + a Sin[x])]}} In[3]:= eqns /. sol//Simplify Out[3]= {{True, True}} ``` --- Obtain the general solution of a higher-order differential equation: ```wl In[1]:= DSolve[{y''[x] + 4y[x] == 7}, y[x], x] Out[1]= {{y[x] -> (7/4) + C[1] Cos[2 x] + C[2] Sin[2 x]}} ``` Particular solution: ```wl In[2]:= DSolve[{y''[x] + 4y[x] == 7, y[0] == 1, y'[0] == 2}, y[x], x] Out[2]= {{y[x] -> (1/4) (7 - 3 Cos[2 x] + 4 Sin[2 x])}} ``` --- Solve a system of differential equations: ```wl In[1]:= eqns = {Derivative[1][y][x] == y[x] - a z[x], Derivative[1][z][x] == y[x] - z[x], y[0] == 1, z[0] == 4}; In[2]:= sol = DSolve[eqns, {y, z}, x] Out[2]= {{y -> Function[{x}, (E^-Sqrt[1 - a] x (-1 + Sqrt[1 - a] + 4 a + E^2 Sqrt[1 - a] x + Sqrt[1 - a] E^2 Sqrt[1 - a] x - 4 a E^2 Sqrt[1 - a] x)/2 Sqrt[1 - a])], z -> Function[{x}, (E^-Sqrt[1 - a] x (3 + 4 Sqrt[1 - a] - 3 E^2 Sqrt[1 - a] x + 4 Sqrt[1 - a] E^2 Sqrt[1 - a] x)/2 Sqrt[1 - a])]}} ``` Plot the solution: ```wl In[3]:= ParametricPlot[Evaluate[Table[{y[x], z[x]} /. sol, {a, 3, 8}]], {x, 0, 10}] Out[3]= [image] ``` Verify the solution: ```wl In[4]:= eqns /. sol//Simplify Out[4]= {{True, True, True, True}} ``` --- Solve a system of differential-algebraic equations: ```wl In[1]:= sol = DSolve[{y'[x] - 3z[x] == Sin[x], y[x] + z[x] == 1 / 5, y[Pi / 2] == 1 / 2}, {y, z}, x] Out[1]= {{y -> Function[{x}, (1/10) (2 - Cos[x] + 3 Sin[x])], z -> Function[{x}, (1/10) (Cos[x] - 3 Sin[x])]}} In[2]:= Plot[Evaluate[{y[x], z[x], y[x] + z[x]} /. sol], {x, 2, 16}, PlotLegends -> {y[x], z[x], y[x] + z[x]}] Out[2]= [image] ``` --- Solve a system of differential equations using a one-argument specification: ```wl In[1]:= DSolve[{y'[x] == z[x], z'[x] == -y[x]}] Out[1]= {{y[x] -> C[1] Cos[x] + C[2] Sin[x], z[x] -> C[2] Cos[x] - C[1] Sin[x]}} ``` --- Solve a partial differential equation: ```wl In[1]:= eqn = D[u[x, y], x] + 3D[u[x, y], y] + u[x, y] == 1; In[2]:= sol = DSolve[eqn, u, {x, y}] Out[2]= {{u -> Function[{x, y}, E^-x (E^x + C[1][-3 x + y])]}} ``` Obtain a particular solution: ```wl In[3]:= sol /. C[1] -> Sin Out[3]= {{u -> Function[{x, y}, E^-x (E^x + Sin[-3 x + y])]}} ``` Plot the solution: ```wl In[4]:= Plot3D[u[x, y] /. %, {x, -3, 3}, {y, -3, 3}, PlotRange -> All] Out[4]= [image] ``` --- Solve a partial differential equation using a one-argument specification: ```wl In[1]:= DSolve[{D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}]}] Out[1]= {{u[x, t] -> C[1][t - x] + C[2][t + x]}} ``` --- Use different names for the arbitrary constants in the general solution: ```wl In[1]:= DSolve[y''[x] - 4y[x] == 0, y[x], x] Out[1]= {{y[x] -> E^2 x C[1] + E^-2 x C[2]}} In[2]:= DSolve[y''[x] - 4y[x] == 0, y[x], x, GeneratedParameters -> f] Out[2]= {{y[x] -> E^2 x f[1] + E^-2 x f[2]}} ``` --- Solve a delay differential equation: ```wl In[1]:= sol = DSolve[{x'[t] == x[t - 1] ^ 2, x[t /; t ≤ 0] == a t}, x[t], {t, 0, 2}] Out[1]= {{x[t] -> Piecewise[{{a*t, t <= 0}, {(1/3)*(3*a^2*t - 3*a^2*t^2 + a^2*t^3), Inequality[0, Less, t, LessEqual, 1]}, {(1/126)*(42*a^2 - 163*a^4 + 686*a^4*t - 1176*a^4*t^2 + 1064*a^4*t^3 - 553*a^4*t^4 + 168*a^4*t^5 - 28*a^4*t^6 + 2*a^4*t^7), Inequality[1, Less, t, LessEqual, 2]}}, Indeterminate]}} ``` Plot the solution for different values of the parameter: ```wl In[2]:= Plot[Evaluate[Table[x[t] /. sol, {a, 1, 5}]], {t, -1, 2}] Out[2]= [image] ``` --- Solve a hybrid differential equation: ```wl In[1]:= sol = DSolve[{y''[t] == -981 / 100, y[0] == 5, y'[0] == 0, WhenEvent[y[t] == 0, y'[t] -> -(7 / 10)y'[t]]}, y[t], {t, 0, 4}] Out[1]= {{y[t] -> Piecewise[{{5 - (981*t^2)/200, 0 <= t <= (10*Sqrt[10/109])/3}, {(3/200)*(-800 + 34*Sqrt[1090]*t - 327*t^2), Inequality[(10*Sqrt[10/109])/3, Less, t, LessEqual, 8*Sqrt[10/109]]}, {(3*(-13520 + 289*Sqrt[1090]*t - 1635*t^2))/1000, Inequality[8*Sqrt[10/109], Less, t, LessEqual, (169*Sqrt[2/545])/3]}, {(-687154 + 11169*Sqrt[1090]*t - 49050*t^2)/10000, Inequality[(169*Sqrt[2/545])/3, Less, t, LessEqual, 4]}}, Indeterminate]}} ``` This system models a bouncing ball: ```wl In[2]:= Plot[y[t] /. sol, {t, 0, 4}] Out[2]= [image] ``` --- Apply ``N[DSolve[...]]`` to invoke ``NDSolve`` if symbolic solution fails: ```wl In[1]:= DSolve[{y''[x] + x ^ 2y'[x] + Sin[y[x]] == 1, y[0] == 1, y'[0] == 1}, y, {x, 0, 1}] Out[1]= DSolve[{Sin[y[x]] + x^2 Derivative[1][y][x] + Derivative[2][y][x] == 1, y[0] == 1, Derivative[1][y][0] == 1}, y, {x, 0, 1}] In[2]:= N[%] Out[2]= {{y -> InterpolatingFunction[{{0., 1.}}, {5, 7, 2, {30}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False}, {{0., 0.0000966430845726476, 0.0001932861691452952, 0.004773198829667701, 0.009353111490190108, 0.013933024150712514, 0.028261452402542964, 0 ... 06, 0.9168218629734495, 0.9584109314867247, 1.}}, {Developer`PackedArrayForm, {0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42, 45, 48, 51, 54, 57, 60, 63, 66, 69, 72, 75, 78, 81, 84, 87, 90}, CompressedData["«1055»"]}, {Automatic}]}} In[3]:= Table[y[x] /. %[[1]], {x, 0, 1, 0.2}] Out[3]= {1., 1.20237, 1.40565, 1.60295, 1.78648, 1.94882} ``` --- Find the general solution of an ODE with quantities: ```wl In[1]:= DSolve[Derivative[2][x][t] == Quantity[5, ("Meters"/"Seconds"^2)], x[t], t] Out[1]= {{x[t] -> C[1] + C[2] t + (Quantity[5/2, "Meters"/"Seconds"^2]) t^2}} ``` Solve an initial-value problem: ```wl In[2]:= DSolve[{Derivative[2][x][t] == Quantity[5, ("Meters"/"Seconds"^2)], x[Quantity[0, "Seconds"]] == Quantity[3, "Meters"], Derivative[1][x][Quantity[0, "Seconds"]] == Quantity[7, ("Meters"/"Seconds")]}, x[t], t] Out[2]= {{x[t] -> Quantity[3, "Meters"] + (Quantity[7, "Meters"/"Seconds"]) t + (Quantity[5/2, "Meters"/"Seconds"^2]) t^2}} ``` #### Linear Differential Equations (6) Exponential equation: ```wl In[1]:= DSolve[{y'[x] == ay[x], y[0] == 1}, y[x], x] Out[1]= {{y[x] -> E^a x}} ``` Inhomogeneous first-order equation: ```wl In[2]:= DSolve[y'[x] + x y[x] == Cos[x], y[x], x] Out[2]= {{y[x] -> E^-(x^2/2) C[1] + (1/2) E^(1/2) - (x^2/2) Sqrt[(π/2)] (Erfi[(-I + x/Sqrt[2])] + Erfi[(I + x/Sqrt[2])])}} ``` --- Solve a boundary value problem: ```wl In[1]:= DSolve[{y''[x] - x y[x] == 0, y[0] == 1, y[9] == 1}, y, x] Out[1]= {{y -> Function[{x}, (Sqrt[3] AiryAi[x] - AiryBi[x] - 3^2 / 3 AiryAi[x] AiryBi[9] Gamma[(2/3)] + 3^2 / 3 AiryAi[9] AiryBi[x] Gamma[(2/3)]) / (Sqrt[3] AiryAi[9] - AiryBi[9])]}} ``` Plot the solution: ```wl In[2]:= Plot[y[x] /. %, {x, -10, 10}] Out[2]= [image] ``` --- Second-order equation with constant coefficients: ```wl In[1]:= DSolve[y''[x] + 4y'[x] + 5y[x] == 0, y[x], x] Out[1]= {{y[x] -> E^-2 x C[2] Cos[x] + E^-2 x C[1] Sin[x]}} ``` Cauchy–Euler equation: ```wl In[2]:= DSolve[x ^ 2y''[x] + 4x y'[x] + 7y[x] == 0, y[x], x] Out[2]= {{y[x] -> (C[2] Cos[(1/2) Sqrt[19] Log[x]]/x^3 / 2) + (C[1] Sin[(1/2) Sqrt[19] Log[x]]/x^3 / 2)}} ``` Second-order equation with variable coefficients, solved in terms of elementary functions: ```wl In[3]:= DSolve[x y''[x] + 2y'[x] - x y[x] == Sin[x], y, x] Out[3]= {{y -> Function[{x}, (E^-x C[1]/x) + (E^x C[2]/2 x) - (Sin[x]/2 x)]}} ``` Airy's equation: ```wl In[4]:= DSolve[y''[x] - x y[x] == 0, y, x] Out[4]= {{y -> Function[{x}, AiryAi[x] C[1] + AiryBi[x] C[2]]}} ``` Spheroidal equation: ```wl In[5]:= DSolve[(1 - x ^ 2)y''[x] - 2x y'[x] + (49 - 49x ^ 2 - 16 / (1 - x ^ 2) + SpheroidalEigenvalue[5, 4, 7])y[x] == 0, y, x] Out[5]= {{y -> Function[{x}, C[1] SpheroidalPS[5, 4, 7, x] + C[2] SpheroidalQS[5, 4, 7, x]]}} ``` Equations with nonrational coefficients: ```wl In[6]:= DSolve[y''[x] + (8 + 3Cos[x] ^ 2)y[x] == 0, y, x] Out[6]= {{y -> Function[{x}, C[1] MathieuC[(19/2), -(3/4), x] + C[2] MathieuS[(19/2), -(3/4), x]]}} In[7]:= DSolve[y''[x] + (9 + Sech[x] ^ 2)y[x] == 0, y, x] Out[7]= {{y -> Function[{x}, C[1] LegendreP[(1/2) (-1 + Sqrt[5]), 3 I, Tanh[x]] + C[2] LegendreQ[(1/2) (-1 + Sqrt[5]), 3 I, Tanh[x]]]}} In[8]:= DSolve[y''[x] == Exp[-x]y[x], y, x] Out[8]= {{y -> Function[{x}, BesselI[0, 2 Sqrt[E^-x]] C[1] + 2 BesselK[0, 2 Sqrt[E^-x]] C[2]]}} ``` --- Higher-order equations: ```wl In[1]:= DSolve[y'''[x] + 4y'[x] == 5y[x], y, x] Out[1]= {{y -> Function[{x}, E^x C[3] + E^-x / 2 C[2] Cos[(Sqrt[19] x/2)] + E^-x / 2 C[1] Sin[(Sqrt[19] x/2)]]}} ``` Solution in terms of hypergeometric functions: ```wl In[2]:= DSolve[3x ^ 3y''''[x] + 19x ^ 2y'''[x] + (2x - 3x ^ 2)y''[x] + (2 - 9x)y'[x] - 3y[x] == 0, y, x] Out[2]= {{y -> Function[{x}, -(-1)^2 / 3 x^5 / 3 C[2] HypergeometricPFQ[{(8/3)}, {(11/3) - Sqrt[5], (11/3) + Sqrt[5]}, x] + (-1)^-1 - Sqrt[5] x^-1 - Sqrt[5] C[3] HypergeometricPFQ[{-Sqrt[5]}, {1 - 2 Sqrt[5], -(5/3) - Sqrt[5]}, x] + (-1)^-1 + Sqrt[5] x^-1 + Sqrt[5] C[4] HypergeometricPFQ[{Sqrt[5]}, {-(5/3) + Sqrt[5], 1 + 2 Sqrt[5]}, x] + C[1] HypergeometricPFQ[{1, 1}, {-(2/3), 2 - Sqrt[5], 2 + Sqrt[5]}, x]]}} ``` Fourth-order equation solved in terms of Kelvin functions: ```wl In[3]:= DSolve[(525 + x ^ 4)y[x] - 51(-x y'[x] + x ^ 2y''[x]) + 2x ^ 3y'''[x] + x ^ 4y''''[x] == 0, y, x] Out[3]= {{y -> Function[{x}, C[4] KelvinBei[5, x] + C[3] KelvinBer[5, x] + C[2] KelvinKei[5, x] + C[1] KelvinKer[5, x]]}} ``` --- Linear equation with ``q``-rational coefficients: ```wl In[1]:= DSolve[((-1 + 2^x) Derivative[4][y][x]/1 + 2^x) + 2^x y[x] == 0, y, x] Out[1]= {{y -> Function[{x}, C[4] DifferentialRoot[Function[{\[FormalY], \[FormalX]}, {(1 + \[FormalX])*\[FormalY][\[FormalX]] + ((-1 + \[FormalX])*Log[2]^4)*Derivative[1][\[FormalY]][\[FormalX]] + (((7*(-1 + \[FormalX]))*\[FormalX])*Log[2]^4)*Derivative[2][\[FormalY]][\[FormalX]] + (((6*(-1 + \[FormalX]))*\[FormalX]^2)*Log[2]^4)* Derivative[3][\[FormalY]][\[FormalX]] + (((-1 + \[FormalX])*\[FormalX] ... \[FormalX]))*\[FormalX]^2)*Log[2]^4)* Derivative[3][\[FormalY]][\[FormalX]] + (((-1 + \[FormalX])*\[FormalX]^3)*Log[2]^4)*Derivative[4][\[FormalY]][\[FormalX]] == 0, \[FormalY][2] == 1, Derivative[1][\[FormalY]][2] == 0, Derivative[2][\[FormalY]][2] == 0, Derivative[3][\[FormalY]][2] == 0}], Function[\[FormalX], {{Re[\[FormalX]] <= 0, Im[\[FormalX]] == 0}}]][2^x]]}} ``` --- Higher-order inhomogeneous equation with constant coefficients: ```wl In[1]:= DSolve[Derivative[4][y][x] + 5y'''[x] + 2y[x] == x E^3x Cos[3x], y[x], x] Out[1]= {{y[x] -> E^x Root[2 + 5*#1^3 + #1^4 & , 1, 0] C[1] + E^x Root[2 + 5*#1^3 + #1^4 & , 2, 0] C[2] + E^x Root[2 + 5*#1^3 + #1^4 & , 3, 0] C[3] + E^x Root[2 + 5*#1^3 + #1^4 & , 4, 0] C[4] - (1/22404634562)E^3 x (-26914788 Cos[3 x] + 31328936 x Cos[3 x] - 8230653 Sin[3 x] - 14288535 x Sin[3 x])}} ``` #### Nonlinear Differential Equations (6) Solve a Riccati equation: ```wl In[1]:= DSolve[y'[x] + x y'[x] ^ 2 == 1, y, x] Out[1]= {{y -> Function[{x}, C[1] + (1/2) (2 Sqrt[1 + 4 x] - 2 ArcTanh[Sqrt[1 + 4 x]] - Log[x])]}, {y -> Function[{x}, C[1] + (1/2) (-2 Sqrt[1 + 4 x] + 2 ArcTanh[Sqrt[1 + 4 x]] - Log[x])]}} ``` --- Implicit solution for an Abel equation: ```wl In[1]:= DSolve[y'[x] == y[x] ^ 3 - ((x + 1)y[x] ^ 2) / x, y[x], x]//Quiet Out[1]= Solve[(E^-x + (1/y[x])/x) + C[1] + ExpIntegralEi[-x + (1/y[x])] == 0, y[x]] ``` --- Homogeneous equation: ```wl In[1]:= DSolve[y'[x] - Sqrt[(y[x] / x)] == y[x] / x, y, x] Out[1]= {{y -> Function[{x}, (1/4) (x C[1]^2 + 2 x C[1] Log[x] + x Log[x]^2)]}} ``` --- Solution in terms of ``WeierstrassP`` : ```wl In[1]:= DSolve[y''[x] == y[x] ^ 2 + 1, y, x] Out[1]= {{y -> Function[{x}, 6^1 / 3 WeierstrassP[(x + C[1]/6^1 / 3), {-2 6^1 / 3, C[2]}]]}} ``` --- Solution in terms of hyperbolic functions: ```wl In[1]:= DSolve[y''[x] == 3y[x] y'[x] + (3y[x] ^ 2 + 4y[x] + 1), y[x], x][[1]]//Simplify Out[1]= {y[x] -> ((I Sqrt[6] Sqrt[E^-2 x] Sqrt[C[1]] - 6 C[2]) Cosh[Sqrt[(3/2)] Sqrt[E^-2 x] Sqrt[C[1]]] + (-3 I + 2 Sqrt[6] Sqrt[E^-2 x] Sqrt[C[1]] C[2]) Sinh[Sqrt[(3/2)] Sqrt[E^-2 x] Sqrt[C[1]]]) / (6 C[2] Cosh[Sqrt[(3/2)] Sqrt[E^-2 x] Sqrt[C[1]]] + 3 I Sinh[Sqrt[(3/2)] Sqrt[E^-2 x] Sqrt[C[1]]])} ``` --- Solve a nonlinear boundary value problem: ```wl In[1]:= DSolve[{y''[x] == y[x] ^ 3, y[0] == 5, y'[0] == 25 / Sqrt[2]}, y[x], x] Out[1]= {{y[x] -> -(10/-2 + 5 Sqrt[2] x)}} ``` #### Systems of Differential Equations (8) Diagonal linear system: ```wl In[1]:= DSolve[{y'[x] == x ^ 2y[x], z'[x] == 5z[x]}, {y, z}, x] Out[1]= {{y -> Function[{x}, E^(x^3/3) C[1]], z -> Function[{x}, E^5 x C[2]]}} ``` --- Inhomogeneous linear system with constant coefficients: ```wl In[1]:= With[{θ = Pi / 3}, DSolve[{y'[x] == Cos[θ]y[x] - Sin[θ]z[x] + 1, z'[x] == Sin[θ]y[x] + Cos[θ]z[x], y[0] == 1, z[0] == 1}, {y, z}, x]] Out[1]= {{y -> Function[{x}, (1/2) (3 E^x / 2 Cos[(Sqrt[3] x/2)] - Cos[(Sqrt[3] x/2)]^2 - 2 E^x / 2 Sin[(Sqrt[3] x/2)] + Sqrt[3] E^x / 2 Sin[(Sqrt[3] x/2)] - Sin[(Sqrt[3] x/2)]^2)], z -> Function[{x}, (1/2) (2 E^x / 2 Cos[(Sqrt[3] x/2)] - Sqrt[3] E^x / 2 Cos[(Sqrt[3] x/2)] + Sqrt[3] Cos[(Sqrt[3] x/2)]^2 + 3 E^x / 2 Sin[(Sqrt[3] x/2)] + Sqrt[3] Sin[(Sqrt[3] x/2)]^2)]}} In[2]:= Plot[Evaluate[{y[x], z[x]} /. First[%]], {x, 0, 9}] Out[2]= [image] ``` --- Linear system with rational function coefficients: ```wl In[1]:= DSolve[{Derivative[1][y][x] == -(4 y[x]/x) + (4 z[x]/x), Derivative[1][z][x] == ((4/x) - (x/4)) y[x] - (4 z[x]/x)}, {y[x], z[x]}, x] Out[1]= {{y[x] -> (BesselJ[4, x] C[1]/x^4) + (BesselY[4, x] C[2]/x^4), z[x] -> ((BesselJ[3, x]/8 x^3) - (BesselJ[5, x]/8 x^3)) C[1] + ((BesselY[3, x]/8 x^3) - (BesselY[5, x]/8 x^3)) C[2]}} ``` --- Nonlinear system: ```wl In[1]:= DSolve[ {y'[x] == Exp[z[x]] + 1, z'[x] == y[x] - x}, {y, z}, x]//Quiet Out[1]= {{z -> Function[{x}, Log[C[1] + C[1] Tan[(1/2) (Sqrt[2] x Sqrt[C[1]] + 2 Sqrt[2] Sqrt[C[1]] C[2])]^2]], y -> Function[{x}, x + Sqrt[2] Sqrt[C[1]] Tan[(1/2) (Sqrt[2] x Sqrt[C[1]] + 2 Sqrt[2] Sqrt[C[1]] C[2])]]}} ``` --- Solve a linear system of ODEs using vector variables: ```wl In[1]:= m = {{4, -6}, {1, -1}}; In[2]:= DSolve[x'[t] == m.x[t], x[t]∈Vectors[2], t] Out[2]= {{x[t] -> {E^t (-2 + 3 E^t) C[1] - 6 E^t (-1 + E^t) C[2], E^t (-1 + E^t) C[1] - E^t (-3 + 2 E^t) C[2]}}} ``` Alternatively, define $x$ as a ``VectorSymbol``: ```wl In[3]:= x = VectorSymbol["x", 2] Out[3]= VectorSymbol["x", 2] In[4]:= DSolve[x'[t] == m.x[t], x[t], t] Out[4]= {{VectorSymbol["x", 2][t] -> {E^t (-2 + 3 E^t) C[1] - 6 E^t (-1 + E^t) C[2], E^t (-1 + E^t) C[1] - E^t (-3 + 2 E^t) C[2]}}} ``` --- Solve a linear system of four ODEs using matrix variables: ```wl In[1]:= m = {{0, 1}, {-1, 0}}; x0 = {{1, 2}, {3, 4}}; In[2]:= sol = DSolve[{x'[t] == m.x[t], x[0] == x0}, Element[x[t], Matrices[{2, 2}]], t]; x[t] /. sol[[1]]//MatrixForm Out[2]//MatrixForm= (⁠| | | | ----------------- | --------------------- | | Cos[t] + 3 Sin[t] | 2 (Cos[t] + 2 Sin[t]) | | 3 Cos[t] - Sin[t] | 2 (2 Cos[t] - Sin[t]) |⁠) ``` Alternatively, define $x$ as a ``MatrixSymbol`` : ```wl In[3]:= x = MatrixSymbol["x", {2, 2}] Out[3]= MatrixSymbol["x", {2, 2}] In[4]:= sol = DSolve[{x'[t] == m.x[t], x[0] == x0}, x[t], t]; x[t] /. sol[[1]]//MatrixForm Out[4]//MatrixForm= (⁠| | | | ----------------- | --------------------- | | Cos[t] + 3 Sin[t] | 2 (Cos[t] + 2 Sin[t]) | | 3 Cos[t] - Sin[t] | 2 (2 Cos[t] - Sin[t]) |⁠) ``` Plot the solution: ```wl In[5]:= Plot[Evaluate[Flatten[x[t] /. sol[[1]]]], {t, 0, 10}] Out[5]= [image] ``` --- Solve inhomogeneous linear system of ODEs with constant coefficients: ```wl In[1]:= m = {{2, -1}, {1, -2}}; In[2]:= g = {(2t - 5)Cos[3t], (3t + 2)Sin[3t]} E^-2t; In[3]:= DSolve[x'[t] - m.x[t] == g, Element[x[t], Vectors[2]], t] Out[3]= {{x[t] -> {(1/6) E^-Sqrt[3] t (3 - 2 Sqrt[3] + 3 E^2 Sqrt[3] t + 2 Sqrt[3] E^2 Sqrt[3] t) C[1] - (E^-Sqrt[3] t (-1 + E^2 Sqrt[3] t) C[2]/2 Sqrt[3]) + E^-2 t ((991/1352) - (27 t/52)) Cos[3 t] + E^-2 t (-(209/676) + (9 t/26)) Sin[3 t], (E^-Sqrt[3] t (-1 + E^2 Sqrt[3] t) C[1]/2 Sqrt[3]) - (1/6) E^-Sqrt[3] t (-3 - 2 Sqrt[3] - 3 E^2 Sqrt[3] t + 2 Sqrt[3] E^2 Sqrt[3] t) C[2] + E^-2 t (-(105/169) - (29 t/26)) Cos[3 t] + E^-2 t ((833/1352) - (9 t/52)) Sin[3 t]}}} ``` --- Solve a coupled system of linear and nonlinear ODEs: ```wl In[1]:= DSolve[{Derivative[1][y][x] == -(4 y[x]/x) + (4 z[x]/x), Derivative[1][z][x] == ((4/x) - (x/4)) y[x] - (4 z[x]/x), Derivative[1][y1][x] == 1 + E^z1[x], Derivative[1][z1][x] == -x + y1[x], w1'[x] + 2 E ^ w2[x] == 0, w2'[x] == -w1[x]}, {y[x], z[x], y1[x], z1[x], w1[x], w2[x]}, x]//Quiet Out[1]= {{y[x] -> (BesselJ[4, x] C[1]/x^4) + (BesselY[4, x] C[2]/x^4), z[x] -> ((BesselJ[3, x] - BesselJ[5, x]) C[1]/8 x^3) + ((BesselY[3, x] - BesselY[5, x]) C[2]/8 x^3), z1[x] -> Log[C[3] + C[3] Tan[(1/2) (Sqrt[2] x Sqrt[C[3]] + 2 Sqrt[2] C[2] Sqrt[C[3]])]^2], y1[x] -> x + Sqrt[2] Sqrt[C[3]] Tan[(1/2) (Sqrt[2] x Sqrt[C[3]] + 2 Sqrt[2] C[2] Sqrt[C[3]])], w2[x] -> Log[C[5] + C[5] Tan[x Sqrt[C[5]] - 2 C[2] Sqrt[C[5]]]^2], w1[x] -> -2 Sqrt[C[5]] Tan[x Sqrt[C[5]] - 2 C[2] Sqrt[C[5]]]}} ``` #### Differential-Algebraic Equations (3) Solve a system of linear differential-algebraic equations: ```wl In[1]:= DSolve[{2 y'[x] + z'[x] == 4y[x] + x, y[x] - z[x] == 1}, {y[x], z[x]}, x] Out[1]= {{y[x] -> (1/3) + 6 ((1/288) (-25 - 12 x) + E^4 x / 3 C[1]), z[x] -> -(2/3) + 6 ((1/288) (-25 - 12 x) + E^4 x / 3 C[1])}} ``` --- Solve a boundary value problem: ```wl In[1]:= DSolve[{y'[x] - 4z[x] == Cos[x], y[x] + z[x] == 1 / 2, y[Pi / 2] == 1 / 2}, {y, z}, x] Out[1]= {{y -> Function[{x}, (1/34) E^-4 x (-2 E^2 π + 17 E^4 x + 8 E^4 x Cos[x] + 2 E^4 x Sin[x])], z -> Function[{x}, -(1/17) E^-4 x (-E^2 π + 4 E^4 x Cos[x] + E^4 x Sin[x])]}} In[2]:= Plot[Evaluate[{y[x], z[x], y[x] + z[x]} /. %], {x, 2, 16}] Out[2]= [image] ``` --- An index-2 differential-algebraic equation: ```wl In[1]:= DSolve[{x'[t] == x[t] + 2y[t], y[t] == x[t] + 2 ^ t, x[t] + z[t] == 0}, {x[t], y[t], z[t]}, t] Out[1]= {{x[t] -> (1/8) (-E^3 t C[1] + (16 E^3 t + t (-3 + Log[2])/-3 + Log[2])), y[t] -> 2^t + (1/8) (-E^3 t C[1] + (16 E^3 t + t (-3 + Log[2])/-3 + Log[2])), z[t] -> (1/8) (E^3 t C[1] - (16 E^3 t + t (-3 + Log[2])/-3 + Log[2]))}} ``` #### Delay Differential Equations (2) Solve a delay differential equation: ```wl In[1]:= sol = DSolve[{x'[t] == x[t - 1] ^ 2, x[t /; t ≤ 0] == a t}, x[t], {t, 0, 2}] Out[1]= {{x[t] -> Piecewise[{{a*t, t <= 0}, {(1/3)*(3*a^2*t - 3*a^2*t^2 + a^2*t^3), Inequality[0, Less, t, LessEqual, 1]}, {(1/126)*(42*a^2 - 163*a^4 + 686*a^4*t - 1176*a^4*t^2 + 1064*a^4*t^3 - 553*a^4*t^4 + 168*a^4*t^5 - 28*a^4*t^6 + 2*a^4*t^7), Inequality[1, Less, t, LessEqual, 2]}}, Indeterminate]}} ``` Plot the solution for different values of the parameter ``a`` : ```wl In[2]:= Plot[Evaluate[Table[x[t] /. sol[[1]], {a, 1, 5}]], {t, -1, 2}] Out[2]= [image] ``` --- Solve a system of delay differential equations: ```wl In[1]:= eqns = {x'[t] == a y[t - 1] + y[t - 3], y'[t] == x[t - 1], x[t /; t ≤ 0] == t, y[t /; t ≤ 0] == t ^ 2}; In[2]:= sol = DSolve[eqns, {x[t], y[t]}, {t, 0, 5}] Out[2]= {{x[t] -> Piecewise[{{t, t <= 0}, {(1/3)*(27*t + 3*a*t - 9*t^2 - 3*a*t^2 + t^3 + a*t^3), Inequality[0, Less, t, LessEqual, 1]}, {(1/6)*(-2*a + 54*t + 9*a*t - 18*t^2 - 6*a*t^2 + 2*t^3 + a*t^3), Inequality[1, Less, t, LessEqual, 2]}, {( ... *t - 12552*a^2*t + 4320*t^2 + 26550*a*t^2 + 7290*a^2*t^2 - 600*t^3 - 7560*a*t^3 - 2200*a^2*t^3 + 30*t^4 + 1080*a*t^4 + 360*a^2*t^4 - 72*a*t^5 - 30*a^2*t^5 + 2*a*t^6 + a^2*t^6), Inequality[4, Less, t, LessEqual, 5]}}, Indeterminate]}} ``` Plot the solution for different values of the parameter ``a`` : ```wl In[3]:= ParametricPlot[Evaluate[Table[{x[t], y[t]} /. sol, {a, -1, 2, 1 / 3}]], {t, -2, 5}, Exclusions -> None] Out[3]= [image] ``` #### Piecewise Differential Equations (4) Using a piecewise forcing function: ```wl In[1]:= DSolve[{y''[x] - y[x] == Max[x, x ^ 2], y[0] == 1, y'[0] == 2}, y[x], x]//Simplify Out[1]= {{y[x] -> Piecewise[{{(1/2)*(E^(-x) + 5*E^x - 2*(2 + x^2)), x <= 0}, {-E^(-x) + 2*E^x - x, Inequality[0, Less, x, LessEqual, 1]}}, (1/2)*(E^(1 - x) + 3*E^(-1 + x) - 2/E^x + 4*E^x - 2*(2 + x^2))]}} ``` A differential equation with a piecewise coefficient: ```wl In[2]:= DSolve[y'[x] + Clip[x] ^ 2y[x] == 0, y[x], x] Out[2]= {{y[x] -> E^Piecewise[{{-x, x <= -1}, {2/3 - x^3/3, Inequality[-1, Less, x, LessEqual, 1]}}, 4/3 - x] C[1]}} ``` --- A nonlinear piecewise-defined differential equation: ```wl In[1]:= sol = DSolve[{y'[x] == Piecewise[{{Cos[x] ^ 2, x > 2}}, y[x] ^ 2 / 3], y[0] == 1}, y, x] Out[1]= {{y -> Function[{x}, Piecewise[{{-(3/(-3 + x)), x <= 2}}, (1/4)*(8 + 2*x - Sin[4] + Sin[2*x])]]}} In[2]:= Plot[y[x] /. sol, {x, 0, 7}] Out[2]= [image] ``` --- Differential equations involving generalized functions: ```wl In[1]:= DSolve[y''[x] - y[x] == HeavisideTheta[x], y[x], x] Out[1]= {{y[x] -> E^x C[1] + E^-x C[2] + (1/2) E^-x (-1 + E^x)^2 HeavisideTheta[x]}} ``` A simple impulse response or Green's function: ```wl In[2]:= DSolve[y'[x] + 7y[x] == DiracDelta[x], y[x], x] Out[2]= {{y[x] -> E^-7 x C[1] + E^-7 x HeavisideTheta[x]}} ``` --- Solve a piecewise differential equation on the interval ``{0, 1}`` : ```wl In[1]:= DSolve[y'[x] == UnitStep[x], y, {x, 0, 1}] Out[1]= {{y -> Function[{x}, x + C[1]]}} ``` Solutions on other intervals: ```wl In[2]:= DSolve[y'[x] == UnitStep[x], y, {x, -1, 0}] Out[2]= {{y -> Function[{x}, C[1]]}} In[3]:= DSolve[y'[x] == UnitStep[x], y, {x, -1, 1}] Out[3]= {{y -> Function[{x}, C[1] + x UnitStep[x]]}} ``` #### Hybrid Differential Equations (8) Solve a first-order differential equation with a time-dependent event: ```wl In[1]:= sol = DSolve[{x'[t] == x[t], x[0] == 1, WhenEvent[Mod[t, 1] == 0, x[t] -> x[t] + 10]}, x, {t, 0, 3}] Out[1]= {{x -> Function[{t}, Piecewise[{{E^t, 0 <= t <= 1}, {E^(-1 + t)*(10 + E), Inequality[1, Less, t, LessEqual, 2]}, {E^(-2 + t)*(10 + 10*E + E^2), Inequality[2, Less, t, LessEqual, 3]}}, Indeterminate]]}} ``` Plot the solution: ```wl In[2]:= Plot[x[t] /. sol[[1]], {t, 0, 4}, Exclusions -> None] Out[2]= [image] ``` --- Solve a second-order differential equation with a time-dependent event: ```wl In[1]:= sol = DSolve[{x''[t] == x[t] + t, x[0] == 1, x'[0] == 3, WhenEvent[Mod[2t, 1] == 0, {x[t] -> x[t] + t, x'[t] -> -x'[t]}]}, x[t], {t, 0, 4}] Out[1]= {{x[t] -> Piecewise[{{-(3/(E^t*2)) + (5*E^t)/2 - t, 0 <= t <= 1/2}, {(1/4)*(-3*E^(1/2 - t) + 10*E^(1 - t) - 6*E^(-1 + t) + 5*E^(-(1/2) + t) - 4*t), Inequality[1/2, Less, t, LessEqual, 1]}, {-2*E^(1 - t) + (5/4)*E^(3/2 - t) - (3/4)*E^(-(3/ ... *(-6*E^(3 - t) + 21*E^(7/2 - t) - 3*E^(-(7/2) + t) + 34*E^(-3 + t) - 4*t), Inequality[3, Less, t, LessEqual, 7/2]}, {(17*E^(4 - t))/2 - (3*E^(-4 + t))/2 + 8*E^(-(7/2) + t) - t, Inequality[7/2, Less, t, LessEqual, 4]}}, Indeterminate]}} ``` Plot the solution: ```wl In[2]:= Plot[x[t] /. sol[[1]], {t, 0, 4}, Exclusions -> None] Out[2]= [image] ``` --- Solve a system of differential equations with a state-dependent event: ```wl In[1]:= eqns = {y'[t] == z[t], z'[t] == -10, y[0] == 1, z[0] == 0}; events = {WhenEvent[y[t] == 0, {y[t] -> 0, z[t] -> -(70 / 100)z[t]}]}; sol = DSolve[Join[eqns, events], {y[t], z[t]}, {t, 0, 2}] Out[1]= {{y[t] -> Piecewise[{{1 - 5*t^2, 0 <= t <= 1/Sqrt[5]}, {-(12/5) + (17*t)/Sqrt[5] - 5*t^2, Inequality[1/Sqrt[5], Less, t, LessEqual, 12/(5*Sqrt[5])]}, {-(1014/125) + (289*t)/(10*Sqrt[5]) - 5*t^2, Inequality[12/(5*Sqrt[5]), Less, t, LessEqual, ... LessEqual, 169/(50*Sqrt[5])]}, {3723/(100*Sqrt[5]) - 10*t, Inequality[169/(50*Sqrt[5]), Less, t, LessEqual, 2033/(500*Sqrt[5])]}, {43061/(1000*Sqrt[5]) - 10*t, Inequality[2033/(500*Sqrt[5]), Less, t, LessEqual, 2]}}, Indeterminate]}} ``` Plot the solution for ``y`` : ```wl In[2]:= Plot[y[t] /. sol[[1]], {t, 0, 2}] Out[2]= [image] ``` --- Stop the integration when an event occurs: ```wl In[1]:= DSolve[{x'[t] == x[t], x[0] == 1, WhenEvent[t == 2, "StopIntegration"]}, x, {t, 0, 3}] Out[1]= {{x -> Function[{t}, Piecewise[{{E^t, 0 <= t <= 2}}, Indeterminate]]}} ``` --- Remove an event after it has occurred once: ```wl In[1]:= sol = DSolve[{x''[t] + x[t] == 0, x[0] == 0, x'[0] == 1, WhenEvent[x[t] == 0, {x'[t] -> -x'[t], "RemoveEvent"}]}, x, {t, 0, 10}] Out[1]= {{x -> Function[{t}, Piecewise[{{Sin[t], 0 <= t <= Pi}, {-Sin[t], Inequality[Pi, Less, t, LessEqual, 10]}}, Indeterminate]]}} In[2]:= Plot[x[t] /. sol, {t, 0, 10}] Out[2]= [image] ``` --- Specify that a variable maintains its value between events: ```wl In[1]:= sol = DSolve[{x'[t] == a[t], x[0] == 0, a[0] == 1, WhenEvent[Mod[x[t], 1] == 0, a[t] -> -a[t]]}, {x, a}, {t, 0, 4}, DiscreteVariables -> a] Out[1]= {{x -> Function[{t}, Piecewise[{{t, 0 <= t <= 1}, {2 - t, Inequality[1, Less, t, LessEqual, 2]}, {-2 + t, Inequality[2, Less, t, LessEqual, 3]}, {4 - t, Inequality[3, Less, t, LessEqual, 4]}}, Indeterminate]], a -> Function[{t}, Piecewise[{{1, 0 <= t <= 1}, {-1, Inequality[1, Less, t, LessEqual, 2]}, {1, Inequality[2, Less, t, LessEqual, 3]}, {-1, Inequality[3, Less, t, LessEqual, 4]}}, Indeterminate]]}} In[2]:= Plot[x[t] /. sol, {t, 0, 4}] Out[2]= [image] ``` --- Prescribe a different action at each event: ```wl In[1]:= sol = DSolve[{x'[t] == y[t], y'[t] == 1, x[0] == 0, y[0] == 0, WhenEvent[t == 1, x[t] -> x[t] + 1], WhenEvent[t == 2, y[t] -> y[t] + 2]}, {x, y}, {t, 0, 3}] Out[1]= {{x -> Function[{t}, Piecewise[{{t^2/2, 0 <= t <= 1}, {(1/2)*(2 + t^2), Inequality[1, Less, t, LessEqual, 2]}, {(1/2)*(-6 + 4*t + t^2), Inequality[2, Less, t, LessEqual, 3]}}, Indeterminate]], y -> Function[{t}, Piecewise[{{t, 0 <= t <= 1 || Inequality[1, Less, t, LessEqual, 2]}, {2 + t, Inequality[2, Less, t, LessEqual, 3]}}, Indeterminate]]}} In[2]:= ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 3}] Out[2]= [image] ``` --- Several events with the same action: ```wl In[1]:= sol = DSolve[{x'[t] == 1, x[0] == 0, WhenEvent[{t == 1, t == 2}, x[t] -> x[t] + 1]}, x, {t, 0, 3}] Out[1]= {{x -> Function[{t}, Piecewise[{{t, 0 <= t <= 1}, {1 + t, Inequality[1, Less, t, LessEqual, 2]}, {2 + t, Inequality[2, Less, t, LessEqual, 3]}}, Indeterminate]]}} In[2]:= Plot[x[t] /. sol, {t, 0, 3}] Out[2]= [image] ``` #### Sturm-Liouville Problems (6) Solve an eigenvalue problem with Dirichlet conditions: ```wl In[1]:= sol = DSolve[{y''[x] + λ y[x] == 0, y[0] == 0, y[π] == 0}, y[x], x] Out[1]= {{y[x] -> Piecewise[{{C[1]*Sin[x*Sqrt[λ]], Element[\[FormalN], Integers] && \[FormalN] >= 1 && λ == \[FormalN]^2}}, 0]}} ``` Make a table of the first 5 eigenfunctions: ```wl In[2]:= eigfuns = Table[y[x] /. sol[[1]] //. {\[FormalN] -> i, λ -> \[FormalN]^2} /. {C[1] -> 1}, {i, 5}] Out[2]= {Sin[x], Sin[2 x], Sin[3 x], Sin[4 x], Sin[5 x]} ``` Plot the eigenfunctions: ```wl In[3]:= Plot[Evaluate[eigfuns], {x, 0, Pi}] Out[3]= [image] ``` --- Solve an eigenvalue problem with Neumann conditions: ```wl In[1]:= sol = DSolve[{y''[x] + λ y[x] == 0, y'[0] == 0, y'[π] == 0}, y[x], x] Out[1]= {{y[x] -> Piecewise[{{C[1]*Cos[x*Sqrt[λ]], Element[\[FormalN], Integers] && \[FormalN] >= 0 && λ == \[FormalN]^2}}, 0]}} ``` Make a table of the first 5 eigenfunctions: ```wl In[2]:= eigfuns = Table[y[x] /. sol[[1]] //. {\[FormalN] -> i, λ -> \[FormalN]^2} /. {C[1] -> 1}, {i, 0, 4}] Out[2]= {1, Cos[x], Cos[2 x], Cos[3 x], Cos[4 x]} ``` Plot the eigenfunctions: ```wl In[3]:= Plot[Evaluate[eigfuns], {x, 0, Pi}] Out[3]= [image] ``` --- Solve an eigenvalue problem with mixed boundary conditions: ```wl In[1]:= sol = DSolve[{y''[x] + λ y[x] == 0, y[0] == 0, y'[π] == 0}, y[x], x] Out[1]= {{y[x] -> Piecewise[{{C[1]*Sin[x*Sqrt[λ]], Element[\[FormalN], Integers] && \[FormalN] >= 1 && λ == (-(1/2) + \[FormalN])^2}}, 0]}} ``` Make a table of the first 5 eigenfunctions: ```wl In[2]:= eigfuns = Table[y[x] /. sol[[1]] //. {\[FormalN] -> i, λ -> (-(1/2) + \[FormalN])^2} /. {C[1] -> 1}, {i, 5}] Out[2]= {Sin[(x/2)], Sin[(3 x/2)], Sin[(5 x/2)], Sin[(7 x/2)], Sin[(9 x/2)]} ``` Plot the eigenfunctions: ```wl In[3]:= Plot[Evaluate[eigfuns], {x, 0, Pi}] Out[3]= [image] ``` --- Solve an eigenvalue problem with a Robin condition at the left end of the interval: ```wl In[1]:= sol = DSolve[{y''[x] + λ y[x] == 0, y[0] + y'[0] == 0, y[1] == 0}, y[x], x] Out[1]= {{y[x] -> Piecewise[{{C[1]*((-Sqrt[λ])*Cos[x*Sqrt[λ]] + Sin[x*Sqrt[λ]]), Sqrt[λ]*Cos[Sqrt[λ]] == Sin[Sqrt[λ]]}}, 0]}} ``` Locate the roots of the transcendental eigenvalue equation in the range ``10 < λ < 80`` : ```wl In[2]:= roots = λ /. Solve[Sqrt[λ] Cos[Sqrt[λ]] == Sin[Sqrt[λ]] && 10 < λ < 80, λ];data = Table[{λ, 0}, {λ, roots}]; In[3]:= Plot[Sqrt[λ] Cos[Sqrt[λ]] - Sin[Sqrt[λ]], {λ, 10, 80}, Epilog -> {PointSize[Large], Red, Point[data]}] Out[3]= [image] ``` Obtain the eigenfunctions within this range by using ``Assumptions`` : ```wl In[4]:= DSolve[{y''[x] + λ y[x] == 0, y[0] + y'[0] == 0, y[1] == 0}, y[x], x, Assumptions -> 1 < λ < 80] Out[4]= {{y[x] -> Piecewise[{{C[1]*((-Sqrt[λ])*Cos[x*Sqrt[λ]] + Sin[x*Sqrt[λ]]), λ == Root[{-Sin[Sqrt[#1]] + Cos[Sqrt[#1]]*Sqrt[#1] & , 20.19072855642662997452564796127538748749`20.298070440614282}] || λ == Root[{-Sin[Sqrt[#1]] + Cos[Sqrt[#1]]*Sqrt[#1] & , 59.6795159441094188805499668598608664638`20.30138232495084}]}}, 0]}} ``` --- Solve an eigenvalue problem with Robin boundary conditions at both ends: ```wl In[1]:= DSolve[{y''[x] + 2λ y'[x] + λ ^ 2 y[x] == 0, y[1] + y'[1] == 0, 3y[2] + 2y'[2] == 0}, y[x], x] Out[1]= {{y[x] -> Piecewise[{{((2 + x*(-1 + λ) - λ)*C[1])/(E^(x*λ)*(-1 + λ)), λ == 1/2 || λ == 2}}, 0]}} ``` --- Solve an eigenvalue problem for the Airy operator: ```wl In[1]:= eqns = {-y''[x] + x y[x] == λ y[x], y[0] == 0, y[1] == 0}; In[2]:= sol = y[x] /. DSolve[eqns, y[x], x, Assumptions -> 1 < λ < 200][[1]]; ``` Eigenvalues for the problem: ```wl In[3]:= eigvals = {ToRules[sol[[1, 1, 2]]]} Out[3]= {{λ -> Root[{AiryAi[-#1]*AiryBi[1 - #1] - AiryAi[1 - #1]*AiryBi[-#1] & , 10.3685071617373903698`7.5}]}, {λ -> Root[{AiryAi[-#1]*AiryBi[1 - #1] - AiryAi[1 - #1]*AiryBi[-#1] & , 39.9787447898614269817`7.5}]}, {λ -> Root[{AiryAi[-#1]*AiryBi[1 - #1] - AiryAi[1 - #1]*AiryBi[-#1] & , 89.3266345424780715158`9.5}]}, {λ -> Root[{AiryAi[-#1]*AiryBi[1 - #1] - AiryAi[1 - #1]*AiryBi[-#1] & , 158.41378981430977769698117916625548925432`12.}]}} ``` Eigenfunctions for this problem are given by the following: ```wl In[4]:= eigfuns = (sol[[1, 1, 1]] /. {C[1] -> 1}) Out[4]= AiryBi[x - λ] - (AiryAi[x - λ] AiryBi[-λ]/AiryAi[-λ]) ``` Plot the eigenfunctions for the range ``1 < λ < 200`` : ```wl In[5]:= Plot[Evaluate[eigfuns /. eigvals], {x, 0, 1}] Out[5]= [image] ``` #### Integral Equations (6) Solve a Volterra integral equation: ```wl In[1]:= eqn = y[x] == x^3 + λSubsuperscript[∫, 0, x](t - x)y[t]\[DifferentialD]t; In[2]:= sol = DSolve[eqn, y[x], x] Out[2]= {{y[x] -> (6 x/λ) - (6 Sin[x Sqrt[λ]]/λ^3 / 2)}} ``` Plot the solution: ```wl In[3]:= Plot[Table[y[x] /. sol[[1]], {λ, 1, 3, 0.5}]//Evaluate, {x, 0, 20}] Out[3]= [image] ``` --- Solve a Fredholm integral equation: ```wl In[1]:= eqn = y[x] == Sin[7 x] - (x/a) + Subsuperscript[∫, 0, (π/2)]x y[t]\[DifferentialD]t; In[2]:= sol = DSolve[eqn, y[x], x] Out[2]= {{y[x] -> (56 x - 8 a x - 56 a Sin[7 x] + 7 a π^2 Sin[7 x]/7 a (-8 + π^2))}} ``` Plot the solution: ```wl In[3]:= Plot[Table[y[x] /. sol[[1]], {a, -1, 4, 0.7}]//Evaluate, {x, 0, Pi / 2}] Out[3]= [image] ``` --- Solve an integro-differential equation: ```wl In[1]:= eqn = Derivative[1][y][x] == 1 + Sin[a x] + Subsuperscript[∫, 0, x]y[t]\[DifferentialD]t; ``` Obtain the general solution: ```wl In[2]:= sol1 = DSolve[eqn, y[x], x] Out[2]= {{y[x] -> E^-x (-1 + (1 + E^2 x) C[1]) - (a Cos[a x]/1 + a^2)}} ``` Specify an initial condition to obtain a particular solution: ```wl In[3]:= init = y[0] == -1; In[4]:= sol2 = DSolve[{eqn, init}, y[x], x] Out[4]= {{y[x] -> (E^-x (-2 + a - 2 a^2 + a E^2 x - 2 a E^x Cos[a x])/2 (1 + a^2))}} ``` Plot the solution: ```wl In[5]:= Plot[Table[y[x] /. sol2[[1]], {a, -1, 4, 0.7}]//Evaluate, {x, 0, 3}] Out[5]= [image] ``` --- Solve a singular Abel integral equation: ```wl In[1]:= eqn = Sin[a x] == Subsuperscript[∫, 0, x](y[t]/Sqrt[x - t])\[DifferentialD]t; In[2]:= sol = DSolve[eqn, y[x], x] Out[2]= {{y[x] -> Sqrt[a] Sqrt[(2/π)] (Cos[a x] FresnelC[Sqrt[(2/π)] Sqrt[a x]] + FresnelS[Sqrt[(2/π)] Sqrt[a x]] Sin[a x])}} ``` Plot the solution: ```wl In[3]:= Plot[Table[y[x] /. sol[[1]], {a, 1, 3, 0.7}]//Evaluate, {x, 0, 9}] Out[3]= [image] ``` --- Solve a weakly singular Volterra integral equation: ```wl In[1]:= eqn = y[x] == x^a - Subsuperscript[∫, 0, x](y[t]/Sqrt[x - t])\[DifferentialD]t; In[2]:= sol = DSolve[eqn, y[x], x] Out[2]= {{y[x] -> x^a + (E^π x π^-a Gamma[1 + a] Gamma[(1/2) + a, π x]/Gamma[(1/2) + a]) - E^π x π^-a Gamma[1 + a, π x]}} ``` Plot the solution: ```wl In[3]:= Plot[Table[y[x] /. sol[[1]], {a, 1, 4, 0.7}]//Evaluate, {x, 0, 2}] Out[3]= [image] ``` --- Solve a homogeneous Fredholm equation of the second kind: ```wl In[1]:= eqn = y[x] == λ Subsuperscript[∫, 0, π]Cos[x + t]y[t]\[DifferentialD]t; In[2]:= sol = DSolve[eqn, y[x], x] Out[2]= {{y[x] -> Piecewise[{{(2*C[2]*Sin[x])/Pi, λ == -(2/Pi)}, {(2*C[1]*Cos[x])/Pi, λ == 2/Pi}}, 0]}} ``` Plot the solution: ```wl In[3]:= Plot[Table[y[x] /. sol[[1]] /. { C[1] -> 1, C[2] -> 1}, {λ, {-2 / Pi, 2 / Pi}}]//Evaluate, {x, 0, 3π}] Out[3]= [image] ``` #### First-Order Partial Differential Equations (7) General solution for a linear first-order partial differential equation: ```wl In[1]:= DSolve[3D[u[x, y], x] + 5D[u[x, y], y] == x , u, {x, y}] Out[1]= {{u -> Function[{x, y}, (1/6) (x^2 + 6 C[1][(1/3) (-5 x + 3 y)])]}} ``` The solution with a particular choice of the arbitrary function ``C[1]`` : ```wl In[2]:= (u[x, y] /. %[[1]]) /. {C[1][e_] :> Sin[e] + Cos[e]} Out[2]= (1/6) (x^2 + 6 (Cos[(1/3) (-5 x + 3 y)] + Sin[(1/3) (-5 x + 3 y)])) In[3]:= Plot3D[%, {x, -5, 5}, {y, -5, 5}] Out[3]= [image] ``` --- Initial value problem for a linear first-order partial differential equation: ```wl In[1]:= DSolve[{x D[u[x, y], y] + y D[u[x, y], x] == -4x y u[x, y], u[x, 0] == E ^ (-x ^ 2)}, u, {x, y}] Out[1]= {{u -> Function[{x, y}, E^-x^2 - y^2]}} In[2]:= Plot3D[u[x, y] /. %[[1]], {x, -2, 2}, {y, -2, 2}, ColorFunction -> "TemperatureMap"] Out[2]= [image] ``` --- Initial-boundary value problem for a linear first-order partial differential equation: ```wl In[1]:= DSolve[{D[u[t, x], t] + D[u[t, x], x] == 0, u[t, 0] == 0, u[0, x] == Sin[x]}, u, {t, x}] Out[1]= {{u -> Function[{t, x}, -HeavisideTheta[-t + x] Sin[t - x]]}} In[2]:= Plot3D[u[t, x] /. %[[1]]//Evaluate, {t, 0, 3}, {x, 0, 3}, Exclusions -> None] Out[2]= [image] ``` --- General solution for the transport equation: ```wl In[1]:= DSolve[D[u[t, x], t] + c D[u[t, x], x] == 0, u, {t, x}] Out[1]= {{u -> Function[{t, x}, C[1][-c t + x]]}} ``` Initial value problem: ```wl In[2]:= DSolve[{D[u[t, x], t] + c D[u[t, x], x] == 0, u[0, x] == E ^ (-x ^ 2)}, u, {t, x}] Out[2]= {{u -> Function[{t, x}, E^-(-c t + x)^2]}} ``` Plot the traveling wave with speed ``c=1`` : ```wl In[3]:= Plot[u[t, x] /. %[[1]] /. {c -> 1, t -> 3}//Evaluate, {x, 0, 10}, PlotRange -> All, Ticks -> False, Filling -> Axis, Epilog -> Arrow[{{4, 4 / 5}, {7, 4 / 5}}]] Out[3]= [image] ``` --- General solution for a quasilinear first-order partial differential equation: ```wl In[1]:= DSolve[2D[u[x, y], x] + 5D[u[x, y], y] == u[x, y] ^ 2 + 1, u, {x, y}]//Quiet Out[1]= {{u -> Function[{x, y}, Tan[(1/2) (x + 2 C[1][(1/2) (-5 x + 2 y)])]]}} ``` --- Initial value problem for a scalar conservation law: ```wl In[1]:= DSolve[{D[u[x, y], x] + u[x, y] D[u[x, y], y] == 0, u[x, 0] == 1 / (x + 1)}, u, {x, y}] Out[1]= {{u -> Function[{x, y}, (1 + y/1 + x)]}} In[2]:= Plot3D[u[x, y] /. %[[1]], {x, 0, 3}, {y, 0, 2}] Out[2]= [image] ``` --- Complete integral for a nonlinear first-order Clairaut equation: ```wl In[1]:= DSolve[u[x, y] == x D[u[x, y], x] + y D[u[x, y], y] + Sin[D[u[x, y], x] + D[u[x, y], y]] , u, {x, y}]//Quiet Out[1]= {{u -> Function[{x, y}, x C[1] + y C[2] + Sin[C[1] + C[2]]]}} ``` #### Hyperbolic Partial Differential Equations (11) Initial value problem for the wave equation: ```wl In[1]:= weqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}]; In[2]:= ic = {u[x, 0] == E ^ (-x ^ 2), Derivative[0, 1][u][x, 0] == 1}; In[3]:= DSolve[{weqn, ic}, u, {x, t}] Out[3]= {{u -> Function[{x, t}, Piecewise[{{(1/2)*(E^(-(-t + x)^2) + E^(-(t + x)^2)) + t, t >= 0}}, Indeterminate]]}} ``` The wave propagates along a pair of characteristic directions: ```wl In[4]:= Plot3D[Evaluate[u[x, t] /. %[[1]]], {x, -7, 7}, {t, 0, 4}, PlotRange -> All] Out[4]= [image] ``` --- Initial value problem for the wave equation with piecewise initial data: ```wl In[1]:= weqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}]; In[2]:= ic = {u[x, 0] == UnitBox[x] + UnitTriangle[x / 3], Derivative[0, 1][u][x, 0] == 0}; In[3]:= DSolve[ {weqn, ic}, u[x, t], {x, t}]//Simplify//TraditionalForm Out[3]//TraditionalForm= $$\left\{\left\{u(x,t)\to \begin{array}{ll} \{ & \begin{array}{ll} \frac{1}{2} \left(\Pi (t-x)+\Pi (t+x)+\Lambda \left(\frac{x-t}{3}\right)+\Lambda \left(\frac{t+x}{3}\right)\right) & t\geq 0 \\ \text{Indeterminate} & \text{True} \\ \end{array} \\ \end{array} \right\}\right\}$$ ``` Discontinuities in the initial data are propagated along the characteristic directions: ```wl In[4]:= Plot3D[Evaluate[u[x, t] /. %[[1]]], {x, -7, 7}, {t, 0, 4}, PlotRange -> All, PlotPoints -> 250, Exclusions -> None] Out[4]= [image] ``` --- Initial value problem with a pair of decaying exponential functions as initial data: ```wl In[1]:= weqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}]; In[2]:= ic = {u[x, 0] == E ^ (-(x - 6) ^ 2) + E ^ (-(x + 6) ^ 2), Derivative[0, 1][u][x, 0] == 1 / 2}; In[3]:= sol = DSolve[ {weqn, ic}, u[x, t], {x, t}] Out[3]= {{u[x, t] -> Piecewise[{{(1/2)*(E^(-(-6 - t + x)^2) + E^(-(6 - t + x)^2) + E^(-(-6 + t + x)^2) + E^(-(6 + t + x)^2)) + t/2, t >= 0}}, Indeterminate]}} ``` Visualize the solution: ```wl In[4]:= Plot3D[Evaluate[u[x, t] /. sol[[1]]], {x, -30, 30}, {t, 0, 20}, PlotRange -> All, Exclusions -> All, PlotPoints -> 30, WorkingPrecision -> 20]//Quiet Out[4]= [image] ``` --- Initial value problem for an inhomogeneous wave equation: ```wl In[1]:= weqn = {D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}] + m}; In[2]:= ic = {u[x, 0] == Sin[x] - Cos[3 * x] / E ^ (Abs[x] / 6), Derivative[0, 1][u][x, 0] == 0}; In[3]:= sol = DSolve[{weqn, ic}, u[x, t], {x, t}] Out[3]= {{u[x, t] -> Piecewise[{{(m*t^2)/2 + (1/2)*((-E^((-(1/6))*Abs[-t + x]))*Cos[3*(-t + x)] - Cos[3*(t + x)]/E^((1/6)*Abs[t + x]) - Sin[t - x] + Sin[t + x]), t >= 0}}, Indeterminate]}} ``` Visualize the solution for different values of ``m`` : ```wl In[4]:= Table[Plot3D[Evaluate[u[x, t] /. sol[[1]]], {x, -7, 7}, {t, 0, 4}, PlotRange -> All, Ticks -> False], {m, {-1, 0, 1}}] Out[4]= [image] ``` --- Initial value problem for the wave equation with a Dirichlet condition on the half-line: ```wl In[1]:= weqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}]; In[2]:= ic = {u[x, 0] == Piecewise[{{Sin[x] ^ 2, Pi < x < 2Pi}}], Derivative[0, 1][u][x, 0] == 0}; In[3]:= bc = u[0, t] == 0; In[4]:= sol = DSolve[{weqn, ic, bc}, u, {x, t}]; ``` Visualize the solution: ```wl In[5]:= Plot3D[Evaluate[u[x, t] /. sol[[1]]], {x, 0, 12}, {t, 0, 10}, Exclusions -> None, PlotRange -> All, PlotPoints -> 120] Out[5]= [image] ``` The wave bifurcates starting from the initial data: ```wl In[6]:= Grid[Partition[Table[Plot[Evaluate[u[x, t] /. sol[[1]]], {x, 0, 20}, Exclusions -> None, PlotRange -> All], {t, 0, 12}], 3], ItemSize -> 10] Out[6]= | | | | | ------- | ------- | ------- | | [image] | [image] | [image] | | [image] | [image] | [image] | | [image] | [image] | [image] | | [image] | [image] | [image] | ``` --- Initial value problem for the wave equation with a Neumann condition on the half-line: ```wl In[1]:= weqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}]; In[2]:= ic = {u[x, 0] == Sin[x] ^ 3, Derivative[0, 1][u][x, 0] == 1 - E ^ (-x / 10)}; In[3]:= bc = Derivative[1, 0][u][0, t] == 1; In[4]:= sol = DSolve[{weqn, ic , bc }, u, {x, t}]; ``` Visualize the solution: ```wl In[5]:= Plot3D[Evaluate[u[x, t] /. sol[[1]]], {x, 0, 20}, {t, 0, 2}, Exclusions -> None, PlotRange -> All] Out[5]= [image] ``` In this example, the wave evolves to a non-oscillating function: ```wl In[6]:= Grid[Partition[Table[Plot[Evaluate[u[x, t] /. sol[[1]]], {x, 0, 20}, Exclusions -> None], {t, 0, 12}], 3], ItemSize -> 8] Out[6]= [image] ``` --- Dirichlet problem for the wave equation on a finite interval: ```wl In[1]:= weqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}]; In[2]:= ic = {u[x, 0] == x ^ 2(Pi - x), Derivative[0, 1][u][x, 0] == 0}; In[3]:= bc = {u[0, t] == 0, u[Pi, t] == 0}; ``` The solution is an infinite trigonometric series: ```wl In[4]:= sol = DSolve[{weqn, ic, bc}, u, {x, t}] Out[4]= {{u -> Function[{x, t}, Inactive[Sum][-((4*(1 + 2*(-1)^K[1])*Cos[t*K[1]]*Sin[x*K[1]])/K[1]^3), {K[1], 1, Infinity}]]}} ``` Extract the first three terms from the ``Inactive`` sum: ```wl In[5]:= asol = TruncateSum[u[x, t] /. sol[[1]], 3] Out[5]= 4 Cos[t] Sin[x] - (3/2) Cos[2 t] Sin[2 x] + (4/27) Cos[3 t] Sin[3 x] In[6]:= Plot3D[asol, {x, 0, Pi}, {t, 0, 4Pi}, ColorFunction -> Hue] Out[6]= [image] ``` The wave executes periodic motion in the vertical direction: ```wl In[7]:= Plot[Evaluate[Table[asol, {t, 0, 7}]], {x, 0, Pi}] Out[7]= [image] ``` --- Dirichlet problem for the wave equation in a rectangle: ```wl In[1]:= weqn = D[u[x, y, t], {t, 2}] == Laplacian[u[x, y, t], {x, y}]; In[2]:= ic = {u[x, y, 0] == (1 / 10)(x - x ^ 2)(2y - y ^ 2), Derivative[0, 0, 1][u][x, y, 0] == 0}; In[3]:= bc = {u[x, 0, t] == 0, u[0, y, t] == 0, u[1, y, t] == 0, u[x, 2, t] == 0}; ``` The solution is a doubly infinite trigonometric series: ```wl In[4]:= (sol = FullSimplify[u[x, y, t] /. DSolve[{weqn, ic, bc}, u, {x, y, t}][[1]], K[1]∈Integers && K[1] ≥ 1 && K[3]∈Integers && K[3] ≥ 1])//TraditionalForm Out[4]//TraditionalForm= $$\underset{K[1]=1}{\overset{\infty }{\sum }}\underset{K[3]=1}{\overset{\infty }{\sum }}\frac{32 \left((-1)^{K[1]}-1\right) \left((-1)^{K[3]}-1\right) \cos \left(\frac{1}{2} \pi t \sqrt{4 K[1]^2+K[3]^2}\right) \sin (\pi x K[1]) \sin \left(\frac{1}{2} \pi y K[3]\right)}{5 \pi ^6 K[1]^3 K[3]^3}$$ ``` Extract a few terms from the ``Inactive`` sum: ```wl In[5]:= h[x_, y_, t_] = TruncateSum[sol, 3] Out[5]= (128 Cos[(1/2) Sqrt[5] π t] Sin[π x] Sin[(π y/2)]/5 π^6) + (128 Cos[(1/2) Sqrt[37] π t] Sin[3 π x] Sin[(π y/2)]/135 π^6) + (128 Cos[(1/2) Sqrt[13] π t] Sin[π x] Sin[(3 π y/2)]/135 π^6) + (128 Cos[(3/2) Sqrt[5] π t] Sin[3 π x] Sin[(3 π y/2)]/3645 π^6) ``` The two-dimensional wave executes periodic motion in the vertical direction: ```wl In[6]:= Animate[Plot3D[h[x, y, t], {x, 0, 1}, {y, 0, 2}, Ticks -> False, PlotRange -> {-1 / 36, 1 / 36}, MeshStyle -> Red, PerformanceGoal -> "Quality"], {t, 0, 8}, SaveDefinitions -> True, DefaultDuration -> 12] Out[6]= DynamicModule[«9»] ``` --- Dirichlet problem for the wave equation in a circular disk: ```wl In[1]:= reqn = r D[u[r, t], {t, 2}] == D[r D[u[r, t], r], r]; In[2]:= ic = {u[r, 0] == 1, Derivative[0, 1][u][r, 0] == r / 3}; In[3]:= bc = u[1, t] == 0; ``` The solution is an infinite Bessel series: ```wl In[4]:= (sol = u[r, t] /. DSolve[{reqn, ic, bc}, u[r, t], {r, t}][[1]]//FullSimplify)//TraditionalForm Out[4]//TraditionalForm= $$\underset{K[1]=1}{\overset{\infty }{\sum }}\frac{2 J_0\left(r j_{0,K[1]}\right) \left(\sin \left(t j_{0,K[1]}\right) _1F_2\left(\frac{3}{2};1,\frac{5}{2};-\frac{1}{4} \left(j_{0,K[1]}\right){}^2\right)+9 J_1\left(j_{0,K[1]}\right) \cos \left(t j_{0,K[1]}\right)\right)}{9 j_{0,K[1]} \left(J_0\left(j_{0,K[1]}\right){}^2+J_1\left(j_{0,K[1]}\right){}^2\right)}$$ ``` Extract the first three terms from the ``Inactive`` sum: ```wl In[5]:= h[r_, t_] = TruncateSum[sol, 3]//N; In[6]:= Plot3D[h[r, 0.7] /. {r -> Sqrt[x ^ 2 + y ^ 2]}, {x, y}∈Disk[], PlotRange -> All, Ticks -> None] Out[6]= [image] ``` --- General solution for a second-order hyperbolic partial differential equation: ```wl In[1]:= DSolve[12 D[u[x, t], {x, 2}] == D[u[x, t], {t, 2}] + D[u[x, t], x, t], u, {x, t}] Out[1]= {{u -> Function[{x, t}, C[1][t - (x/4)] + C[2][t + (x/3)]]}} ``` Hyperbolic partial differential equation with non-rational coefficients: ```wl In[2]:= DSolve[D[u[x, y], {x, 2}] - 2Sin[x] D[u[x, y], x, y] - Cos[x] ^ 2 D[u[x, y], {y, 2}] - Cos[x]D[u[x, y], y] == 0, u, {x, y} ] Out[2]= {{u -> Function[{x, y}, C[1][x + y - Cos[x]] + C[2][-x + y - Cos[x]]]}} ``` Inhomogeneous hyperbolic partial differential equation with constant coefficients: ```wl In[3]:= DSolve[3D[u[x, t], {x, 2}] - D[u[x, t], {t, 2}] + D[u[x, t], x, t] == 1, u, {x, t}] Out[3]= {{u -> Function[{x, t}, (x^2/6) + C[1][t - (1/6) (1 + Sqrt[13]) x] + C[2][t - (1/6) (1 - Sqrt[13]) x]]}} ``` --- Initial value problem for an inhomogeneous linear hyperbolic system with constant coefficients: ```wl In[1]:= eqns = {D[u[x, t], t] == D[v[x, t], x] + 1, D[v[x, t], t] == -D[u[x, t], x] - 1}; In[2]:= ic = {u[x, 0] == Cos[x] ^ 2, v[x, 0] == Sin[x]}; In[3]:= sol = DSolve[{eqns, ic}, {u[x, t], v[x, t]}, {x, t}]//FullSimplify Out[3]= {{u[x, t] -> (1/2) + t + (1/2) Cos[2 x] Cosh[2 t] + Cos[x] Sinh[t], v[x, t] -> -t + Cosh[t] Sin[x] (1 + 2 Cos[x] Sinh[t])}} In[4]:= Plot3D[({u[x, t], v[x, t]} /. sol[[1]])//Evaluate, {x, 0, 4}, {t, 0, 3}, PlotRange -> {-70, 120}] Out[4]= [image] ``` #### Parabolic Partial Differential Equations (7) Initial value problem for the heat equation: ```wl In[1]:= heqn = D[u[x, t], t] == D[u[x, t], {x, 2}]; In[2]:= ic = u[x, 0] == E ^ (-x ^ 2); In[3]:= sol = DSolve[{heqn, ic }, u[x, t], {x, t}] Out[3]= {{u[x, t] -> ConditionalExpression[1/(E^(x^2/(1 + 4*t))*(Sqrt[4 + 1/t]*Sqrt[t])), Re[1/t] >= -4]}} In[4]:= Plot3D[Evaluate[u[x, t] /. sol[[1]]], {x, -5, 5}, {t, 0, 4}, PlotRange -> All] Out[4]= [image] ``` Visualize the diffusion of heat with the passage of time: ```wl In[5]:= Plot[Evaluate[Table[u[x, t] /. sol[[1]], {t, 0, 4}]], {x, -5, 5}, PlotRange -> All, Filling -> Axis]//Quiet Out[5]= [image] ``` --- Initial value problem for the heat equation with piecewise initial data: ```wl In[1]:= heqn = D[u[x, t], t] == D[u[x, t], {x, 2}]; In[2]:= ic = u[x, 0] == UnitBox[x]; ``` The solution is given in terms of the error function ``Erf`` : ```wl In[3]:= sol = DSolve[{heqn, ic }, u[x, t], {x, t}] Out[3]= {{u[x, t] -> (1/2) (Erf[(1 - 2 x/4 Sqrt[t])] + Erf[(1 + 2 x/4 Sqrt[t])])}} ``` Discontinuities in the initial data are smoothed instantly: ```wl In[4]:= Plot3D[Evaluate[u[x, t] /. sol[[1]]], {x, -2, 2}, {t, 0, 1}, PlotRange -> All, PlotPoints -> 250] Out[4]= [image] ``` --- Initial value problem for an inhomogeneous heat equation: ```wl In[1]:= heqn = D[u[x, t], t] == D[u[x, t], {x, 2}] + m; In[2]:= ic = u[x, 0] == Sin[x]; In[3]:= sol = DSolve[{heqn, ic }, u[x, t], {x, t}] Out[3]= {{u[x, t] -> m t + E^-t Sin[x]}} ``` Visualize the growth of the solution for different values of the parameter ``m`` : ```wl In[4]:= Table[Plot3D[Evaluate[u[x, t] /. sol[[1]]] /. {m -> i}, {x, -7, 7}, {t, 0, 4}, PlotRange -> All, Axes -> False], {i, {2, 0, -2}}] Out[4]= [image] ``` --- Dirichlet problem for the heat equation on a finite interval: ```wl In[1]:= heqn = D[u[x, t], t] == D[u[x, t], {x, 2}]; In[2]:= ic = u[x, 0] == x(3 - x) ^ 2; In[3]:= bc = {u[0, t] == 0, u[3, t] == 0}; ``` The solution is a Fourier sine series: ```wl In[4]:= sol = DSolve[{heqn, ic, bc }, u[x, t], {x, t}] Out[4]= {{u[x, t] -> Inactive[Sum][(108*(2 + (-1)^K[1])*Sin[(1/3)*Pi*x*K[1]])/(E^((1/9)*Pi^2*t*K[1]^2)*(Pi^3*K[1]^3)), {K[1], 1, Infinity}]}} ``` Extract three terms from the ``Inactive`` sum: ```wl In[5]:= asol = TruncateSum[u[x, t] /. sol[[1]], 3] Out[5]= (108 E^-(π^2 t/9) Sin[(π x/3)]/π^3) + (81 E^-(4 π^2 t/9) Sin[(2 π x/3)]/2 π^3) + (4 E^-π^2 t Sin[π x]/π^3) In[6]:= Plot3D[asol//Evaluate, {x, 0, 3}, {t, 0, 3 / 2}, Exclusions -> None, PlotRange -> All] Out[6]= [image] ``` --- Neumann problem for the heat equation on a finite interval: ```wl In[1]:= heqn = D[u[x, t], t] == D[u[x, t], {x, 2}]; In[2]:= ic = u[x, 0] == x(3 - x); In[3]:= bc = {Derivative[1, 0][u][0, t] == 0, Derivative[1, 0][u][3, t] == 0}; ``` The solution is a Fourier cosine series: ```wl In[4]:= sol = u[x, t] /. DSolve[{heqn, ic, bc }, u[x, t], {x, t}][[1]] Out[4]= (3/2) + (2/3) Inactive[Sum][-((27*(1 + (-1)^K[1])*Cos[(1/3)*Pi*x*K[1]])/(E^((1/9)*Pi^2*t*K[1]^2)*(Pi^2*K[1]^2))), {K[1], 1, Infinity}] ``` Extract a few terms from the ``Inactive`` sum: ```wl In[5]:= asol = TruncateSum[sol, 4] Out[5]= (3/2) + (2/3) (-(27 E^-(4 π^2 t/9) Cos[(2 π x/3)]/2 π^2) - (27 E^-(16 π^2 t/9) Cos[(4 π x/3)]/8 π^2)) In[6]:= Plot3D[asol//Evaluate, {x, 0, 3}, {t, 0, 1}, Exclusions -> None, PlotRange -> All] Out[6]= [image] ``` The solution approaches $\frac{3}{2}$ as $t\to \infty$ : ```wl In[7]:= asol /. t -> ∞ Out[7]= (3/2) ``` Visualize the evolution of the solution to its steady state: ```wl In[8]:= Plot[Evaluate[Table[asol, {t, 0, 1, 0.2}]], {x, 0, 3}, PlotRange -> All, Filling -> Axis] Out[8]= [image] ``` --- Dirichlet problem for the heat equation in a disk: ```wl In[1]:= rheqn = r D[u[r, t], t] == D[r D[u[r, t], r], r]; In[2]:= ic = u[r, 0] == 1 - r; In[3]:= bc = u[1, t] == 0; ``` The solution is an infinite Bessel series: ```wl In[4]:= (sol = u[r, t] /. DSolve[{rheqn, ic, bc}, u[r, t], {r, t}][[1]])//TraditionalForm Out[4]//TraditionalForm= $$\underset{K[1]=1}{\overset{\infty }{\sum }}\frac{\pi e^{-t \left(j_{0,K[1]}\right){}^2} J_0\left(r j_{0,K[1]}\right) \left(\pmb{H}_0\left(j_{0,K[1]}\right) J_1\left(j_{0,K[1]}\right)-\pmb{H}_1\left(j_{0,K[1]}\right) J_0\left(j_{0,K[1]}\right)\right)}{\left(j_{0,K[1]}\right){}^2 J_1\left(j_{0,K[1]}\right){}^2}$$ ``` Extract a few terms from the ``Inactive`` sum: ```wl In[5]:= h[r_, t_] = TruncateSum[sol, 3]//N; ``` Visualize the individual terms of the solution at time ``t = 0.1`` : ```wl In[6]:= Table[Plot3D[h[r, 0.1][[i]] /. {r -> Sqrt[x ^ 2 + y ^ 2]}//Evaluate, {x, y}∈Disk[], Ticks -> None], {i, 3}] Out[6]= [image] ``` --- Boundary value problem for the Black–Scholes equation: ```wl In[1]:= BlackScholesPDE = D[v[t, s], t] + 1 / 2σ ^ 2 s ^ 2 D[v[t, s], {s, 2}] + (r - q) s D[v[t, s], s] - r v[t, s] == 0; In[2]:= bc = v[T, s] == ψ[s]; In[3]:= DSolve[{BlackScholesPDE, bc}, v[t, s], {t, s}] Out[3]= {{v[t, s] -> (E^r (t - T) Subsuperscript[∫, -∞, ∞]E^-(((-t + T) (-q + r - (σ^2/2)) - K[1] + Log[s])^2/2 (-t + T) σ^2) ψ[E^K[1]]\[DifferentialD]K[1]/Sqrt[2 π] Sqrt[(-t + T) σ^2])}} ``` #### Elliptic Partial Differential Equations (9) Dirichlet problem for the Laplace equation in the upper half-plane: ```wl In[1]:= leqn = Laplacian[u[x, y], {x, y}] == 0; In[2]:= bc = u[x, 0] == UnitBox[x]; In[3]:= sol = DSolve[{leqn, bc}, u[x, y], {x, y}] Out[3]= {{u[x, y] -> Piecewise[ {{Piecewise[{{ArcTan[(1/2 - x)/y] + ArcTan[(1/2 + x)/y], y > 0 || x > 1/2 || x < -(1/2)}}, 0]/Pi, y >= 0}}, Indeterminate]}} ``` Discontinuities in the boundary data are smoothed out: ```wl In[4]:= Plot3D[u[x, y] /. sol[[1]], {x, -3, 3}, {y, 0, 5}, PlotRange -> All] Out[4]= [image] ``` --- Dirichlet problem for the Laplace equation in the right half-plane: ```wl In[1]:= leqn = Laplacian[u[x, y], {x, y}] == 0; In[2]:= bc = u[0, y] == Sinc[y]; In[3]:= sol = DSolve[{leqn, bc}, u[x, y], {x, y}] Out[3]= {{u[x, y] -> Piecewise[{{(x + (x*Cos[y] - y*Sin[y])*(-Cosh[x] + Sinh[x]))/(x^2 + y^2), x >= 0}}, Indeterminate]}} ``` Visualize the solution: ```wl In[4]:= Plot3D[u[x, y] /. sol[[1]], {y, -10, 10}, {x, 0, 5}, PlotRange -> All] Out[4]= [image] ``` --- Dirichlet problem for the Laplace equation in the first quadrant: ```wl In[1]:= leqn = Laplacian[u[x, y], {x, y}] == 0; In[2]:= bc = {u[x, 0] == (-1 / ((x - 2) ^ 2 + 3)), u[0, y] == (1 / ((y - 3) ^ 2 + 1))}; In[3]:= sol = DSolve[{leqn, bc}, u[x, y], {x, y}]; ``` Visualize the solution: ```wl In[4]:= Plot3D[u[x, y] /. sol[[1]], {x, 0, 10}, {y, 0, 7}, PlotRange -> All] Out[4]= [image] ``` --- Neumann problem for the Laplace equation in the upper half-plane: ```wl In[1]:= leqn = Laplacian[u[x, y], {x, y}] == 0; In[2]:= bc = Derivative[0, 1][u][x, 0] == UnitBox[x]; In[3]:= sol = DSolve[{leqn, bc}, u[x, y], {x, y}] Out[3]= {{u[x, y] -> Piecewise[{{(1/(4*Pi))*(-4 + 4*y*ArcCot[(2*y)/(1 - 2*x)] + 4*y*ArcTan[(1/2 + x)/y] - 2*Log[4] + Log[(1 - 2*x)^2 + 4*y^2] - 2*x*Log[(1 - 2*x)^2 + 4*y^2] + Log[(1 + 2*x)^2 + 4*y^2] + 2*x*Log[(1 + 2*x)^2 + 4*y^2]), y >= 0}}, Indeterminate]}} ``` Visualize the solution: ```wl In[4]:= Plot3D[u[x, y] /. sol[[1]], {x, -3, 3}, {y, 0, 5}, PlotRange -> All] Out[4]= [image] ``` --- Dirichlet problem for the Laplace equation in a rectangle: ```wl In[1]:= leqn = Laplacian[u[x, y], {x, y}] == 0; In[2]:= bc = {u[x, 0] == x ^ 2(1 - x), u[x, 2] == 0, u[0, y] == 0, u[1, y] == 0}; ``` The solution is an infinite trigonometric series: ```wl In[3]:= sol = FullSimplify[u[x, y] /. DSolve[{leqn, bc}, u[x, y], {x, y}][[1]]] Out[3]= Inactive[Sum][-((4*(1 + 2*(-1)^K[1])*Csch[2*Pi*K[1]]*Sin[Pi*x*K[1]]*Sinh[Pi*(2 - y)*K[1]])/ (Pi^3*K[1]^3)), {K[1], 1, Infinity}] ``` Extract a few terms from the ``Inactive`` sum: ```wl In[4]:= asol = TruncateSum[sol, 4] Out[4]= (4 Csch[2 π] Sin[π x] Sinh[π (2 - y)]/π^3) - (3 Csch[4 π] Sin[2 π x] Sinh[2 π (2 - y)]/2 π^3) + (4 Csch[6 π] Sin[3 π x] Sinh[3 π (2 - y)]/27 π^3) - (3 Csch[8 π] Sin[4 π x] Sinh[4 π (2 - y)]/16 π^3) ``` Visualize the solution: ```wl In[5]:= Plot3D[asol, {x, 0, 1}, {y, 0, 2}, PlotRange -> All] Out[5]= [image] ``` --- Dirichlet problem for the Laplace equation in a disk: ```wl In[1]:= leqn = Laplacian[u[r, θ], {r, θ}, "Polar"] == 0; In[2]:= bc = u[3, θ] == Sin[6θ]; In[3]:= sol = DSolve[{leqn, bc}, u[r, θ], {r, θ}] Out[3]= {{u[r, θ] -> (1/729) r^6 Sin[6 θ]}} ``` Visualize the solution: ```wl In[4]:= Plot3D[Evaluate[u[r, θ] /. sol[[1]] /. {r -> Sqrt[x ^ 2 + y ^ 2], θ -> ArcTan[y / x]}], {x, y}∈Disk[{0, 0}, 3], PlotRange -> All, Exclusions -> None] Out[4]= [image] ``` --- Dirichlet problem for the Laplace equation in an annulus: ```wl In[1]:= leqn = D[u[r, θ], {r, 2}] + (1 / r) D[u[r, θ], r] + (1 / r ^ 2) D[u[r, θ], {θ, 2}] == 0; In[2]:= bc = {u[1, θ] == 0, u[2, θ] == 5}; In[3]:= sol = DSolve[{leqn, bc}, u[r, θ], {r, θ}] Out[3]= {{u[r, θ] -> Piecewise[{{(5*Log[r])/Log[2], 1 <= r <= 2}}, Indeterminate]}} ``` Visualize the solution: ```wl In[4]:= Plot3D[Evaluate[u[r, θ] /. sol[[1]] /. {r -> Sqrt[x ^ 2 + y ^ 2], θ -> ArcTan[y / x]}], {x, y}∈Annulus[{0, 0}, {1, 2}], PlotRange -> All, Exclusions -> None] Out[4]= [image] ``` --- Dirichlet problem for the Poisson equation in a rectangle: ```wl In[1]:= peqn = Laplacian[u[x, y], {x, y}] == 6x - 6y; In[2]:= bc = {u[x, 0] == 1 + 11x + x ^ 3, u[x, 2] == -7 + 11x + x ^ 3, u[0, y] == 1 - y ^ 3, u[4, y] == 109 - y ^ 3}; In[3]:= sol = DSolve[{peqn, bc}, u[x, y], {x, y}] Out[3]= {{u[x, y] -> 1 + 11 x + x^3 - y^3}} ``` --- Dirichlet problem for the Helmholtz equation in a rectangle: ```wl In[1]:= heqn = {Laplacian[u[x, y], {x, y}] + 5u[x, y] == 0}; In[2]:= bc = {u[x, 0] == UnitTriangle[x - 2], u[x, 2] == 0, u[0, y] == 0, u[4, y] == 0}; In[3]:= (sol = u[x, y] /. DSolve[{heqn, bc}, u[x, y], {x, y}][[1]])//TraditionalForm Out[3]//TraditionalForm= $$\underset{K[1]=1}{\overset{\infty }{\sum }}\frac{1}{\pi ^2 K[1]^2}64 \sin ^3\left(\frac{1}{8} \pi K[1]\right) \left(\cos \left(\frac{1}{8} \pi K[1]\right)+\cos \left(\frac{3}{8} \pi K[1]\right)\right) \text{csch}\left(\frac{1}{2} \sqrt{\pi ^2 K[1]^2-80}\right) \sin \left(\frac{1}{4} \pi x K[1]\right) \sinh \left(\frac{1}{4} (2-y) \sqrt{\pi ^2 K[1]^2-80}\right)$$ ``` Extract a finite number of terms from the ``Inactive`` sum: ```wl In[4]:= fsol = TruncateSum[sol, 30]; ``` Visualize the approximate solution: ```wl In[5]:= Plot3D[fsol//Evaluate, {x, 0, 4}, {y, 0, 2}, PlotRange -> All] Out[5]= [image] ``` #### General Partial Differential Equations (6) A potential-free Schrödinger equation with Dirichlet boundary conditions: ```wl In[1]:= eqn = I ℏ D[ψ[x, t], t] == (-ℏ^2/2m) D[ψ[x, t], {x, 2}]; DSolve[{eqn, ψ[a, t] == 0, ψ[b, t] == 0}, ψ[x, t], {x, t}] Out[1]= {{ψ[x, t] -> Inactive[Sum][(C[K[1]]*Sin[(Pi*(-a + x)*K[1])/(-a + b)])/ E^((I*Pi^2*t*ℏ*K[1]^2)/(2*(-a + b)^2*m)), {K[1], 1, Infinity}]}} ``` Extract the first four terms in the solution: ```wl In[2]:= u[x_, t_] = TruncateSum[ψ[x, t] /. First[%], 4] Out[2]= E^-(I π^2 t ℏ/2 (-a + b)^2 m) C[1] Sin[(π (-a + x)/-a + b)] + E^-(2 I π^2 t ℏ/(-a + b)^2 m) C[2] Sin[(2 π (-a + x)/-a + b)] + E^-(9 I π^2 t ℏ/2 (-a + b)^2 m) C[3] Sin[(3 π (-a + x)/-a + b)] + E^-(8 I π^2 t ℏ/(-a + b)^2 m) C[4] Sin[(4 π (-a + x)/-a + b)] ``` For any choice of the four constants ``C[k]``, ``ψ`` obeys the equation and boundary conditions: ```wl In[3]:= Simplify[{eqn, ψ[a, t] == 0, ψ[b, t] == 0} /. ψ -> u] Out[3]= {True, True, True} ``` --- Initial value problem for a Schrödinger equation with Dirichlet boundary conditions: ```wl In[1]:= eqn = I D[ψ[x, t], t] == -2 D[ψ[x, t], {x, 2}]; f[x_] := -350 + 155 x - 22 x^2 + x^3 sol = DSolve[{eqn, ψ[5, t] == 0, ψ[10, t] == 0, ψ[x, 2] == f[x]}, ψ, {x, t}] Out[1]= {{ψ -> Function[{x, t}, Inactive[Sum][(100*(7 + 8*(-1)^K[1])*Sin[(1/5)*Pi*(-5 + x)*K[1]])/ (E^((2/25)*I*Pi^2*(-2 + t)*K[1]^2)*(Pi^3*K[1]^3)), {K[1], 1, Infinity}]]}} ``` Define a family of partial sums of the solution: ```wl In[2]:= u[k_Integer] = ψ /. TruncateSum[First[sol], k] Out[2]= Function[{x, t}, Underoverscript[∑, K[1] = 1, k](100 (7 + 8 (-1)^K[1]) E^-(2/25) I π^2 (-2 + t) K[1]^2 Sin[(1/5) π (-5 + x) K[1]]/π^3 K[1]^3)] ``` For each ``k``, ``uk`` satisfies the differential equation: ```wl In[3]:= Table[eqn /. ψ -> u[k], {k, 4}]//Simplify Out[3]= {True, True, True, True} ``` The boundary conditions are also satisfied: ```wl In[4]:= Table[{u[k][5, t] == 0, u[k][10, t] == 0}, {k, 4}] Out[4]= {{True, True}, {True, True}, {True, True}, {True, True}} ``` The initial condition is only satisfied for ``u∞``, but there is rapid convergence at ``t == 2`` : ```wl In[5]:= Grid@Partition[Table[Plot[{u[k][x, 2], f[x]}, {x, 5, 10}, ImageSize -> 200], {k, 4}], 2] Out[5]= | | | | ------- | ------- | | [image] | [image] | | [image] | [image] | ``` --- Solve a Schrödinger equation with potential over the reals: ```wl In[1]:= eqn = I D[ψ[x, t], t] == -D[ψ[x, t], {x, 2}] + 2 x^2ψ[x, t]; In[2]:= sol = DSolve[{eqn, ψ[-∞, t] == 0, ψ[∞, t] == 0}, ψ[x, t], {x, t}] Out[2]= {{ψ[x, t] -> Inactive[Sum][E^(-(x^2/Sqrt[2]) - 2*I*Sqrt[2]*t*(1/2 + K[1]))*C[K[1]]*HermiteH[K[1], 2^(1/4)*x], {K[1], 0, Infinity}]}} ``` Extract the first two terms in the solution: ```wl In[3]:= u[x_, t_] = TruncateSum[ψ[x, t] /. First[sol], 2] Out[3]= E^-I Sqrt[2] t - (x^2/Sqrt[2]) C[0] + 2 2^1 / 4 E^-3 I Sqrt[2] t - (x^2/Sqrt[2]) x C[1] ``` For any values of the constants ``C[0]`` and ``C[1]``, the equation is satisfied: ```wl In[4]:= Simplify[eqn /. ψ -> u] Out[4]= True ``` The conditions at infinity are also satisfied: ```wl In[5]:= Limit[u[x, t], x -> ∞] Out[5]= 0 In[6]:= Limit[u[x, t], x -> -∞] Out[6]= 0 ``` Although the function is time dependent, its $L^2$-norm is constant: ```wl In[7]:= Integrate[Abs[u[x, t]]^2, {x, -∞, ∞}, Assumptions -> t∈Reals] Out[7]= (Sqrt[π] (Abs[C[0]]^2 + 2 Abs[C[1]]^2)/2^1 / 4) ``` --- Initial value problem for Burgers' equation with viscosity ``ν`` : ```wl In[1]:= BurgersEqn = D[u[x, t], t] + u[x, t]D[u[x, t], x] == ν D[u[x, t], {x, 2}]; In[2]:= sol = DSolve[{BurgersEqn, u[x, 0] == UnitBox[x]}, u[x, t], {x, t}];(dsol = u[x, t] /. sol[[1]]//FullSimplify)//TraditionalForm Out[2]//TraditionalForm= $$\frac{e^{\frac{t+1}{4 \nu }} \left(\text{erf}\left(\frac{2 t-2 x+1}{4 \sqrt{\nu t}}\right)-\text{erf}\left(\frac{2 t-2 x-1}{4 \sqrt{\nu t}}\right)\right)}{e^{\frac{t+1}{4 \nu }} \left(\text{erfc}\left(\frac{2 t-2 x-1}{4 \sqrt{\nu t}}\right)-\text{erfc}\left(\frac{2 t-2 x+1}{4 \sqrt{\nu t}}\right)\right)+e^{\frac{x}{2 \nu }} \left(\text{erfc}\left(\frac{1-2 x}{4 \sqrt{\nu t}}\right)+e^{\left.\frac{1}{2}\right/\nu } \text{erfc}\left(\frac{2 x+1}{4 \sqrt{\nu t}}\right)\right)}$$ ``` Plot the solution at different times for ``ν = 1 / 40`` : ```wl In[3]:= Plot[Table[dsol /. {ν -> 1 / 40}, {t, 1 / 100, 7}]//Evaluate, {x, -2, 3.3}, PlotRange -> All, Filling -> Axis, WorkingPrecision -> 20] Out[3]= [image] ``` Plot the solution for different values of ``ν`` : ```wl In[4]:= Table[Plot3D[dsol//Evaluate, {x, -2, 2}, {t, 0, 5}, PlotRange -> All, Ticks -> False], {ν, {2 / 3, 1 / 6, 1 / 14}}] Out[4]= [image] ``` --- Boundary value problem for the Tricomi equation: ```wl In[1]:= TricomiEqn = D[u[x, y], {x, 2}] + y D[u[x, y], {y, 2}] == 0; In[2]:= DSolve[{TricomiEqn, u[x, 0] == 0, Derivative[0, 1][u][x, 0] == x ^ 2}, u[x, y], {x, y}] Out[2]= {{u[x, y] -> -y (-x^2 + y)}} ``` Visualize the solution: ```wl In[3]:= Plot3D[u[x, y] /. %[[1]]//Evaluate, {x, -3, 3}, {y, 0, 5}] Out[3]= [image] ``` --- Traveling wave solution for the Korteweg–de Vries (KdV) equation: ```wl In[1]:= KdV = {D[u[x, t], t] + D[u[x, t], {x, 3}] + 6u[x, t]D[u[x, t], x] == 0}; In[2]:= sol = DSolve[KdV, u[x, t], {x, t}]//Quiet Out[2]= {{u[x, t] -> -(-8 C[1]^3 + C[2] + 12 C[1]^3 Tanh[x C[1] + t C[2] + C[3]]^2/6 C[1])}} ``` Obtain a particular solution for a fixed choice of arbitrary constants: ```wl In[3]:= psol = u[x, t] /. sol[[1]] /. {C[1] -> 1, C[2] -> -4, C[3] -> 1 / 4}//Simplify Out[3]= 2 - 2 Tanh[(1/4) - 4 t + x]^2 ``` The wave travels to the right without changing its shape: ```wl In[4]:= Plot[Evaluate[Table[psol, {t, 1, 12, 3}]], {x, -7, 35}, PlotRange -> All, Filling -> Axis, Epilog -> {Arrow[{{5, 1.6}, {9, 1.6}}], Arrow[{{17, 1.6}, {21, 1.6}}], Arrow[{{29, 1.6}, {33, 1.6}}]}] Out[4]= [image] ``` #### Partial Differential Equations over Regions (3) Dirichlet problem for the Laplace equation in a rectangle: ```wl In[1]:= leqn = Laplacian[u[x, y], {x, y}] == 0; In[2]:= dcond = DirichletCondition[u[x, y] == Piecewise[{{UnitTriangle[2 x - 1], y == 0 || y == 2}}, 0], True]; In[3]:= Ω = Rectangle[{0, 0}, {1, 2}]; In[4]:= (sol = FullSimplify[DSolve[{leqn, dcond}, u[x, y], {x, y}∈Ω]])//TraditionalForm Out[4]//TraditionalForm= $$\left\{\left\{u(x,y)\to \underset{K[1]=1}{\overset{\infty }{\sum }}\frac{8 \cosh (\pi (y-1) K[1]) \text{sech}(\pi K[1]) \sin \left(\frac{1}{2} \pi K[1]\right) \sin (\pi x K[1])}{\pi ^2 K[1]^2}\right\}\right\}$$ ``` Extract the first ``100`` terms from the ``Inactive`` sum: ```wl In[5]:= asol = TruncateSum[u[x, y] /. sol[[1]], 100]; ``` Visualize the solution on the rectangle: ```wl In[6]:= Plot3D[asol//Evaluate, {x, y}∈Ω, PlotRange -> All, PlotStyle -> Hue[0.2]] Out[6]= [image] ``` --- Dirichlet problem for the Laplace equation in a disk: ```wl In[1]:= leqn = Laplacian[u[x, y], {x, y}] == 0; In[2]:= dcond = DirichletCondition[u[x, y] == Sin[6ArcTan[y / x]], True]; In[3]:= Ω = Disk[{0, 0}, 3]; In[4]:= sol = DSolve[{leqn, dcond}, u[x, y], {x, y}∈Ω] Out[4]= {{u[x, y] -> (1/729) (x^2 + y^2)^3 Sin[6 ArcTan[(y/x)]]}} ``` Visualize the solution in the disk: ```wl In[5]:= Plot3D[u[x, y] /. sol[[1]]//Evaluate, {x, y}∈Ω, PlotRange -> All, PlotStyle -> Hue[0.5]] Out[5]= [image] ``` --- Dirichlet problem for the Laplace equation in a right half-plane: ```wl In[1]:= leqn = Laplacian[u[x, y], {x, y}] == 0; In[2]:= dcond = DirichletCondition[u[x, y] == 1 / (1 + x ^ 2), True]; In[3]:= Ω = HalfPlane[{{0, 0}, {1, 0}}, {0, 1}]; In[4]:= sol = DSolve[{leqn, dcond}, u[x, y], {x, y}∈Ω] Out[4]= {{u[x, y] -> (1 + y/x^2 + (1 + y)^2)}} ``` Visualize the solution in the half-plane: ```wl In[5]:= Plot3D[u[x, y] /. sol[[1]]//Evaluate, {x, y}∈Ω, PlotRange -> All, PlotStyle -> Green] Out[5]= [image] ``` #### Fractional Differential Equations (8) Solve a fractional differential equation of order 1/2: ```wl In[1]:= sol = DSolve[CaputoD[y[x], {x, 1 / 2}] - 3 y[x] == 0, y[x], x] Out[1]= {{y[x] -> E^9 x C[1] Erfc[-3 Sqrt[x]]}} ``` Verify the solution: ```wl In[2]:= CaputoD[y[x], {x, 1 / 2}] == 3y[x] /. sol//FullSimplify Out[2]= {True} ``` Add an initial condition: ```wl In[3]:= sol = DSolve[{CaputoD[y[x], {x, 1 / 2}] - 3 y[x] == 0, y[0] == 5}, y[x], x] Out[3]= {{y[x] -> 5 E^9 x Erfc[-3 Sqrt[x]]}} ``` Plot the solution: ```wl In[4]:= Plot[Evaluate[y[x] /. sol[[1]]], {x, 0, 1}] Out[4]= [image] ``` --- Solve a fractional differential equation containing ``CaputoD`` of order 0.7: ```wl In[1]:= sol = DSolve[CaputoD[y[x], {x, 0.7}] + 3 y[x] == 0, y[x], x] Out[1]= {{y[x] -> 1. C[1] MittagLefflerE[0.7, -3. x^0.7]}} ``` Verify the solution: ```wl In[2]:= CaputoD[y[x], {x, 0.7}] + 3y[x] == 0 /. sol//FullSimplify Out[2]= {True} ``` Add initial conditions and plot the solution: ```wl In[3]:= sol = DSolve[{CaputoD[y[x], {x, 0.7}] + 3 y[x] == 0, y[0] == 4}, y[x], x]; Plot[y[x] /. sol[[1]], {x, 0, 1}] Out[3]= [image] ``` --- Solve a fractional differential equation containing two Caputo derivatives of different orders: ```wl In[1]:= sol = DSolve[{CaputoD[y[x], {x, 3 / 4}] == -2CaputoD[y[x], {x, 1 / 2}] + y[x], y[0] == 1}, y[x], x] Out[1]= {{y[x] -> (1/x^3 / 4)2 (-(x^1 / 4/Sqrt[π]) - MittagLefflerE[(1/4), (1/4), -x^1 / 4] - (2 MittagLefflerE[(1/4), (1/4), (1/2) (-1 + Sqrt[5]) x^1 / 4]/-5 + Sqrt[5]) + (2 MittagLefflerE[(1/4), (1/4), -(1/2) (1 + Sqrt[5]) x^1 / 4]/5 + Sqrt[5])) + (1/5 x^3 / 4 Gamma[(1/4)])(-5 + 5 Gamma[(1/4)] MittagLefflerE[(1/4), (1/4), -x^1 / 4] + Sqrt[5] Gamma[(1/4)] MittagLefflerE[(1/4), (1/4), (1/2) (-1 + Sqrt[5]) x^1 / 4] - Sqrt[5] Gamma[(1/4)] MittagLefflerE[(1/4), (1/4), -(1/2) (1 + Sqrt[5]) x^1 / 4])}} ``` Plot the solution: ```wl In[2]:= Plot[y[x] /. sol[[1]], {x, 0, 1}] Out[2]= [image] ``` --- Solve a non-homogeneous fractional differential equation of order 1/7: ```wl In[1]:= sol = DSolve[{CaputoD[y[x], {x, 1 / 7}] == 8 (1 - y[x]), y[0] == 2}, y, x] Out[1]= {{y -> Function[{x}, 1 + MittagLefflerE[(1/7), -8 x^1 / 7]]}} ``` Verify the solution: ```wl In[2]:= {CaputoD[y[x], {x, 1 / 7}] == 8 (1 - y[x]), y[0] == 2} /. sol//Simplify Out[2]= {{True, True}} ``` --- Solve a system of two fractional differential equations: ```wl In[1]:= eqns = {CaputoD[x1[t], {t, 0.95}] == 2x1[t] - x2[t], CaputoD[x2[t], {t, 0.95}] == 4x1[t] - 3x2[t], x1[0] == 1.2, x2[0] == 4.2}; In[2]:= sol = DSolve[eqns, {x1, x2}, t] Out[2]= {{x1 -> Function[{t}, 1. MittagLefflerE[0.95, -2 t^0.95] + 0.2 MittagLefflerE[0.95, t^0.95]], x2 -> Function[{t}, 4. MittagLefflerE[0.95, -2 t^0.95] + 0.2 MittagLefflerE[0.95, t^0.95]]}} ``` Verify the solution: ```wl In[3]:= eqns /. sol//Simplify Out[3]= {{True, True, True, True}} ``` Parametrically plot the solution: ```wl In[4]:= ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol], {t, 0, 2}] Out[4]= [image] ``` --- Solve a system of two fractional ODEs in vector form: ```wl In[1]:= m = {{0, 1}, {-1, 0}}; a = 16 / 17; v = {-3, 5}; sol = DSolve[{CaputoD[x[t], {t, a}] == m. x[t], x[0] == v}, Element[x[t], Vectors[2]], t] Out[1]= {{x[t] -> {(-(3/2) + (5 I/2)) MittagLefflerE[(16/17), -I t^16 / 17] - ((3/2) + (5 I/2)) MittagLefflerE[(16/17), I t^16 / 17], ((5/2) + (3 I/2)) MittagLefflerE[(16/17), -I t^16 / 17] + ((5/2) - (3 I/2)) MittagLefflerE[(16/17), I t^16 / 17]}}} ``` Plot the solution: ```wl In[2]:= Plot[Evaluate[x[t] /. sol[[1]]], {t, 0, 10}] Out[2]= [image] ``` Parametrically plot the solution: ```wl In[3]:= ParametricPlot[Evaluate[x[t] /. sol[[1]]], {t, 0, 20}] Out[3]= [image] ``` --- Solve a system of three fractional differential equations in vector form: ```wl In[1]:= m = {{-1, 0, 0}, {2, 1, -9}, {3, 6, 1}}; a = 0.95; v = {-3, 5, 0}; In[2]:= sol = DSolve[{CaputoD[x[t], {t, a}] == m. x[t], x[0] == v}, Element[x[t], Vectors[3]], t] Out[2]= {{x[t] -> {(-3. + 0. I) MittagLefflerE[0.95, -t^0.95], (1.60345  + 0. I) MittagLefflerE[0.95, -t^0.95] + (1.69828  - 0.190047 I) MittagLefflerE[0.95, (1 - 3 I Sqrt[6]) t^0.95] + (1.69828  + 0.190047 I) MittagLefflerE[0.95, (1 + 3 I Sqrt[6]) t^0.95], (-0.310345 + 0. I) MittagLefflerE[0.95, -t^0.95] + (0.155172  + 1.38664 I) MittagLefflerE[0.95, (1 - 3 I Sqrt[6]) t^0.95] + (0.155172  - 1.38664 I) MittagLefflerE[0.95, (1 + 3 I Sqrt[6]) t^0.95]}}} ``` Plot the solution: ```wl In[3]:= Plot[Evaluate[x[t] /. sol[[1]]], {t, 0, 1}] Out[3]= [image] ``` --- Solve a fractional wave equation: ```wl In[1]:= a = 19 / 10; fweqn = CaputoD[y[t, x], {t, a}] == D[y[t, x], {x, 2}]; ic = {y[0, x] == Sin[Pi x], Derivative[1, 0][y][0, x] == 0, y[t, 0] == 0, y[t, 1] == 0}; In[2]:= sol = DSolve[{fweqn, ic}, y, {t, x}] Out[2]= {{y -> Function[{t, x}, MittagLefflerE[(19/10), -π^2 t^19 / 10] Sin[π x]]}} ``` Plot the solution: ```wl In[3]:= Plot3D[Evaluate[y[t, x] /. sol], {x, -5, 5}, {t, 0, 1}, PlotPoints -> 25] Out[3]= [image] ``` ### Generalizations & Extensions (1) No boundary condition; gives two generated parameters: ```wl In[1]:= DSolve[y''[x] + 4y[x] == 0, y, x] Out[1]= {{y -> Function[{x}, C[1] Cos[2 x] + C[2] Sin[2 x]]}} ``` One boundary condition: ```wl In[2]:= DSolve[{y''[x] + 4y[x] == 0, y[0] == 1}, y, x] Out[2]= {{y -> Function[{x}, Cos[2 x] + C[2] Sin[2 x]]}} ``` Two boundary conditions: ```wl In[3]:= DSolve[{y''[x] + 4y[x] == 0, y[0] == 1, y'[0] == 4}, y, x] Out[3]= {{y -> Function[{x}, Cos[2 x] + 2 Sin[2 x]]}} ``` ### Options (7) #### Assumptions (2) Solve an eigenvalue problem for a linear second-order differential equation: ```wl In[1]:= DSolve[{y''[x] + λy[x] == 0, y[0] == 0, y[π] == 0}, y[x], x] Out[1]= {{y[x] -> Piecewise[{{C[1]*Sin[x*Sqrt[λ]], Element[\[FormalN], Integers] && \[FormalN] >= 1 && λ == \[FormalN]^2}}, 0]}} ``` Use ``Assumptions`` to specify a range for the parameter ``λ`` : ```wl In[2]:= DSolve[{y''[x] + λy[x] == 0, y[0] == 0, y[π] == 0}, y[x], x, Assumptions -> 0 < λ < 20] Out[2]= {{y[x] -> Piecewise[{{C[1]*Sin[x*Sqrt[λ]], λ == 1 || λ == 4 || λ == 9 || λ == 16}}, 0]}} ``` --- Specify that the independent variable ``x`` is real: ```wl In[1]:= eqn = {y'[x] == Abs[x + Abs[x]] y[x], y[0] == 1}; In[2]:= sol = DSolve[eqn, y, x, Assumptions -> x∈Reals] Out[2]= {{y -> Function[{x}, Piecewise[{{1, x <= 0}}, E^x^2]]}} ``` Obtain the same result using an interval specification for ``x`` : ```wl In[3]:= DSolve[eqn, y, {x, -∞, ∞}] Out[3]= {{y -> Function[{x}, Piecewise[{{1, x <= 0}}, E^x^2]]}} ``` #### DiscreteVariables (1) Specify that a variable maintains its value between events: ```wl In[1]:= sol = DSolve[{x'[t] == a[t], x[0] == 0, a[0] == 2, WhenEvent[Mod[x[t], 1] == 0, a[t] -> -a[t]]}, {x, a}, {t, 0, 2}, DiscreteVariables -> a] Out[1]= {{x -> Function[{t}, Piecewise[{{2*t, 0 <= t <= 1/2}, {2 - 2*t, Inequality[1/2, Less, t, LessEqual, 1]}, {2*(-1 + t), Inequality[1, Less, t, LessEqual, 3/2]}, {4 - 2*t, Inequality[3/2, Less, t, LessEqual, 2]}}, Indeterminate]], a -> Function[{t}, Piecewise[{{2, 0 <= t <= 1/2}, {-2, Inequality[1/2, Less, t, LessEqual, 1]}, {2, Inequality[1, Less, t, LessEqual, 3/2]}, {-2, Inequality[3/2, Less, t, LessEqual, 2]}}, Indeterminate]]}} In[2]:= Plot[x[t] /. sol, {t, 0, 2}] Out[2]= [image] ``` #### GeneratedParameters (2) Use differently named constants: ```wl In[1]:= DSolve[y''[x] == y[x], y[x], x, GeneratedParameters -> d] Out[1]= {{y[x] -> E^x d[1] + E^-x d[2]}} ``` Use subscripted constants: ```wl In[2]:= DSolve[y''[x] == y[x], y[x], x, GeneratedParameters -> (Subscript[c, #]&)] Out[2]= {{y[x] -> E^x Subscript[c, 1] + E^-x Subscript[c, 2]}} ``` --- Generate uniquely named constants of integration: ```wl In[1]:= DSolve[y''[x] == y[x], y[x], x, GeneratedParameters -> Unique[C]] Out[1]= {{y[x] -> E^x C$751514[1] + E^-x C$751514[2]}} ``` The constants of integration are unique across different invocations of ``DSolve`` : ```wl In[2]:= DSolve[y''[x] == y[x], y[x], x, GeneratedParameters -> Unique[C]] Out[2]= {{y[x] -> E^x C$751520[1] + E^-x C$751520[2]}} ``` #### Method (1) Solve a linear ordinary differential equation: ```wl In[1]:= DSolve[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1}, y[x], x] Out[1]= {{y[x] -> Sin[x]}} ``` Obtain a solution in terms of ``DifferentialRoot`` : ```wl In[2]:= DSolve[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1}, y[x], x, Method -> "Holonomic"] Out[2]= {{y[x] -> DifferentialRoot[Function[{\[FormalY], \[FormalX]}, {\[FormalY][\[FormalX]] + Derivative[2][\[FormalY]][\[FormalX]] == 0, \[FormalY][0] == 0, Derivative[1][\[FormalY]][0] == 1}]][x]}} ``` #### SingularSolutions (1) By default, DSolve returns the general solution for this ODE: ```wl In[1]:= DSolve[y'[x] ^ 2 - y'[x] * x + y[x] == 0, y[x], x] Out[1]= {{y[x] -> x C[1] - C[1]^2}} ``` Use ``IncludeSingularSolutions`` to compute singular solutions along with the general solution: ```wl In[2]:= DSolve[y'[x] ^ 2 - y'[x] * x + y[x] == 0, y[x], x, IncludeSingularSolutions -> True] Out[2]= {{y[x] -> x C[1] - C[1]^2}, {y[x] -> (x^2/4)}} ``` ### Applications (36) #### Ordinary Differential Equations (12) Solve a logistic (Riccati) equation: ```wl In[1]:= DSolve[{y'[x] == y[x](1 - y[x] / 27), y[0] == a}, y, x]//Quiet Out[1]= {{y -> Function[{x}, (27 a E^x/27 - a + a E^x)]}} ``` Plot the solution for different initial values: ```wl In[2]:= Plot[Evaluate[y[x] /. % /. {{a -> 1 / 13}, {a -> 1 / 2}, {a -> 4}}], {x, 0, 18}] Out[2]= [image] ``` --- Solve a linear pendulum equation: ```wl In[1]:= DSolve[{y''[x] + y[x] == 0, y[0] == 1, y'[0] == 1 / 3}, y , x] Out[1]= {{y -> Function[{x}, (1/3) (3 Cos[x] + Sin[x])]}} In[2]:= Plot[y[x] /. %, {x, 0, 17}] Out[2]= [image] ``` --- Displacement of a linear, damped pendulum: ```wl In[1]:= DSolve[{y''[x] + 3y'[x] + 40 y[x] == 0, y[0] == 1, y'[0] == 1 / 3}, y, x] Out[1]= {{y -> Function[{x}, (1/453) E^-3 x / 2 (453 Cos[(Sqrt[151] x/2)] + 11 Sqrt[151] Sin[(Sqrt[151] x/2)])]}} In[2]:= Plot[y[x] /. %, {x, 0, 4}, PlotRange -> All] Out[2]= [image] ``` --- Study the phase portrait of a dynamical system: ```wl In[1]:= DSolve[{x'[t] == -2y[t] + 3x[t], y'[t] == 15x[t] - y[t], x[0] == a, y[0] == b}, {x, y}, t]; In[2]:= ParametricPlot[Evaluate[Table[{x[t], y[t]} /. % /. {a -> 1 / (13 + m), b -> 1 / (15 + m)}, {m, 0, 20, 7}]], {t, -2, 2}] Out[2]= [image] ``` --- Find a power series solution when the exact solution is known: ```wl In[1]:= DSolve[{y'[x] + Exp[x]y[x] == 1, y[0] == 3}, y , x] Out[1]= {{y -> Function[{x}, E^-E^x (3 E - ExpIntegralEi[1] + ExpIntegralEi[E^x])]}} In[2]:= Series[y[x] /. %[[1]], {x, 0, 7}] Out[2]= SeriesData[x, 0, {3, -2, Rational[-1, 2], Rational[1, 3], Rational[1, 6], Rational[-1, 120], Rational[-11, 360], Rational[-1, 105]}, 0, 8, 1] ``` --- Verify that ``y = 0`` and ``y = 1`` are equilibrium solutions of the logistic equation: ```wl In[1]:= logisticeqn = y'[x] == r y[x](1 - y[x]); In[2]:= DSolve[{logisticeqn, y[0] == 0}, y[x], x] Out[2]= {{y[x] -> 0}} In[3]:= DSolve[{logisticeqn, y[0] == 1}, y[x], x] Out[3]= {{y[x] -> 1}} ``` --- Model a block on a moving conveyor belt anchored to a wall by a spring. Compare positions and velocities for the different values of the parameters of the system (mass of the block, belt speed, friction coefficient, spring constant): [image] ```wl In[1]:= beltv[t_] = vb; spring[x_] = k(l - x); friction[v_] := -a(v - beltv[t]); ``` Newton's equation for the block: ```wl In[2]:= sys := {m x''[t] == spring[x[t]] + friction[x'[t]], x[0] == l, x'[0] == 0}; ``` Solve for position and velocity: ```wl In[3]:= sol = DSolve[sys, x, t]; In[4]:= {pos[m_, k_, a_, vb_, t_], vel[m_, k_, a_, vb_, t_]} = {x[t], x'[t]} /. sol[[1]]//FullSimplify Out[4]= {(1/2 k Sqrt[a^2 - 4 k m])E^-((a + Sqrt[a^2 - 4 k m]) t/2 m) (a (a - Sqrt[a^2 - 4 k m]) vb - a E^(Sqrt[a^2 - 4 k m] t/m) (a + Sqrt[a^2 - 4 k m]) vb + 2 E^((a + Sqrt[a^2 - 4 k m]) t/2 m) Sqrt[a^2 - 4 k m] (k l + a vb)), (2 a E^-(a t/2 m) vb Sinh[(Sqrt[a^2 - 4 k m] t/2 m)]/Sqrt[a^2 - 4 k m])} ``` The block stabilizes just above the spring's natural length of ``l`` : ```wl In[5]:= l = 1; Manipulate[{Plot[{pos[m, k, a, vb, t]}, {t, 0, 2}, PlotRange -> All, PlotLabel -> "Position"], Plot[{vel[m, k, a, vb, t], vb}, {t, 0, 2}, PlotRange -> All, PlotLabel -> "Velocity"]}, {{m, 1, "mass of the block"}, 1, 10, 1}, {{k, 100, "spring constant"}, 100, 1000, 100}, {{vb, 1, "belt velocity"}, 1, 10, 1}, {{a, 5, "friction coefficient"}, 5, 45, 10}, SaveDefinitions -> True] Out[5]= DynamicModule[«9»] ``` --- Use component laws together with Kirchhoff's laws for connections to simulate the response of an RLC circuit to a step impulse in the voltage $v_1$ applied at time $t=0.01$ : [image] The component laws read $R I_R(t)=V_R(t)$, $L I_L'(t)=V_L(t)$ $I_C(t)=C V_C'(t)$, $V_C(t)+V_L(t)+V_R(t)-V_1(t)=0$ and $I_L(t)=I_R(t)=I_C(t)$ : ```wl In[1]:= ClearAll[r, l, c]; vR[t_] := r iL[t]; vL[t_] := v1[t] - vR[t] - vC[t]; sys = {l iL'[t] == vL[t], c vC'[t] == iL[t]}; ``` Simulate a step response: ```wl In[2]:= v1[t_] := UnitStep[t - .01]; ic = {iL[0] == 0, vC[0] == 0}; In[3]:= sol = DSolve[Flatten[{sys, ic}], {iL, vC}, t]; In[4]:= {iL1[r_, c_, l_, t_], vC1[r_, c_, l_, t_]} = {iL[t], vC[t]} /. sol[[1]]; ``` Visualize the response for different values of $R$, $L$ and $C$ : In[5]:= Manipulate[{Plot[Evaluate[{v1[t],Re[vC1[r,c,l,t]]}],{t,0,.03},Filling->{1->0},Ticks->{Range[0.01,0.03,0.01],Automatic},PlotRange->All,PlotLabel->"Subscript[V, C][t]"],Plot[Evaluate[Re[iL1[r,c,l,t]]],{t,0,.03},Ticks->{Range[0.01,0.03,0.01],Automatic},PlotRange->All,PlotLabel->"Subscript[I, L][t]"]},{{r,1,"R"},1,10,1},{{l,10^-2,"L"},10^-2,5 10^-2,10^-2},{{c,10^-4,"C"},10^-4,5 10^-4,10^-4},SaveDefinitions->True] Out[5]= Manipulate[{Plot[Evaluate[{v1[t], Re[vC1[r, c, l, t]]}], {t, 0, 0.03}, Filling -> {1 -> 0}, Ticks -> {Range[0.01, 0.03, 0.01], Automatic}, PlotRange -> All, PlotLabel -> "\!\(\*SubscriptBox[\(V\), \(C\)]\)[t]"], Plot[Evaluate[Re[iL1[r ... (Sqrt[c]\*l))\*E^((0.005\*Sqrt[-4.\*l + (1.\*c)\*r^2])/(Sqrt[c]\*l)))\* (-4.\*l + (1.\*c)\*r^2))\*((1.\*Sqrt[c])\*r - 1.\*Sqrt[-4.\*l + (1.\*c)\*r^2]))\*((1.\*Sqrt[c])\*r + 1.\*Sqrt[-4.\*l + (1.\*c)\*r^2]))}] --- Model the change in height of water in two cylindrical tanks as water flows from one tank to another through a pipe: [image] Use pressure relations $p=\rho g h$ and mass conservation: ```wl In[1]:= p1[t_] := ρ g h1[t]; p2[t_] := ρ g h2[t]; ``` Model the flow across the pipe with the Hagen–Poiseuille relation: ```wl In[2]:= flowRate = (p1[t] - p2[t])(π pipeDia ^ 4) / (128 μ pipeLen); massConservation = {a1 h1'[t] == -flowRate, a2 h2'[t] == flowRate}; ic = {h1[0] == 1, h2[0] == 0}; ``` Simulate the system: ```wl In[3]:= sol = DSolve[Flatten[{massConservation, ic}], {h1, h2}, t] Out[3]= {{h1 -> Function[{t}, (a1 + a2 E^((-a1 - a2) g π pipeDia^4 t ρ/128 a1 a2 pipeLen μ)/a1 + a2)], h2 -> Function[{t}, -(a1 (-1 + E^((-a1 - a2) g π pipeDia^4 t ρ/128 a1 a2 pipeLen μ))/a1 + a2)]}} In[4]:= {h1[t_], h2[t_]} = {h1[t], h2[t]} /. sol[[1]]; ``` Visualize the solution for particular values of the parameters: ```wl In[5]:= params = {pipeLen -> 0.1, pipeDia -> 0.2, a1 -> 1, a2 -> 1, ρ -> 0.2, μ -> 2 10 ^ -3, g -> 9.81}; ``` In[6]:= Plot[Evaluate[{h1[t],h2[t]}/.params],{t,0,15},PlotRange->All,PlotLegends->{"Subscript[h, 1][t]","Subscript[h, 2][t]"}] Out[6]= Subscript[h, 1][t] Subscript[h, 2][t] --- Solve the equation of a fractional harmonic oscillator of order 1.9: ```wl In[1]:= DSolve[{CaputoD[x[t], {t, 1.9}] + x[t] == 0, x[0] == 1, x'[0] == 0}, x[t], t] Out[1]= {{x[t] -> 1. MittagLefflerE[1.9, -1. t^1.9]}} ``` Plot the solution: ```wl In[2]:= fracsol = Plot[x[t] /. %, {t, 0, 20}] Out[2]= [image] ``` A fractional harmonic oscillator behaves like a damped harmonic oscillator: ```wl In[3]:= dampedOsc = DSolve[{x''[t] + a x'[t] + x[t] == 0, x[0] == 1, x'[0] == 0}, x[t], t]//FullSimplify Out[3]= {{x[t] -> E^-(a t/2) (Cosh[(1/2) Sqrt[-4 + a^2] t] + (a Sinh[(1/2) Sqrt[-4 + a^2] t]/Sqrt[-4 + a^2]))}} ``` Fit the data for a fractional oscillator with the solution of a damped oscillator to determine the corresponding value for damping factor $a$ : ```wl In[4]:= data = Cases[InputForm[fracsol], Line[m_] :> m, -1][[1]]; fit = FindFit[data, x[t] /. dampedOsc, a, t] Out[4]= {a -> 0.152103} ``` Compare the results: ```wl In[5]:= dampedsol = Plot[x[t] /. dampedOsc /. fit, {t, 0, 20}]; Show[fracsol, dampedsol] Out[5]= [image] ``` --- Solve the equation of fractional LC electric circuit: ```wl In[1]:= LCeqn[α_, l_, c_] = {CaputoD[i[t], {t, α}] + (1/l c)i[t] == 0, i[0] == 0.1, i'[0] == 0}; In[2]:= LCsol[α_, l_, c_] := DSolve[LCeqn[α, l, c], i, t]; ``` The classical solution can be obtained for $\alpha =2$ : ```wl In[3]:= LCsol[2, l, c] Out[3]= {{i -> Function[{t}, 0.1 Cos[(t/Sqrt[c] Sqrt[l])]]}} ``` Solve the LC circuit equation for $\alpha =1.9$ : ```wl In[4]:= LCsol[1.9, l, c] Out[4]= {{i -> Function[{t}, 0.1 MittagLefflerE[1.9, -(t^1.9/c l)]]}} ``` Plot the solutions for various values of equation order $\alpha$ : ```wl In[5]:= Plot[Evaluate[i[t] /. {LCsol[2, 1, 1], LCsol[1.9, 1, 1], LCsol[1.8, 1, 1]}], {t, 0, 10}, PlotLegends -> {"α=2", "α=1.9", "α=1.8"}] Out[5]= [image] ``` --- Solve the equation of fractional RC electric circuit: ```wl In[1]:= RCeqn[α_, r_, c_] = {CaputoD[v[t], {t, α}] + (1/r c)v[t] == 0, v[0] == 20}; In[2]:= RCsol[α_, r_, c_] := DSolve[RCeqn[α, r, c], v, t] ``` The classical solution can be obtained for $\alpha =1$ : ```wl In[3]:= RCsol[1, r, c] Out[3]= {{v -> Function[{t}, 20 E^-(t/c r)]}} ``` Solve the RC circuit equation for $\alpha =0.9$ : ```wl In[4]:= RCsol[0.9, r, c] Out[4]= {{v -> Function[{t}, 20. MittagLefflerE[0.9, -(t^0.9/c r)]]}} ``` Plot the solutions for various values of equation order $\alpha$ : ```wl In[5]:= Plot[Evaluate[v[t] /. {RCsol[1, 10, 1], RCsol[0.9, 10, 1], RCsol[0.8, 10, 1]}], {t, 0, 100}, PlotLegends -> {"α=1", "α=0.9", "α=0.8"}] Out[5]= [image] ``` #### Hybrid Differential Equations (8) Model a damped oscillator that gets a kick at regular time intervals: ```wl In[1]:= system = {x''[t] + .1x'[t] + x[t] == 0, x[0] == 1, x'[0] == 0}; control = WhenEvent[Mod[t, 1] == 0, x'[t] -> x'[t] + 1]; In[2]:= sol = DSolve[{system, control}, x, {t, 0, 50}]; In[3]:= Plot[Evaluate[{x[t], x'[t]} /. sol], {t, 0, 50}] Out[3]= [image] ``` --- Rectify a sine wave: ```wl In[1]:= sol = DSolve[{x''[t] + x[t] == 0, x[0] == 0, x'[0] == 1, WhenEvent[x[t] == 0, x'[t] -> -x'[t]]}, x, {t, 0, 4π}] Out[1]= {{x -> Function[{t}, Piecewise[{{Sin[t], 0 <= t <= Pi}, {-Sin[t], Inequality[Pi, Less, t, LessEqual, 2*Pi]}, {Sin[t], Inequality[2*Pi, Less, t, LessEqual, 3*Pi]}, {-Sin[t], Inequality[3*Pi, Less, t, LessEqual, 4*Pi]}}, Indeterminate]]}} In[2]:= Plot[{Sin[t], (x[t] /. sol)}, {t, 0, 4π}] Out[2]= [image] ``` --- Model a ball bouncing down steps: ```wl In[1]:= c = .75; sol = DSolve[{y''[t] == -9.8, y[0] == 13.5, y'[0] == 5, a[0] == 13, WhenEvent[y[t] - a[t] == 0, y'[t] -> -c y'[t]], WhenEvent[Mod[t, 1], a[t] -> a[t] - 1]}, {y, a}, {t, 0, 8}, DiscreteVariables -> {a}] ; In[2]:= Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 8}, Filling -> {2 -> 0}, Exclusions -> None] Out[2]= [image] ``` --- In a square box, model a ball that changes direction upon impact with the side walls: ```wl In[1]:= sol = DSolve[{x'[t] == a[t], y'[t] == b[t], x[0] == 0, y[0] == 0, a[0] == 1, b[0] == 17 / 12, WhenEvent[x[t] ^ 2 == 1, a[t] -> -a[t]], WhenEvent[y[t] ^ 2 == 1, b[t] -> -b[t]]}, {x, y}, {t, 0, 50}, DiscreteVariables -> {a, b}]; In[2]:= ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 50}, Frame -> True, FrameTicks -> None, PlotRange -> 1, Axes -> False] Out[2]= [image] ``` --- Simulate the system $y'=y+u$ stabilized with a discrete-time controller $u(k \tau )=-2 y(k \tau )$ : ```wl In[1]:= system = {y'[t] == y[t] + u[t], y[0] == 1, u[0] == 0}; control = WhenEvent[Mod[t, 1] == 0, u[t] -> -2y[t]]; ``` Simulate and visualize: ```wl In[2]:= DSolve[{system, control}, {y, u}, {t, 0, 20}, DiscreteVariables -> u]; In[3]:= Plot[Evaluate[{u[t], y[t]} /. %], {t, 0, 10}, PlotRange -> {-6, 6}, Exclusions -> None] Out[3]= [image] ``` --- Control a double integrator using a dead-beat, discrete-time controller: ```wl In[1]:= system = {x''[t] == u[t], x[0] == 1, x'[0] == 1, u[0] == -7}; ``` Use a dead-beat digital feedback controller $u(k \tau )=-1\left/\tau ^2\right.x(k \tau )-3/(2 \tau ) x'(k \tau )$ : ```wl In[2]:= τ = 0.5; control = WhenEvent[Mod[t, τ] == 0, u[t] -> -(1 / τ ^ 2)x[t] - (3 / (2τ))x'[t]]; ``` Simulate and visualize: ```wl In[3]:= DSolve[{system, control}, {x, u}, {t, 0, 2}, DiscreteVariables -> u]; In[4]:= Plot[Evaluate[{x[t], x'[t], u[t]} /. %], {t, 0, 2}, PlotRange -> All, Exclusions -> None] Out[4]= [image] ``` --- Model the position of a moving body with 1 kg mass: ```wl In[1]:= system = {x''[t] == u[t], x'[0] == x[0] == 0, u[0] == 1}; ``` Use a sampled proportional-derivative (PD) controller to keep the position constant: ```wl In[2]:= kp = 1; td = 1;τ = 1; xref = 1; control = WhenEvent[Mod[t, τ] == 0, u[t] -> kp(xref - x[t] - td x'[t])]; ``` Simulate and visualize: ```wl In[3]:= DSolve[{system, control}, {x, u}, {t, 0, 12}, DiscreteVariables -> u]; In[4]:= Plot[Evaluate[{xref, x[t], u[t]} /. %], {t, 0, 12}, Exclusions -> None] Out[4]= [image] ``` --- Set the physical variable ``a`` to ``0`` whenever it becomes negative: ```wl In[1]:= sol = DSolve[{a'[t] == -1 / 11u[t], a[0] == 1, u[0] == 1, WhenEvent[a[t] < 0, u[t] -> 0]}, {a, u}, {t, 0, 15}, DiscreteVariables -> {u}] Out[1]= {{a -> Function[{t}, Piecewise[{{1 - t/11, 0 <= t <= 11}, {0, Inequality[11, Less, t, LessEqual, 15]}}, Indeterminate]], u -> Function[{t}, Piecewise[{{1, 0 <= t <= 11}, {0, Inequality[11, Less, t, LessEqual, 15]}}, Indeterminate]]}} In[2]:= Plot[Evaluate[a[t] /. sol], {t, 0, 15}, PlotRange -> All] Out[2]= [image] ``` #### Delay Differential Equations (2) Study the onset of chaos in a dynamical system governed by a delay differential equation: ```wl In[1]:= sol = DSolve[{x'[t] == Sign[x[t - 2]] - x[t - 2], x[t /; t <= 0] == 0.1}, x, {t, 0, 30}]; In[2]:= ParametricPlot[Evaluate[{x[t - 2], x[t]} /. sol], {t, 2, 30}] Out[2]= [image] ``` --- Study the stability of the solutions for a linear delay differential equation: ```wl In[1]:= f[λ_, μ_] := x[t] /. DSolve[{x'[t] == λ x[t] + μ x[t - 1], x[t /; t ≤ 0] == 1 - t}, x[t], {t, 0, 10}][[1]] In[2]:= Plot[Evaluate[{f[1 / 2, -1], f[-7 / 2, 4], f[-5, 4]}], {t, 0, 15}, Exclusions -> None] Out[2]= [image] ``` #### Integral Equations (4) The tautochrone problem requires finding the curve down which a bead placed anywhere will fall to the bottom in the same amount of time. Expressing the total fall time in terms of the arc length of the curve and the speed ``v`` yields the Abel integral equation $T=\int dt=\int 1/v \, ds$. Defining the unknown function $h$ by the relationship $ds=h(y) dy$ and using the conservation of energy equation $\left.v^2\right/2=g \left(y_{\max }-y\right)$ yields the explicit equation: ```wl In[1]:= abeleqn = T == Subsuperscript[∫, 0, y](h[z]/Sqrt[2g]Sqrt[y - z])\[DifferentialD]z; ``` Use ``DSolve`` to solve the integral equation: ```wl In[2]:= DSolve[abeleqn, h[y], y] Out[2]= {{h[y] -> (Sqrt[2] g T/π Sqrt[g y])}} ``` Using the relationship $ds=\sqrt{dx^2+dy^2}$, solve for $x'(y)$ : ```wl In[3]:= dxdy = Sqrt[h[y]^2 - 1] /. First[%] Out[3]= Sqrt[-1 + (2 g T^2/π^2 y)] ``` Starting the curve from the origin and integrating—with assumptions that ensure the integrand is real-valued—yields $x$ as a function of $y$ : ```wl In[4]:= x[y_] = Integrate[dxdy, {y, 0, y}, Assumptions -> (2 g T^2 /π^2y) > 1 && y > 0] Out[4]= (π Sqrt[y (2 g T^2 - π^2 y)] + 2 g T^2 ArcTan[(π/Sqrt[-π^2 + (2 g T^2/y)])]/π^2) ``` Substituting values for $g$ and $T$, use ``ParametricPlot`` to display the maximal tautochrone curve: ```wl In[5]:= Show[ParametricPlot[{{x[y], y}, {-x[y], y}} /. {g -> 9.8, T -> 2}, {y, 0, (2(9.8)2^2/π^2)}], ImageSize -> Medium] Out[5]= [image] ``` Making the change of variables $y=\frac{g T^2 (1-\cos (\theta ))}{\pi ^2}$ gives a simple, nonsingular parametrization of the curve with $-\pi <\theta <\pi$ : ```wl In[6]:= c[θ_] = (g T^2/π^2){Sin[θ] + θ, 1 - Cos[θ]} ; ``` Combining the conservation of energy equation and the chain rule $\frac{ds}{dt}=\frac{ds}{d\theta }\frac{d\theta }{dt}$ produces the following differential equation for $\theta$ as a function of $t$ : ```wl In[7]:= FullSimplify[ (Sqrt[2g (Last[c[θMax]] - Last[c[θ]])] /Sqrt[c'[θ].c'[θ]]) , g > 0 && T > 0] Out[7]= (π Sqrt[Cos[θ] - Cos[θMax]]/T Sqrt[1 + Cos[θ]]) In[8]:= teqn = θ'[t] == (% /. {θ -> θ[t], T -> 2}); ``` This equation can be solved numerically for different starting points: ```wl In[9]:= Θ[t_] = Table[NDSolveValue[{teqn, θ[0] == θMax}, θ[t], {t, 0, 4}], {θMax, {-2 Pi / 3, -Pi / 2, -Pi / 4, -Pi / 12}}]; ``` Plotting the solutions shows that all reach the bottom, $\theta =0$, at the two-second mark: ```wl In[10]:= Plot[Evaluate[Θ[t]], {t, 0, 4}, PlotTheme -> "Scientific"] Out[10]= [image] ``` Visualize the motion along the tautochrone: ```wl In[11]:= {cShape[t_], cBeads[t_]} = {c[t], c /@ Θ[t]} /. {T -> 2, g -> 9.8}; Animate@@{ParametricPlot[cShape[s] , {s, -Pi, Pi}, ImageSize -> Medium, Epilog -> {AbsolutePointSize[8], Point[cBeads[t]]}], {t, 0, 4, Appearance -> "Labeled"}, SaveDefinitions -> True, AnimationDirection -> ForwardBackward} Out[11]= DynamicModule[«8»] ``` --- A spring-mass system, with an attached mass of $$3g$$, has a spring constant of $$3\text{dynes}/\text{cm}$$ and a damping coefficient of $$2g/s$$. At time $t=0$, the mass is pushed down and then released with a velocity 28 cm/s downward. A force of $50 \cos (2 t)\text{dynes}$ acts downward on the mass for $t\geq 0$. Find the velocity of the system as a function of time. The integral equation for the velocity is given by the following: ```wl In[1]:= ieqn = 3Derivative[1][v][t] + 2 v[t] + 40 Subsuperscript[∫, 0, t]v[m]\[DifferentialD]m == 50 Cos[2 t]; ``` Initial value of the velocity: ```wl In[2]:= ic = v[0] == 28; ``` Solve the integro-differential equation using ``DSolve`` : ```wl In[3]:= sol = DSolve[{ieqn, ic}, v[t], t] Out[3]= {{v[t] -> (1/14) (7 Cos[2 t] + 385 E^-t / 3 Cos[(Sqrt[119] t/3)] - 49 Sin[2 t] + 5 Sqrt[119] E^-t / 3 Sin[(Sqrt[119] t/3)])}} ``` Plot the velocity as a function of time: ```wl In[4]:= Plot[v[t] /. sol[[1]], {t, 0, 8}, PlotRange -> All, Filling -> Axis] Out[4]= [image] ``` Visualize the motion of the physical spring-mass system: [image] --- An LRC circuit has a voltage source given by $v(t)=2 \cos (t)V$. The resistance in the circuit is 0.2 $\Omega$, the inductance is $$5H$$, and the capacitance is $$0.7F$$. Initially, the current in the resistor is $$1A$$. Find the current as a function of time: [image] The integral equation for the current is given by: ```wl In[1]:= ieqn = 5 Derivative[1][i][t] + 0.2 i[t] + (1/0.7)Subsuperscript[∫, 0, t]i[m]\[DifferentialD]m == 2 Cos[t]; ``` Initial value of the current: ```wl In[2]:= ic = i[0] == 1; ``` Solve the integro-differential equation using ``DSolve`` : ```wl In[3]:= sol = DSolve[{ieqn, ic}, i[t], t]//FullSimplify//Chop Out[3]= {{i[t] -> 0.031262 Cos[1. t] + E^-0.02 t (0.968738 Cos[0.534148 t] - 0.334878 Sin[0.534148 t]) + 0.558249 Sin[1. t]}} ``` Plot the current as a function of time: ```wl In[4]:= Plot[i[t] /. sol[[1]], {t, 0, 80}, PlotRange -> All, Filling -> Axis] Out[4]= [image] ``` --- A linear Volterra integral equation is equivalent to an initial value problem for a linear differential equation. Verify this relationship for the following Volterra equation: ```wl In[1]:= eqn = y[x] == x^3 + λSubsuperscript[∫, 0, x](t - x)y[t]\[DifferentialD]t; ``` Solve the integral equation using ``DSolve`` : ```wl In[2]:= sol = DSolve[eqn, y[x], x] Out[2]= {{y[x] -> (6 x/λ) - (6 Sin[x Sqrt[λ]]/λ^3 / 2)}} ``` Set up the corresponding differential equation: ```wl In[3]:= deqn = D[eqn, {x, 2}] Out[3]= Derivative[2][y][x] == 6 x - λ y[x] ``` Add two initial conditions since the differential equation is of second order: ```wl In[4]:= init = {(eqn /. {x -> 0}), (D[eqn, x] /. {x -> 0})} Out[4]= {y[0] == 0, Derivative[1][y][0] == 0} ``` The solution of the initial value problem agrees with that of the integral equation: ```wl In[5]:= DSolve[{deqn, init}, y[x], x]//Simplify Out[5]= {{y[x] -> (6 x/λ) - (6 Sin[x Sqrt[λ]]/λ^3 / 2)}} ``` Plot the solution: ```wl In[6]:= Plot[Table[y[x] /. %[[1]], {λ, 1, 3, 0.5}]//Evaluate, {x, 0, 20}, Filling -> Axis] Out[6]= [image] ``` #### Classical Partial Differential Equations (5) Model the vibrations of a string with fixed length, say ``π``, using the wave equation: ```wl In[1]:= eqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}]; ``` Specify that the ends of the string remain fixed during the vibrations: ```wl In[2]:= bc = {u[0, t] == 0, u[π, t] == 0}; ``` Obtain the fundamental and higher harmonic modes of oscillation: ```wl In[3]:= ic = {u[x, 0] == Sin[m x], Derivative[0, 1][u][x, 0] == 0}; In[4]:= dsol = DSolve[{eqn, bc, ic}, u, {x, t}, Assumptions -> Element[m, Integers]] Out[4]= {{u -> Function[{x, t}, Cos[m t] Sin[m x]]}} ``` Visualize the vibrations of the string for these modes: ```wl In[5]:= sol1[x_, t_] = u[x, t] /. dsol[[1]] Out[5]= Cos[m t] Sin[m x] In[6]:= Table[Show[Plot[Table[sol1[x, t], {t, 0, 4}]//Evaluate, {x, 0, Pi}, Ticks -> False]], {m, 4}] Out[6]= {[image], [image], [image], [image]} ``` In general, the solution is composed of an infinite number of harmonics: ```wl In[7]:= ic = {u[x, 0] == x(π - x) ^ 2, Derivative[0, 1][u][x, 0] == 0}; In[8]:= dsol = DSolve[{eqn, bc, ic}, u, {x, t}] /. {K[1] -> m} Out[8]= {{u -> Function[{x, t}, Inactive[Sum][(4*(2 + (-1)^m)*Cos[t*m]*Sin[x*m])/m^3, {m, 1, Infinity}]]}} ``` Extract four terms from the ``Inactive`` sum: ```wl In[9]:= sol[x_, t_] = TruncateSum[u[x, t] /. dsol[[1]], 4] Out[9]= 4 Cos[t] Sin[x] + (3/2) Cos[2 t] Sin[2 x] + (4/27) Cos[3 t] Sin[3 x] + (3/16) Cos[4 t] Sin[4 x] ``` Visualize the vibration of the string: ```wl In[10]:= Animate[Plot[sol[x, t], {x, 0, π}, PerformanceGoal -> "Quality", PlotRange -> {-5, 5}, ImageSize -> Medium], {t, 0, 12}, SaveDefinitions -> True] Out[10]= DynamicModule[«9»] ``` --- Model the oscillations of a circular membrane of radius ``1`` using the wave equation in 2D: ```wl In[1]:= eqn = r D[u[r, t], {t, 2}] == D[r D[u[r, t], r], r]; ``` Specify that the boundary of the membrane remain fixed: ```wl In[2]:= bc = u[1, t] == 0; ``` Initial condition for the problem: ```wl In[3]:= ic = {u[r, 0] == 0, Derivative[0, 1][u][r, 0] == 1}; ``` Obtain a solution in terms of Bessel functions: ```wl In[4]:= (dsol = DSolve[{eqn, bc, ic}, u[r, t], {r, t}])//TraditionalForm Out[4]//TraditionalForm= $$\left\{\left\{u(r,t)\to \underset{K[1]=1}{\overset{\infty }{\sum }}\frac{2 J_0\left(r j_{0,K[1]}\right) J_1\left(j_{0,K[1]}\right) \sin \left(t j_{0,K[1]}\right)}{\left(J_0\left(j_{0,K[1]}\right){}^2+J_1\left(j_{0,K[1]}\right){}^2\right) \left(j_{0,K[1]}\right){}^2}\right\}\right\}$$ ``` Extract a finite number of terms from the ``Inactive`` sum: ```wl In[5]:= h[r_, t_] = TruncateSum[u[r, t] /. dsol[[1]], 3]//N; ``` Visualize the oscillations of the membrane: ```wl In[6]:= Animate[Plot3D[h[r, t] /. {r -> Sqrt[x ^ 2 + y ^ 2]}, {x, y}∈Disk[], PlotRange -> {-1, 1}, Ticks -> None, Mesh -> True, MeshStyle -> {Red, Blue}, PlotStyle -> Yellow], {t, 0, 4}, SaveDefinitions -> True] Out[6]= DynamicModule[«9»] ``` --- Model the flow of heat in a bar of length ``1`` using the heat equation: ```wl In[1]:= heqn = D[u[x, t], t] == D[u[x, t], {x, 2}]; ``` Specify the fixed temperature at both ends of the bar: ```wl In[2]:= bc = {u[0, t] == 20, u[1, t] == 50}; ``` Specify an initial condition: ```wl In[3]:= ic = u[x, 0] == 0; ``` Solve the heat equation subject to these conditions: ```wl In[4]:= sol = DSolve[{heqn, bc, ic}, u[x, t], {x, t}] Out[4]= {{u[x, t] -> 20 + 30 x - (2 Inactive[Sum][((20 - 50*(-1)^K[1])*Sin[Pi*x*K[1]])/(E^(Pi^2*t*K[1]^2)*K[1]), {K[1], 1, Infinity}]/π)}} ``` Extract a few terms from the ``Inactive`` sum: ```wl In[5]:= approxsol = TruncateSum[u[x, t] /. sol[[1]], 3]//Expand Out[5]= 20 + 30 x - (140 E^-π^2 t Sin[π x]/π) + (30 E^-4 π^2 t Sin[2 π x]/π) - (140 E^-9 π^2 t Sin[3 π x]/3 π) ``` Visualize the evolution of the temperature to a steady state: ```wl In[6]:= Plot[Table[approxsol, {t, 0.02, 0.5, 0.07}]//Evaluate, {x, 0, 1}, AxesOrigin -> {0, 0}] Out[6]= [image] ``` Obtain the steady-state solution ``v[x]``, which is independent of time: ```wl In[7]:= ssol = DSolve[{v''[x] == 0, v[0] == 20, v[1] == 50}, v[x], x]//Expand Out[7]= {{v[x] -> 20 + 30 x}} ``` The steady-state solution is a linear function of ``x`` : ```wl In[8]:= Plot[v[x] /. ssol[[1]], {x, 0, 1}] Out[8]= [image] ``` --- Model the flow of heat in a bar of length ``1`` that is insulated at both ends: ```wl In[1]:= heqn = D[u[x, t], t] == D[u[x, t], {x, 2}]; ``` Specify that no heat flow through the ends of the bar: ```wl In[2]:= bc = {Derivative[1, 0][u][0, t] == 0, Derivative[1, 0][u][1, t] == 0}; ``` Specify an initial condition: ```wl In[3]:= ic = u[x, 0] == 20 + 80x; ``` Solve the heat equation subject to these conditions: ```wl In[4]:= sol = DSolve[{heqn, bc, ic}, u[x, t], {x, t}] Out[4]= {{u[x, t] -> 60 + 2 Inactive[Sum][(80*(-1 + (-1)^K[1])*Cos[Pi*x*K[1]])/(E^(Pi^2*t*K[1]^2)*(Pi^2*K[1]^2)), {K[1], 1, Infinity}]}} ``` Extract a few terms from the ``Inactive`` sum: ```wl In[5]:= approxsol = TruncateSum[u[x, t] /. sol[[1]], 4]//Expand Out[5]= 60 - (320 E^-π^2 t Cos[π x]/π^2) - (320 E^-9 π^2 t Cos[3 π x]/9 π^2) ``` Visualize the evolution of the temperature to the steady-state value of ``60°`` : ```wl In[6]:= Plot[Table[approxsol, {t, 0.02, 0.9, 0.07}]//Evaluate, {x, 0, 1}, AxesOrigin -> {0, 0}, PlotRange -> All] Out[6]= [image] ``` --- Construct a complex analytic function, starting from the values of its real and imaginary parts on the $x$ axis. The real and imaginary parts ``u`` and ``v`` satisfy the Cauchy–Riemann equations: ```wl In[1]:= creqns = {D[u[x, y], x] == D[v[x, y], y], D[v[x, y], x] == -D[u[x, y], y]}; ``` Prescribe the values of ``u`` and ``v`` on the $x$ axis: ```wl In[2]:= xvals = {u[x, 0] == x ^ 3, v[x, 0] == 0}; ``` Solve the Cauchy–Riemann equations: ```wl In[3]:= sol = DSolve[{creqns, xvals}, {u, v}, {x, y}] Out[3]= {{u -> Function[{x, y}, x^3 - 3 x y^2], v -> Function[{x, y}, 3 x^2 y - y^3]}} ``` Verify that the solutions are harmonic functions: ```wl In[4]:= Laplacian[{u[x, y], v[x, y]} /. sol[[1]], {x, y}] Out[4]= {0, 0} ``` Visualize the streamlines and equipotentials generated by the solution: ```wl In[5]:= ContourPlot[{u[x, y], v[x, y]} /. sol[[1]], {x, -5, 5}, {y, -5, 5}, ContourStyle -> {Red, Blue}, PlotTheme -> "Marketing"] Out[5]= [image] ``` Construct an analytic function from the solution: ```wl In[6]:= f[x_, y_] = u[x, y] + I v[x, y] /. sol[[1]] Out[6]= x^3 - 3 x y^2 + I (3 x^2 y - y^3) ``` This represents the function $f(z)=z^3$ : ```wl In[7]:= (f[x, y]//Factor) /. {x + I y -> z} Out[7]= z^3 ``` #### General Partial Differential Equations (5) Study the evolution of a smooth solution for Burgers' equation to a shock wave in the limit when the viscosity parameter becomes infinitely small: ```wl In[1]:= eqns = { D[u[x, t], t] + u[x, t]D[u[x, t], x] == ϵ D[u[x, t], {x, 2}] , u[x, 0] == Piecewise[{{1, x < 0}}]}; In[2]:= dsol = DSolve[eqns, u, {x, t}] Out[2]= {{u -> Function[{x, t}, (1/1 + (E^-(t - 2 x/4 ϵ) (1 + Erf[(x/2 Sqrt[t ϵ])])/1 + Erf[(t - x/2 Sqrt[t ϵ])]))]}} ``` The solution is smooth for any positive value of ``ϵ`` : ```wl In[3]:= Plot3D[u[x, t] /. dsol[[1]] /. {ϵ -> 1 / 10}, {x, -2, 2}, {t, 0.001, 5}] Out[3]= [image] ``` The solution develops a shock discontinuity in the limit when ``ϵ`` approaches ``0`` : ```wl In[4]:= Table[Plot3D[u[x, t] /. dsol[[1]] /. {ϵ -> 1 / 1000}, {x, -2, 2}, {t, 0.001, 5}, Exclusions -> All, Ticks -> None], {ϵ, {1 / 10, 1 / 100, 1 / 1000}}]//Quiet Out[4]= [image] ``` --- An electron constrained to move in a one-dimensional box of length ``d`` is governed by the free Schrödinger equation with Dirichlet conditions at the endpoints: ```wl In[1]:= eqn = I ℏ D[ψ[x, t], t] == (-ℏ^2/2m) D[ψ[x, t], {x, 2}]; bcs = {ψ[0, t] == 0, ψ[d, t] == 0}; ``` The general solution to this equation is a sum of trigonometric-exponential terms: ```wl In[2]:= DSolve[Join[{eqn}, bcs], ψ[x, t], {x, t}] Out[2]= {{ψ[x, t] -> Inactive[Sum][(C[K[1]]*Sin[(Pi*x*K[1])/d])/E^((I*Pi^2*t*ℏ*K[1]^2)/(2*d^2*m)), {K[1], 1, Infinity}]}} ``` Each term in the sum is called a stationary state, as using the sine as initial conditions leads to the positional probability density $\rho =\psi ^*\psi$ being time independent. For example: ```wl In[3]:= sol1 = DSolve[Join[{eqn, ψ[x, 0] == Sqrt[(2/d)]Sin[Pi(x/d)]}, bcs], ψ[x, t], {x, t}] Out[3]= {{ψ[x, t] -> Sqrt[2] Sqrt[(1/d)] E^-(I π^2 t ℏ/2 d^2 m) Sin[(π x/d)]}} ``` The resulting probability distribution is independent of time: ```wl In[4]:= ρ1[x_, t_] = Simplify[ComplexExpand[Conjugate[ψ[x, t]]ψ[x, t] /. First[sol1]], d > 0] Out[4]= (2 Sin[(π x/d)]^2/d) ``` The normalization of the initial data was chosen so that the integral of the density (the total probability of finding the particle somewhere) is ``1`` : ```wl In[5]:= Integrate[ρ1[x, t], {x, 0, d}] Out[5]= 1 ``` Using any other initial condition, even one as simple as a sum of two stationary states, will lead to a complicated, time-dependent density: ```wl In[6]:= sol2 = DSolve[Join[{eqn, ψ[x, 0] == (1/Sqrt[d])Sin[Pi(x/d)] + (1/Sqrt[d])Sin[2Pi(x/d)]}, bcs], ψ[x, t], {x, t}] Out[6]= {{ψ[x, t] -> (E^-(2 I π^2 t ℏ/d^2 m) (E^(3 I π^2 t ℏ/2 d^2 m) + 2 Cos[(π x/d)]) Sin[(π x/d)]/Sqrt[d])}} ``` This density is not stationary, as ``t`` appears in the second and third cosines: ```wl In[7]:= ρ2[x_, t_] = Simplify[ComplexExpand[Conjugate[ψ[x, t]]ψ[x, t] /. First[sol2]], d > 0] Out[7]= ((3 + 2 Cos[(2 π x/d)] + 2 Cos[(π (2 d m x - 3 π t ℏ)/2 d^2 m)] + 2 Cos[(π (2 d m x + 3 π t ℏ)/2 d^2 m)]) Sin[(π x/d)]^2/d) ``` Although the probability density is time dependent, its integral is still the constant ``1`` : ```wl In[8]:= Integrate[ρ2[x, t], {x, 0, d}] Out[8]= 1 ``` Entering the mass of the electron and the value of $\hbar$ in SI units and setting ``d`` to a typical interatomic distance of 1 nm results in the following density function: ```wl In[9]:= ρe[x_, t_] = ρ2[x, t] /. {ℏ -> 1.055*^-34, m -> 9.109 * 10 ^ -31, d -> 10 ^ -9} Out[9]= 1000000000 (3 + 2 Cos[1.724444315286965*^48 (-9.943140748611695`*^-34 t + 1.8218000000000004`*^-39 x)] + 2 Cos[1.724444315286965*^48 (9.943140748611695`*^-34 t + 1.8218000000000004`*^-39 x)] + 2 Cos[2000000000 π x]) Sin[1000000000 π x]^2 ``` Visualize the function over the spatial domain and one period in time: In[10]:= Plot3D[\[Rho]e[x 10^-9,t 10^-15]10^-9,{x,0,1},{t,0,3.66},PlotTheme->{"Grid","LargeLabels"},AxesLabel->{"x (nm)","t (fs)", "\[Rho]e (nm^-1)"},ImageSize->Medium] Out[10]= Viewing the graph as a movie of probability densities, it can be seen that "center" of the electron moves from side to side of the box: In[11]:= Animate[Plot[\[Rho]e[x,t],{x,0,10^-9},PerformanceGoal->-"Quality",PlotRange->{0,4\*^9},ImageSize->Medium,AxesLabel->{"x (m)", "\[Rho]e (m^-1)"}],{t,0,3.66\*^-15},SaveDefinitions->True] Out[11]= Manipulate[Plot[\[Rho]e[x, t], {x, 0, 10^(-9)}, PerformanceGoal -> -"Quality", PlotRange -> {0, 4000000000}, ImageSize -> Medium, AxesLabel -> {"\!\(\* StyleBox[\"x\", \"SO\"]\) (m)", "\!\(\* StyleBox[\ \"\[Rho]e\", \"SO\"]\) (\!\(\*Supersc ... 15286965\*^48\* (-9.943140748611695\*^-34\*t + 1.8218\*^-39\*x)] + 2\*Cos[1.724444315286965\*^48\*(9.943140748611695\*^-34\*t + 1.8218\*^-39\*x)] + 2\*Cos[(2000000000\*Pi)\*x]))\* Sin[(1000000000\*Pi)\*x]^2}] --- Find the value of a European vanilla call option if the underlying asset price and the strike price are both \$100, the risk-free rate is 5%, the volatility of the underlying asset is 20%, and the maturity period is 1 year, using the Black–Scholes model: ```wl In[1]:= BlackScholesModel = {D[c[t, s], t] + 1 / 2σ ^ 2 s ^ 2 D[c[t, s], {s, 2}] + r s D[c[t, s], s] - r c[t, s] == 0, c[T, s] == Max[s - k, 0]}; ``` Solve the boundary value problem: ```wl In[2]:= dsol = c[t, s] /. DSolve[BlackScholesModel, c[t, s], {t, s}][[1]] Out[2]= (1/2) E^-r T (-E^r t k Erfc[((t - T) (2 r - σ^2) + 2 Log[k] - 2 Log[s]/2 Sqrt[2] Sqrt[-t + T] σ)] + E^r T s Erfc[((t - T) (2 r + σ^2) + 2 Log[k] - 2 Log[s]/2 Sqrt[2] Sqrt[-t + T] σ)]) ``` Compute the value of the European vanilla option: ```wl In[3]:= dsol /. {t -> 0, s -> 100, k -> 100, σ -> 0.2, T -> 1, r -> 0.05} Out[3]= 10.4506 ``` Compare with the value given by ``FinancialDerivative`` : ```wl In[4]:= FinancialDerivative[{"European", "Call"}, {"StrikePrice" -> 100.00, "Expiration" -> 1}, {"InterestRate" -> 0.05, "Volatility" -> 0.2, "CurrentPrice" -> 100}] Out[4]= 10.4506 ``` --- Recover a function from its gradient vector: ```wl In[1]:= sol = DSolve[{D[f[x, y], x] == x y Cos[x y] + Sin[x y], D[f[x, y], y] == -E^-y + x^2 Cos[x y]}, f, {x, y}] Out[1]= {{f -> Function[{x, y}, E^-y + C[1] + x Sin[x y]]}} ``` The solution represents a family of parallel surfaces: ```wl In[2]:= Plot3D[Table[f[x, y, z] /. sol[[1]] /. {C[1] -> k}, {k, 1, 22, 5}]//Evaluate, {x, 0, 10}, {y, 0, 2}, PlotRange -> All] Out[2]= [image] ``` --- Solve a Cauchy problem to generate Stirling numbers: ```wl In[1]:= sol = u[x, y] /. DSolve[{(1 - x)D[u[x, y], x] == y u[x, y], u[0, y] == 1}, u, {x, y}][[1]] Out[1]= (-1)^y (-1 + x)^-y ``` Use the generating function to obtain Stirling numbers: ```wl In[2]:= Table[(-1) ^ i SeriesCoefficient[ Series[sol, {x, 0, 10}, {y, 0, 4}], {10, i}]10!, {i, 4}] Out[2]= {-362880, 1026576, -1172700, 723680} In[3]:= StirlingS1[10, Range[4]] Out[3]= {-362880, 1026576, -1172700, 723680} ``` ### Properties & Relations (10) Solutions satisfy the differential equation and boundary conditions: ```wl In[1]:= DSolve[{y''[x] - y[x] == 0, y[0] == 1, y'[0] == 4}, y, x] Out[1]= {{y -> Function[{x}, (1/2) E^-x (-3 + 5 E^2 x)]}} In[2]:= Simplify[{y''[x] - y[x] == 0, y[0] == 1, y'[0] == 4} /. %] Out[2]= {{True, True, True}} ``` --- Differential equation corresponding to ``Integrate`` : ```wl In[1]:= DSolve[y'[x] == Exp[x ^ 2], y, x] Out[1]= {{y -> Function[{x}, C[1] + (1/2) Sqrt[π] Erfi[x]]}} In[2]:= Integrate[Exp[x ^ 2], x] Out[2]= (1/2) Sqrt[π] Erfi[x] ``` --- Use ``NDSolve`` to find a numerical solution: ```wl In[1]:= exactsol = DSolve[{y''[x] + y[x] == 0, y[0] == 1, y'[0] == 0}, y, x] Out[1]= {{y -> Function[{x}, Cos[x]]}} In[2]:= Table[y[x] /. exactsol[[1]], {x, -2., 2}] Out[2]= {-0.416147, 0.540302, 1., 0.540302, -0.416147} In[3]:= numsol = NDSolve[{y''[x] + y[x] == 0, y[0] == 1, y'[0] == 0}, y, {x, -2, 2}] Out[3]= {{y -> InterpolatingFunction[{{-2., 2.}}, {5, 7, 2, {87}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False}, CompressedData["«720»"], {Developer`PackedArrayForm, CompressedData["«273»"], CompressedData["«1607»"]}, {Automatic}]}} In[4]:= Table[y[x] /. numsol[[1]], {x, -2., 2}] Out[4]= {-0.416147, 0.540302, 1., 0.540302, -0.416147} ``` --- Use ``DEigensystem`` to find eigenvalues and eigenfunctions: ```wl In[1]:= {ℒ, ℬ} = {-Laplacian[u[x], {x}], DirichletCondition[u[x] == 0, True]}; ``` Find the complete eigensystem: ```wl In[2]:= DSolve[{u''[x] + λ u[x] == 0, u[0] == 0, u[π] == 0}, u[x], x] Out[2]= {{u[x] -> Piecewise[{{C[1]*Sin[x*Sqrt[λ]], Element[\[FormalN], Integers] && \[FormalN] >= 1 && λ == \[FormalN]^2}}, 0]}} ``` Find eigenvalues and eigenfunctions: ```wl In[3]:= DEigensystem[{ℒ, ℬ}, u[x], {x, 0, Pi}, 3] Out[3]= {{1, 4, 9}, {Sin[x], Sin[2 x], Sin[3 x]}} ``` --- Compute an impulse response using ``DSolve`` : ```wl In[1]:= DSolve[y'''[x] - 5y''[x] + 9y'[x] - 5y[x] == DiracDelta''[x] + 2DiracDelta'[x] + DiracDelta[x] && y[-1] == 0 && y'[-1] == 0 && y''[-1] == 0, y[x], x]//FullSimplify Out[1]= {{y[x] -> -E^x HeavisideTheta[x] (-2 + E^x (Cos[x] - 7 Sin[x]))}} ``` The same computation using ``InverseLaplaceTransform`` : ```wl In[2]:= InverseLaplaceTransform[(s ^ 2 + 2s + 1/s ^ 3 - 5s ^ 2 + 9s - 5), s, x]//FullSimplify Out[2]= E^x (2 - E^x (Cos[x] - 7 Sin[x])) ``` --- ``DSolve`` returns a rule for the solution: ```wl In[1]:= DSolve[{y'[x] == y[x], y[0] == 1}, y[x], x] Out[1]= {{y[x] -> E^x}} ``` ``DSolveValue`` returns an expression for the solution: ```wl In[2]:= DSolveValue[{y'[x] == y[x], y[0] == 1}, y[x], x] Out[2]= E^x ``` --- Apply ``N[DSolve[...]]`` to invoke ``NDSolve`` if symbolic solution fails: ```wl In[1]:= DSolve[{y''[x] + x ^ 2 y'[x] + Tan[y[x]] == 1, y[0] == 1, y'[0] == 1}, y, {x, 0, 1}] Out[1]= DSolve[{Tan[y[x]] + x^2 Derivative[1][y][x] + Derivative[2][y][x] == 1, y[0] == 1, Derivative[1][y][0] == 1}, y, {x, 0, 1}] In[2]:= N[%] Out[2]= {{y -> InterpolatingFunction[{{0., 1.}}, {5, 7, 2, {61}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False}, CompressedData["«736»"], {Developer`PackedArrayForm, {0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42, 45, 48, 51, 54, 57, 60, 63, 66, 69, 72, 75, 78, 81, 84, 87, 90, 93, 96, 99, 102, 105, 108, 111, 114, 117, 120, 123, 126, 129, 132, 135, 138, 141, 144, 147, 150, 153, 156, 159, 162, 165, 168, 171, 174, 177, 180, 183}, CompressedData["«2061»"]}, {Automatic}]}} In[3]:= Table[y[x] /. %[[1]], {x, 0, 1, 0.2}] Out[3]= {1., 1.18351, 1.30624, 1.3228, 1.23038, 1.08096} ``` --- Solve the Laplace equation in a rectangle: ```wl In[1]:= peqn = Laplacian[u[x, y], {x, y}] == 0; In[2]:= bc = {u[x, 0] == 0, u[x, 1] == 0, u[0, y] == y ^ 2 - y, u[1, y] == 0}; In[3]:= u[x, y] /. DSolve[{peqn, bc}, u[x, y], {x, y}][[1]] Out[3]= Inactive[Sum][(4*(-1 + (-1)^K[1])*Csch[Pi*K[1]]*Sin[Pi*y*K[1]]*Sinh[Pi*(1 - x)*K[1]])/ (Pi^3*K[1]^3), {K[1], 1, Infinity}] ``` Obtain the same result using the region notation: ```wl In[4]:= Ω = Rectangle[{0, 0}, {1, 1}]; In[5]:= rbc = DirichletCondition[u[x, y] == Piecewise[{{y ^ 2 - y, x == 0}}], True]; In[6]:= u[x, y] /. DSolve[{peqn, rbc}, u, {x, y}∈Ω][[1]] Out[6]= Inactive[Sum][(4*(-1 + (-1)^K[1])*Csch[Pi*K[1]]*Sin[Pi*y*K[1]]*Sinh[Pi*(1 - x)*K[1]])/ (Pi^3*K[1]^3), {K[1], 1, Infinity}] ``` --- Use ``DFixedPoints`` to find the fixed points of a system of two ODEs: ```wl In[1]:= DFixedPoints[{x'[t] == -x[t] - y[t] - 1, y'[t] == 2x[t] - y[t] + 5}, {x, y}, t] Out[1]= {{-2, 1}} ``` Use ``DStabilityConditions`` to analyze the stability of the fixed point: ```wl In[2]:= DStabilityConditions[{x'[t] == -x[t] - y[t] - 1, y'[t] == 2x[t] - y[t] + 5}, {x, y}, t] Out[2]= {{{-2, 1}, True}} ``` Solve the system using the fixed point as initial condition: ```wl In[3]:= sol = DSolve[{x'[t] == -x[t] - y[t] - 1, y'[t] == 2x[t] - y[t] + 5, x[0] == -2, y[0] == 1}, {x[t], y[t]}, t]//Simplify Out[3]= {{x[t] -> -2, y[t] -> 1}} ``` Solve the system for given initial conditions: ```wl In[4]:= sol = DSolve[{x'[t] == -x[t] - y[t] - 1, y'[t] == 2x[t] - y[t] + 5, x[0] == 3, y[0] == -2}, {x[t], y[t]}, t]//Simplify Out[4]= {{x[t] -> -2 + 5 E^-t Cos[Sqrt[2] t] + (3 E^-t Sin[Sqrt[2] t]/Sqrt[2]), y[t] -> E^-t (E^t - 3 Cos[Sqrt[2] t] + 5 Sqrt[2] Sin[Sqrt[2] t])}} ``` Plot the solution: ```wl In[5]:= Plot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 20}] Out[5]= [image] ``` --- Use system modeling for numerical solutions to larger hierarchical models: ```wl In[1]:= sim = SystemModelSimulate[[image]] Out[1]= SystemModelSimulationData["/tmp/WolframSystemModeler-secavarat-14.3/WL/wsml_ControlledMotorGearPend\ ulum_lsim_103635_45794_94847052752793608981826305086_1.mat", "DocumentationExamples.Modeling.ControlledMotorGearPendulum", 741, {0., 20.}, "/tmp/ ... ^-6, "DisplayTimeUnit" -> "s", "Epoch" -> None, "DisplayTimeZone" -> Automatic, "InterpolationPoints" -> 2000, "Method" -> Automatic, "StepSize" -> 0, "SynchronizeWithRealTime" -> False, "SimulationRate" -> 1., "CommunicationRate" -> 20]] ``` Plot the most important values: ```wl In[2]:= SystemModelPlot[sim] Out[2]= [image] ``` ### Possible Issues (9) Results may contain inactive integrals: ```wl In[1]:= sol = DSolve[y'[x] == f[x], y[x], x] Out[1]= {{y[x] -> C[1] + Inactive[Integrate][f[K[1]], {K[1], 1, x}]}} ``` Replacing the function $f$ by an expression still returns inactive integrals: ```wl In[2]:= sol /. f -> Function[x, 0] Out[2]= {{y[x] -> C[1] + Inactive[Integrate][0, {K[1], 1, x}]}} ``` Use ``Activate`` to evaluate integrals: ```wl In[3]:= Activate[%] Out[3]= {{y[x] -> C[1]}} ``` --- Capital $K$ and capital $C$ cannot be used as independent variables: ```wl In[1]:= DSolve[x'[K] == K ^ K, x[K], K] ``` DSolve::indep: K cannot be used as an independent variable. ```wl Out[1]= DSolve[Derivative[1][x][K] == K^K, x[K], K] In[2]:= DSolve[x'[C] == C ^ C, x[C], C] ``` DSolve::indep: C cannot be used as an independent variable. ```wl Out[2]= DSolve[Derivative[1][x][C] == C^C, x[C], C] ``` Replacing them by lowercase $k$ or lowercase $c$ fixes the issue: ```wl In[3]:= DSolve[x'[k] == k ^ k, x[k], k] Out[3]= {{x[k] -> C[1] + Inactive[Integrate][K[1]^K[1], {K[1], 1, k}]}} In[4]:= DSolve[x'[c] == c ^ c, x[c], c] Out[4]= {{x[c] -> C[1] + Inactive[Integrate][K[1]^K[1], {K[1], 1, c}]}} ``` --- Inverse functions may be required to find the solution: ```wl In[1]:= DSolve[{y'[x] ^ 2 == (1 - y[x] ^ 2)(1 - (1 / 2)y[x] ^ 2), y[0] == 0}, y, x] ``` DSolve::ifun: Inverse functions are being used by DSolve, so some solutions may not be found. ```wl Out[1]= {{y -> Function[{x}, -JacobiSN[x, (1/2)]]}, {y -> Function[{x}, JacobiSN[x, (1/2)]]}} ``` --- Definitions for an unknown function may affect the evaluation: ```wl In[1]:= y[x_] := Cos[x] In[2]:= DSolve[{y'[x] == 0, y[0] == 1}, y, x] ``` DSolve::deqn: Equation or list of equations expected instead of True in the first argument {-Sin[x]==0,True}. ```wl Out[2]= DSolve[{-Sin[x] == 0, True}, y, x] ``` Clearing the definition for the unknown function fixes the issue: ```wl In[3]:= Clear[y] In[4]:= DSolve[{y'[x] == 0, y[0] == 1}, y, x] Out[4]= {{y -> Function[{x}, 1]}} ``` --- Here ``DSolve`` gives an error message because $z(x)$ is not included in the second argument: ```wl In[1]:= DSolve[{y'[x] == 3y[x], z'[x] == 2 z[x], z[0] == 1, y[0] == 1}, y[x], x] ``` DSolve::deqx: Supplied equations are not differential or integral equations of the given functions. ```wl Out[1]= DSolve[{Derivative[1][y][x] == 3 y[x], Derivative[1][z][x] == 2 z[x], z[0] == 1, y[0] == 1}, y[x], x] ``` Include $z(x)$ in the second argument to obtain a solution: ```wl In[2]:= DSolve[{y'[x] == 3y[x], z'[x] == 2 z[x], z[0] == 1, y[0] == 1}, {y[x], z[x]}, x] Out[2]= {{y[x] -> E^3 x, z[x] -> E^2 x}} ``` Alternatively, use ``DSolveValue`` to obtain the solution for $y(x)$ alone: ```wl In[3]:= DSolveValue[{y'[x] == 3y[x], z'[x] == 2 z[x], z[0] == 1, y[0] == 1}, y[x], x] Out[3]= E^3 x ``` --- This example returns unevaluated because ``DSolve`` is unable to interpret it as a vector equation: ```wl In[1]:= A = {{2, 0}, {-2, 0}}; In[2]:= B = {3t, 4}; In[3]:= eqn = x'[t] == A.x[t] + B Out[3]= Derivative[1][x][t] == {3 t + {{2, 0}, {-2, 0}}.x[t], 4 + {{2, 0}, {-2, 0}}.x[t]} In[4]:= DSolve[eqn, x∈Vectors[2], t] Out[4]= DSolve[Derivative[1][x][t] == {3 t + {{2, 0}, {-2, 0}}.x[t], 4 + {{2, 0}, {-2, 0}}.x[t]}, x∈Vectors[2, ℂ], t] ``` The example can be solved by rewriting it as follows: ```wl In[5]:= eqn1 = x'[t] + A.x[t] == B Out[5]= {{2, 0}, {-2, 0}}.x[t] + Derivative[1][x][t] == {3 t, 4} In[6]:= DSolve[eqn1, x∈Vectors[2], t] Out[6]= {{x -> Function[{t}, {3 (-(1/4) + (t/2)) + E^-2 t C[1], -3 E^2 t (-(1/4) + (t/2)) + 3 (-1 + E^2 t) (-(1/4) + (t/2)) + 4 t + (3 t^2/2) + E^-2 t (-1 + E^2 t) C[1] + C[2]}]}} In[7]:= eqn1 /. %[[1]]//Simplify Out[7]= True ``` --- ``CompleteIntegral`` finds a complete integral for a nonlinear PDE: ```wl In[1]:= deqn = D[u[x, y], x] + D[u[x, y], y] ^ 2 == 3; In[2]:= CompleteIntegral[deqn, u[x, y], {x, y}] Out[2]= {{u[x, y] -> C[1] + y C[2] + x (3 - C[2]^2)}} ``` ``DSolve`` returns the same solution with a warning message: ```wl In[3]:= DSolve[deqn, u[x, y], {x, y}] ``` DSolve::nlpde: Solution requested to nonlinear partial differential equation. Trying to build a special solution. ```wl Out[3]= {{u[x, y] -> C[1] + y C[2] + x (3 - C[2]^2)}} ``` --- Use ``CompleteIntegral`` to find a complete integral for a linear PDE: ```wl In[1]:= deqn = D[u[x, y], x] + 2D[u[x, y], y] == 1; In[2]:= CompleteIntegral[deqn, u[x, y], {x, y}] Out[2]= {{u[x, y] -> C[1] + x (1 - 2 C[2]) + y C[2]}} ``` ``DSolve`` returns the general solution for this PDE: ```wl In[3]:= DSolve[deqn, u[x, y], {x, y}] Out[3]= {{u[x, y] -> x + C[1][-2 x + y]}} ``` --- Solve a second-order linear ODE with symbolic parameters: ```wl In[1]:= eqn = y''[x] + (k + g1 / (a + x)) y'[x] + (m + g2 / (a + x)) y[x] == 0; In[2]:= sol = DSolve[eqn, y[x], x] Out[2]= {{y[x] -> E^-(k x/2) - (1/2) Sqrt[k^2 - 4 m] x C[1] HypergeometricU[-(2 g2 - g1 k - g1 Sqrt[k^2 - 4 m]/2 Sqrt[k^2 - 4 m]), g1, a Sqrt[k^2 - 4 m] + Sqrt[k^2 - 4 m] x] + E^-(k x/2) - (1/2) Sqrt[k^2 - 4 m] x C[2] LaguerreL[(2 g2 - g1 k - g1 Sqrt[k^2 - 4 m]/2 Sqrt[k^2 - 4 m]), -1 + g1, a Sqrt[k^2 - 4 m] + Sqrt[k^2 - 4 m] x]}} ``` The solution obtained is not the general solution of the ODE for specific values of parameters $g_1=g_2=0$ : ```wl In[3]:= y[x] /. sol[[1]] /. {g1 -> 0, g2 -> 0}//Simplify Out[3]= E^-(1/2) (k + Sqrt[k^2 - 4 m]) x (C[1] + C[2]) ``` The general solution for this case can be found by substituting parameter values in the equation first: ```wl In[4]:= DSolve[eqn /. {g1 -> 0, g2 -> 0}, y[x], x] Out[4]= {{y[x] -> E^(1/2) (-k - Sqrt[k^2 - 4 m]) x C[1] + E^(1/2) (-k + Sqrt[k^2 - 4 m]) x C[2]}} ``` ### Neat Examples (2) Generate a Cornu spiral: ```wl In[1]:= DSolve[{x'[s] == Cos[t[s]], y'[s] == Sin[t[s]], t'[s] == s, x[0] == 0, y[0] == 0, t[0] == 0}, {x, y, t}, s] Out[1]= {{t -> Function[{s}, (s^2/2)], x -> Function[{s}, Sqrt[π] FresnelC[(s/Sqrt[π])]], y -> Function[{s}, Sqrt[π] FresnelS[(s/Sqrt[π])]]}} In[2]:= ParametricPlot[Evaluate[{x[s], y[s]} /. %], {s, -10, 10}] Out[2]= [image] ``` --- Solve the sixth symmetric power of the Legendre differential operator: ```wl In[1]:= sol = DSolve[(-((192n(1 + n)x ^ 5) / (1 - x ^ 2) ^ 6) + (2496n(1 + n)x ^ 3) / (1 - x ^ 2) ^ 5 - (8160n ^ 2(1 + n) ^ 2x ^ 3) / (1 - x ^ 2) ^ 5 - (1632n(1 + n)x) / (1 - x ^ 2) ^ 4 + (7824n ^ 2(1 + n) ^ 2x) / (1 - x ^ 2) ^ 4 - (6912n ^ 3(1 + n) ^ 3x) / (1 - x ^ 2) ^ 4)y[x] + ((64x ^ 6) / (1 - x ^ 2) ^ 6 - (1824x ^ 4) / (1 - x ^ 2) ^ 5 + (8256n(1 + n)x ^ 4) / (1 - x ^ 2) ^ 5 + (2880x ^ 2) / (1 - x ^ 2) ^ 4 - (19200n(1 + n)x ^ 2) / (1 - x ^ 2) ^ 4 + (22896n ^ 2(1 + n) ^ 2x ^ 2) / (1 - x ^ 2) ^ 4 - 272 / (1 - x ^ 2) ^ 3 + (2208n(1 + n)) / (1 - x ^ 2) ^ 3 - (4384n ^ 2(1 + n) ^ 2) / (1 - x ^ 2) ^ 3 + (2304n ^ 3(1 + n) ^ 3) / (1 - x ^ 2) ^ 3) y'[x] + (-((2016x ^ 5) / (1 - x ^ 2) ^ 5) + (9408x ^ 3) / (1 - x ^ 2) ^ 4 - (19488n(1 + n)x ^ 3) / (1 - x ^ 2) ^ 4 - (3696x) / (1 - x ^ 2) ^ 3 + (12768n(1 + n)x) / (1 - x ^ 2) ^ 3 - (9408n ^ 2(1 + n) ^ 2x) / (1 - x ^ 2) ^ 3)y''[x] + ((4816x ^ 4) / (1 - x ^ 2) ^ 4 - (7168x ^ 2) / (1 - x ^ 2) ^ 3 + (9632n(1 + n)x ^ 2) / (1 - x ^ 2) ^ 3 + 616 / (1 - x ^ 2) ^ 2 - (1456n(1 + n)) / (1 - x ^ 2) ^ 2 + (784n ^ 2(1 + n) ^ 2) / (1 - x ^ 2) ^ 2)y'''[x] + (-((2800x ^ 3) / (1 - x ^ 2) ^ 3) + (1400x) / (1 - x ^ 2) ^ 2 - (1400n(1 + n)x) / (1 - x ^ 2) ^ 2) y''''[x] + ((560x ^ 2) / (1 - x ^ 2) ^ 2 - 70 / (1 - x ^ 2) + (56n(1 + n)) / (1 - x ^ 2))y'''''[x] - (42x y''''''[x]) / (1 - x ^ 2) + y'''''''[x] == 0, y, x] Out[1]= {{y -> Function[{x}, C[1] LegendreP[n, x]^6 + C[2] LegendreP[n, x]^5 LegendreQ[n, x] + C[3] LegendreP[n, x]^4 LegendreQ[n, x]^2 + C[4] LegendreP[n, x]^3 LegendreQ[n, x]^3 + C[5] LegendreP[n, x]^2 LegendreQ[n, x]^4 + C[6] LegendreP[n, x] LegendreQ[n, x]^5 + C[7] LegendreQ[n, x]^6]}} ``` ## See Also * [`DSolveValue`](https://reference.wolfram.com/language/ref/DSolveValue.en.md) * [`NDSolve`](https://reference.wolfram.com/language/ref/NDSolve.en.md) * [`AsymptoticDSolveValue`](https://reference.wolfram.com/language/ref/AsymptoticDSolveValue.en.md) * [`WhenEvent`](https://reference.wolfram.com/language/ref/WhenEvent.en.md) * [`DEigensystem`](https://reference.wolfram.com/language/ref/DEigensystem.en.md) * [`DEigenvalues`](https://reference.wolfram.com/language/ref/DEigenvalues.en.md) * [`DFixedPoints`](https://reference.wolfram.com/language/ref/DFixedPoints.en.md) * [`DStabilityConditions`](https://reference.wolfram.com/language/ref/DStabilityConditions.en.md) * [`NDEigensystem`](https://reference.wolfram.com/language/ref/NDEigensystem.en.md) * [`NDEigenvalues`](https://reference.wolfram.com/language/ref/NDEigenvalues.en.md) * [`GreenFunction`](https://reference.wolfram.com/language/ref/GreenFunction.en.md) * [`CompleteIntegral`](https://reference.wolfram.com/language/ref/CompleteIntegral.en.md) * [`Solve`](https://reference.wolfram.com/language/ref/Solve.en.md) * [`RSolve`](https://reference.wolfram.com/language/ref/RSolve.en.md) * [`Integrate`](https://reference.wolfram.com/language/ref/Integrate.en.md) * [`DifferentialRoot`](https://reference.wolfram.com/language/ref/DifferentialRoot.en.md) * [`StreamPlot`](https://reference.wolfram.com/language/ref/StreamPlot.en.md) * [`ItoProcess`](https://reference.wolfram.com/language/ref/ItoProcess.en.md) * [`SystemModelSimulate`](https://reference.wolfram.com/language/ref/SystemModelSimulate.en.md) * [`FractionalD`](https://reference.wolfram.com/language/ref/FractionalD.en.md) * [`CaputoD`](https://reference.wolfram.com/language/ref/CaputoD.en.md) * [`TruncateSum`](https://reference.wolfram.com/language/ref/TruncateSum.en.md) ## Tech Notes * [Differential Equations](https://reference.wolfram.com/language/tutorial/Calculus.en.md#9252) * [Differential Equation Solving with DSolve](https://reference.wolfram.com/language/tutorial/DSolveOverview.en.md) * [Symbolic Solutions of PDEs](https://reference.wolfram.com/language/tutorial/SymbolicSolutionsOfPDEs.en.md) * [Introduction to Fractional Calculus](https://reference.wolfram.com/language/tutorial/FractionalCalculus.en.md) * [Numerical Differential Equation Solving with NDSolve](https://reference.wolfram.com/language/tutorial/NDSolveOverview.en.md) * [Implementation notes: Algebra and Calculus](https://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.en.md#8313) ## Related Guides * [Equation Solving](https://reference.wolfram.com/language/guide/EquationSolving.en.md) * [`Calculus`](https://reference.wolfram.com/language/guide/Calculus.en.md) * [Differential Equations](https://reference.wolfram.com/language/guide/DifferentialEquations.en.md) * [Differential Operators](https://reference.wolfram.com/language/guide/DifferentialOperators.en.md) * [Partial Differential Equations](https://reference.wolfram.com/language/guide/PDEModelingAndAnalysis.en.md) * [Mathematical Data](https://reference.wolfram.com/language/guide/MathematicalData.en.md) * [Fractional Calculus](https://reference.wolfram.com/language/guide/FractionalCalculus.en.md) * [Scientific Models](https://reference.wolfram.com/language/guide/ScientificModels.en.md) * [Integral Transforms](https://reference.wolfram.com/language/guide/IntegralTransforms.en.md) * [Symbolic Vectors, Matrices and Arrays](https://reference.wolfram.com/language/guide/SymbolicArrays.en.md) * [Solvers over Regions](https://reference.wolfram.com/language/guide/GeometricSolvers.en.md) ## History * Introduced in 1991 (2.0) \| Updated in 1996 (3.0) ▪ 2003 (5.0) ▪ [2014 (10.0)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn100.en.md) ▪ [2015 (10.3)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn103.en.md) ▪ [2016 (11.0)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn110.en.md) ▪ [2019 (12.0)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn120.en.md) ▪ [2021 (13.0)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn130.en.md) ▪ [2022 (13.1)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn131.en.md) ▪ [2022 (13.2)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn132.en.md) ▪ [2024 (14.1)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn141.en.md) ▪ [2025 (14.2)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn142.en.md)