Dsolve Unable to Find Explicit Solution

Differential Equations

Robert G. Mortimer , in Mathematics for Physical Chemistry (Fourth Edition), 2013

Symbolic Solution

The statement DSolve is used to carry out a symbolic solution of a differential equation. We illustrate this with an example.

Example 12.14

Use Mathematica to solve the differential equation

$\frac{\mathrm{d}y}{\mathrm{d}x}=\mathit{ay}(x)\text{.}$

We enter the Mathematica statement

DSolve[y'[x]= =a y[x], y[x], x]

and press the "Enter" key. Notice how the statement is written inside the brackets. First comes the equation, with the first derivative denoted by y'. The double equal sign must be used to let Mathematica know that an equation is to be solved. We have used a blank space between the a and the y[x] to indicate multiplication. After a comma comes the specification of the dependent variable, y[x]. Note the use of brackets, not parentheses. The independent variable must be included inside the brackets. After another comma comes the statement of the independent variable. Mathematica returns the output

Out[1]={{y[x] $\to$ e^{a   x} C[1]}}

Note the space between the a and the x and the space between the exponential and the constant C[1] in the output. The constant C[1] is to be determined by initial conditions. An initial condition can be included in the original input statement. For example, if y(0)   =   2, we would enter

DSolve[{y'[x]==a y[x],y[0]==2},y[x],x]

and press the "Enter" key or a "Shift-Return." The output would be

Out[1]={{y[x] $\to$ 2 e^{a x}}}

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780124158092000124

Higher-Order Differential Equations

Martha L. Abell , James P. Braselton , in Differential Equations with Mathematica (Fourth Edition), 2016

4.6.3 Variation of Parameters

Cauchy-Euler equations can be homogenous or nonhomogeneous in which case the method of Variation of Parameters can be used to solve the problem. Before implementing the method to find a particular solution, be sure to write the equation in standard form first.

Example 4.6.4

Solve x ² y″ − xy′ + 5y = x, x > 0.

Solution

We first note that DSolve can be used to find a general solution of the equation directly.

Clear[x, y, gensol]gensol=DSolve[x ^∧ 2y"[x]−xy′[x]+5y[x]==x, y[x], x]

$\{\{y[x]\to xC[2]\text{Cos}[2\text{Log}[x]]+xC[1]\text{Sin}[2\text{Log}[x]]+\frac{1}{4}(2x\text{Cos}{[\text{Log}[x]]}^2$ Cos[2Log[x]] + xSin[2Log[x]]²)}}

Alternatively, we can use Mathematica to help us implement Variation of Parameters. We begin by finding a general solution to the corresponding homogeneous equation x ² y″ − xy′ + 5y = 0 with DSolve.

homsol=DSolve[x ^∧ 2y"[x]−xy′[x]+5y[x]==0, y[x], x]

{{y[x] → xC[2]Cos[2Log[x]] + xC[1]Sin[2Log[x]]}}

We see that a general solution of the corresponding homogeneous equation is $y_h=x\left[c_{1}cos\left(2lnx\right)+c_{2}sin\left(2lnx\right)\right]$ . A fundamental set of solutions for the homogeneous equation is $S=\left\{xcos\left(2lnx\right),xsin\left(2lnx\right)\right\}$

y1[x_]=xCos[2Log[x]];y2[x_]=xSin[2Log[x]];caps={y1[x], y2[x]};

and the Wronskian is W(S) = 2x.

ws=Wronskian[caps, x]

2x

To implement Variation of Parameters, we rewrite the equation in standard form

$\begin{array}{l}{y}^{″}-\frac{1}{x}{y}^{\prime }+\frac{5}{x^2}y=\frac{1}{x}\hfill \end{array}$

by dividing by x ² and identify f(x) = 1/x. We then use Integrate to compute

$\begin{array}{l}{u}_1=\int \frac{-y_2(x)f(x)}{2x}dx\text{and}{u}_2=\int \frac{y_1(x)f(x)}{2x}dx.\hfill \end{array}$

f[x_]=1/x;

u1p=−y2[x]f[x]/ws u2p=y1[x]f[x]/ws

$-\frac{\text{Sin}[2\text{Log}[x]]}{2x}$

$\frac{\text{Cos}[2\text{Log}[x]]}{2x}$

u1[x_]=Integrate[u1p, x]u2[x_]=Integrate[u2p, x]

$\frac{1}{2}\text{Cos}{[\text{Log}[x]]}^2$

$\frac{1}{4}\text{Sin}[2\text{Log}[x]]$

A particular solution of the nonhomogeneous equation is given by y _p = y ₁ u ₁ + y ₂ u ₂

yp[x_]=u1[x]y1[x]+u2[x]y2[x]//Simplify

$\frac{1}{2}x\text{Cos}{[\text{Log}[x]]}^2$

and a general solution is given by y = y _h + y _p .

y[x_]=homsol[[1, 1, 2]]+yp[x]

$\frac{1}{2}x\text{Cos}{[\text{Log}[x]]}^2+xC[2]\text{Cos}[2\text{Log}[x]]+xC[1]\text{Sin}[2\text{Log}[x]]$

As in previous examples, we graph this general solution for various values of the arbitrary constants. See Figure 4-30.

Figure 4-30. Various solutions of a nonhomogeneous Cauchy-Euler equation

toplot=Table[y[x]/.{C[1]→i, C[2]→j}, {i, −2, 2}, {j, −2, 2}];

Plot[toplot, {x, 0, 2}, AxesLabel→{x, y}]

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128047767000048

Applications related to ordinary and partial differential equations

Martha L. Abell , James P. Braselton , in Mathematica by Example (Sixth Edition), 2022

6.1.4 Numerical methods

If numerical results are desired, use NDSolve:

NDSolve[{y'[t]==f[t,y[t]],y[t0]==y0},y[t],{t,a,b}]

attempts to generate a numerical solution of $dy/dt=f(t,y)$ , $y\left(t_0\right)=y_0$ , valid for $a⩽t⩽b$ .

Example 6.7

Consider $dy/dt=(t^2-y^2)\sin y$ , $(0)=-1$ . (a) Determine $y(1)$ . (b) Graph $y(t)$ , $-1⩽t⩽10$ .

Solution

We first remark that DSolve can neither exactly solve the differential equation $y^{\prime }=(t^2-y^2)\sin y$ nor find the solution that satisfies $y(0)=-1$ .

When Mathematica returns a command unevaluated, interpret the result to mean that Mathematica cannot evaluate the command.

$\mathbf{\text{sol}}\mathbf{=}\mathbf{\text{DSolve}}\mathbf{[}{\mathit{y}}^{\mathbf{\prime }}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{=}\mathbf{=}\mathbf{(}{\mathit{t}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{-}\mathit{y}{\mathbf{[}\mathit{t}\mathbf{]}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{)}\mathbf{\text{Sin}}\mathbf{[}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{]}\mathbf{,}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{,}\mathit{t}\mathbf{]}$

$\text{DSolve}[y^{\prime }[t]==\text{Sin}[y[t]](t^2-y{[t]}^2),y[t],t]$

$\mathbf{\text{sol}}\mathbf{=}\mathbf{\text{DSolve}}\mathbf{[}\mathbf{\{}{\mathit{y}}^{\mathbf{\prime }}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{=}\mathbf{=}\mathbf{(}{\mathit{t}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{-}\mathit{y}{\mathbf{[}\mathit{t}\mathbf{]}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{)}\mathbf{\text{Sin}}\mathbf{[}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{]}\mathbf{,}\mathit{y}\mathbf{[}\mathbf{0}\mathbf{]}\mathbf{=}\mathbf{=}\mathbf{-}\mathbf{1}\mathbf{\}}\mathbf{,}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{,}\mathit{t}\mathbf{]}$

$\text{DSolve}[\{y^{\prime }[t]==\text{Sin}[y[t]](t^2-y{[t]}^2),y[0]==-1\},y[t],t]$

However, we obtain a numerical solution valid for $-1⩽t⩽10$ using the NDSolve function.

$\mathbf{\text{sol}}\mathbf{=}\mathbf{\text{NDSolve}}\mathbf{[}\mathbf{\{}{\mathit{y}}^{\mathbf{\prime }}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{=}\mathbf{=}\mathbf{(}{\mathit{t}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{-}\mathit{y}{\mathbf{[}\mathit{t}\mathbf{]}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{)}\mathbf{\text{Sin}}\mathbf{[}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{]}\mathbf{,}\mathit{y}\mathbf{[}\mathbf{0}\mathbf{]}\mathbf{=}\mathbf{=}\mathbf{-}\mathbf{1}\mathbf{\}}\mathbf{,}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{,}$

$\mathbf{\{}\mathit{t}\mathbf{,}\mathbf{-}\mathbf{1}\mathbf{,}\mathbf{10}\mathbf{\}}\mathbf{]}$

$\{\{y[t]\to \text{InterpolatingFunction}[][t]\}\}$

Entering sol /.t->1 evaluates the numerical solution if $t=1$ .

$\mathbf{\text{sol}}\mathbf{\text{/.}}\mathit{t}\mathbf{\to }\mathbf{1}$

$\{\{y[1]\to -0.766013\}\}$

The result means that $y(1)\approx -.766$ . We use the Plot command to graph the solution for $0⩽t⩽10$ in Fig. 6.9.

Figure 6.9. Graph of the solution to $y^{\prime }=(t^2-y^2)\sin y$ , y(0)=−1 (University of Connecticut colors).

$\mathbf{\text{Plot}}\mathbf{[}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{\text{/.}}\mathbf{\text{sol}}\mathbf{,}\mathbf{\{}\mathit{t}\mathbf{,}\mathbf{-}\mathbf{1}\mathbf{,}\mathbf{10}\mathbf{\}}\mathbf{,}$

$\mathbf{\text{PlotStyle}}\mathbf{\to }\mathbf{\text{CMYKColor}}\mathbf{[}\mathbf{1}\mathbf{,}\mathbf{.}\mathbf{76}\mathbf{,}\mathbf{.}\mathbf{12}\mathbf{,}\mathbf{.}\mathbf{7}\mathbf{]}\mathbf{]}$   □

Example 6.8 Logistic equation with predation

Incorporating predation into the logistic equation, $y^{\prime }=\alpha y(1-\frac{1}{K}y)$ , results in $\frac{dy}{dt}=\alpha y(1-\frac{1}{K}y)-P(y)$ , where $P(y)$ is a function of y describing the rate of predation. A typical choice for P is $P(y)=a{y}^2/(b^2+y^2)$ , because $P(0)=0$ and P is bounded above: ${\lim }_{t\to \infty }P(y)<\infty$ . Of course, if ${\lim }_{t\to \infty }y(t)=Y$ , then ${\lim }_{t\to \infty }P(y)=a{Y}^2/(b^2+Y^2)$ . Generally, however, ${\lim }_{t\to \infty }P(y)\ne a$ , because ${\lim }_{t\to \infty }y(t)⩽K\ne \infty$ , for some $K⩾0$ , in the predation situation.

If $\alpha =1$ , $a=5$ , and $b=2$ , graph the direction field associated with the equation as well as various solutions if (a) $K=19$ and (b) $K=20$ .

Solution

(a) We define eqn[k] to be $\frac{dy}{dt}=y(1-\frac{1}{K}y)-\frac{5y^2}{4+y^2}$ .

$\mathbf{\text{eqn}}\mathbf{[}\mathbf{\text{k}}\mathrm{\_}\mathbf{]}\mathbf{=}{\mathit{y}}^{\mathbf{\prime }}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{=}\mathbf{=}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{(}\mathbf{1}\mathbf{-}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{/}\mathit{k}\mathbf{)}\mathbf{-}\mathbf{5}\mathit{y}{\mathbf{[}\mathit{t}\mathbf{]}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{/}\mathbf{(}\mathbf{4}\mathbf{+}\mathit{y}{\mathbf{[}\mathit{t}\mathbf{]}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{)}\mathbf{;}$

We use StreamPlot to graph the direction field in Fig. 6.10 (a), and then the direction field along with the solutions that satisfy $y(0)=.5$ , $y(0)=.2$ , and $y(0)=4$ in Fig. 6.10 (b).

Figure 6.10. (a) Direction field. (b) Direction field with three solutions (University of Idaho colors).

$\mathbf{\text{pvf19}}\mathbf{=}\mathbf{\text{StreamPlot}}\mathbf{[}\mathbf{\{}\mathbf{1}\mathbf{,}\mathit{y}\mathbf{(}\mathbf{1}\mathbf{-}\mathbf{1}\mathbf{/}\mathbf{19}\mathit{y}\mathbf{)}\mathbf{-}\mathbf{5}{\mathit{y}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{/}\mathbf{(}\mathbf{4}\mathbf{+}{\mathit{y}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{)}\mathbf{\}}\mathbf{,}\mathbf{\{}\mathit{t}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{10}\mathbf{\}}\mathbf{,}$

$\mathbf{\{}\mathit{y}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{6}\mathbf{\}}\mathbf{,}\mathbf{\text{StreamStyle}}\mathbf{\to }\mathbf{\{}\mathbf{\text{Fine}}\mathbf{,}\mathbf{\text{CMYKColor}}\mathbf{[}\mathbf{0}\mathbf{,}\mathbf{.}\mathbf{23}\mathbf{,}\mathbf{.}\mathbf{53}\mathbf{,}\mathbf{.}\mathbf{35}\mathbf{]}\mathbf{\}}\mathbf{,}$

$\mathbf{\text{Axes}}\mathbf{\to }\mathbf{\text{Automatic}}\mathbf{,}\mathbf{\text{AxesOrigin}}\mathbf{\to }\mathbf{\{}\mathbf{0}\mathbf{,}\mathbf{0}\mathbf{\}}\mathbf{,}$

$\mathbf{\text{Frame}}\mathbf{\to }\mathbf{\text{False}}\mathbf{]}$

$\mathbf{\text{numsols}}\mathbf{=}\mathbf{\text{Map}}\mathbf{[}\mathbf{\text{NDSolve}}\mathbf{[}\mathbf{\{}\mathbf{\text{eqn}}\mathbf{[}\mathbf{19}\mathbf{]}\mathbf{,}\mathit{y}\mathbf{[}\mathbf{0}\mathbf{]}\mathbf{=}\mathbf{=}\mathbf{\#}\mathbf{\}}\mathbf{,}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{,}\mathbf{\{}\mathit{t}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{10}\mathbf{\}}\mathbf{]}\mathbf{\&}\mathbf{,}$

$\mathbf{\text{Table}}\mathbf{[}\mathit{i}\mathbf{,}\mathbf{\{}\mathit{i}\mathbf{,}\mathbf{0.5}\mathbf{,}\mathbf{6}\mathbf{,}\mathbf{5.5}\mathbf{/}\mathbf{9}\mathbf{\}}\mathbf{]}\mathbf{]}\mathbf{;}$

$\mathbf{\text{solplot}}\mathbf{=}\mathbf{\text{Plot}}\mathbf{[}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{\text{/.}}\mathbf{\text{numsols}}\mathbf{,}\mathbf{\{}\mathit{t}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{10}\mathbf{\}}\mathbf{,}\mathbf{\text{PlotRange}}\mathbf{\to }\mathbf{\text{All}}\mathbf{,}$

$\mathbf{\text{PlotStyle}}\mathbf{\to }\mathbf{\{}\mathbf{\{}\mathbf{\text{CMYKColor}}\mathbf{[}\mathbf{0}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{.}\mathbf{5}\mathbf{]}\mathbf{,}\mathbf{\text{Thickness}}\mathbf{[}\mathbf{.}\mathbf{01}\mathbf{]}\mathbf{\}}\mathbf{\}}\mathbf{]}$

$\mathbf{\text{Show}}\mathbf{[}\mathbf{\text{GraphicsRow}}\mathbf{[}\mathbf{\{}\mathbf{\text{pvf19}}\mathbf{,}\mathbf{\text{Show}}\mathbf{[}\mathbf{\text{pvf19}}\mathbf{,}\mathbf{\text{solplot}}\mathbf{]}\mathbf{\}}\mathbf{]}\mathbf{]}$

In the plot, notice that all nontrivial solutions appear to approach an equilibrium solution. We determine the equilibrium solution by solving $y^{\prime }=0$ ,

$\mathbf{\text{eqn}}\mathbf{[}\mathbf{19}\mathbf{]}\mathbf{[}\mathbf{[}\mathbf{2}\mathbf{]}\mathbf{]}$

$(1-\frac{y[t]}{19})y[t]-\frac{5y{[t]}^2}{4+y{[t]}^2}$

$\mathbf{\text{Solve}}\mathbf{[}\mathbf{\text{eqn}}\mathbf{[}\mathbf{19}\mathbf{.}\mathbf{]}\mathbf{[}\mathbf{[}\mathbf{2}\mathbf{]}\mathbf{]}\mathbf{\text{==}}\mathbf{0}\mathbf{,}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{]}$

$\{\{y[t]\to 0.\},\{y[t]\to 0.923351\},\{y[t]\to 9.03832-0.785875i\},$

$\{y[t]\to 9.03832+0.785875i\}\}$ ,

to see that the equilibrium solution ( $dy/dt=0$ ) is $y\approx 0.923$ .

(b) We carry out similar steps for (b). First, we graph the direction field with StreamPlot in Fig. 6.11 (a).

Figure 6.11. (a) Direction field. (b) Direction field with several solutions.

$\mathbf{\text{pvf20}}\mathbf{=}\mathbf{\text{StreamPlot}}\mathbf{[}\mathbf{\{}\mathbf{1}\mathbf{,}\mathit{y}\mathbf{(}\mathbf{1}\mathbf{-}\mathbf{1}\mathbf{/}\mathbf{20}\mathit{y}\mathbf{)}\mathbf{-}\mathbf{5}{\mathit{y}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{/}\mathbf{(}\mathbf{4}\mathbf{+}{\mathit{y}}^{\mathbf{\wedge }}\mathbf{2}\mathbf{)}\mathbf{\}}\mathbf{,}\mathbf{\{}\mathit{t}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{10}\mathbf{\}}\mathbf{,}$

$\mathbf{\{}\mathit{y}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{20}\mathbf{\}}\mathbf{,}\mathbf{\text{Axes}}\mathbf{\to }\mathbf{\text{Automatic}}\mathbf{,}\mathbf{\text{AxesOrigin}}\mathbf{\to }\mathbf{\{}\mathbf{0}\mathbf{,}\mathbf{0}\mathbf{\}}\mathbf{,}$

$\mathbf{\text{AspectRatio}}\mathbf{\to }\mathbf{1}\mathbf{/}\mathbf{\text{GoldenRatio}}\mathbf{,}$

$\mathbf{\text{StreamStyle}}\mathbf{\to }\mathbf{\{}\mathbf{\text{Fine}}\mathbf{,}\mathbf{\text{CMYKColor}}\mathbf{[}\mathbf{0}\mathbf{,}\mathbf{.}\mathbf{23}\mathbf{,}\mathbf{.}\mathbf{53}\mathbf{,}\mathbf{.}\mathbf{35}\mathbf{]}\mathbf{\}}\mathbf{,}$

$\mathbf{\text{Frame}}\mathbf{\to }\mathbf{\text{False}}\mathbf{]}\mathbf{;}$

We then use Map together with NDSolve to numerically find the solution satisfying $y(0)=i$ , for 20 equally spaced values of i between 0 and 20 and name the resulting list numsols. The functions contained in numsols are graphed with Plot in solplot.

$\mathbf{\text{numsols}}\mathbf{=}\mathbf{\text{Map}}\mathbf{[}\mathbf{\text{NDSolve}}\mathbf{[}\mathbf{\{}\mathbf{\text{eqn}}\mathbf{[}\mathbf{20}\mathbf{]}\mathbf{,}\mathit{y}\mathbf{[}\mathbf{0}\mathbf{]}\mathbf{=}\mathbf{=}\mathbf{\#}\mathbf{\}}\mathbf{,}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{,}\mathbf{\{}\mathit{t}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{10}\mathbf{\}}\mathbf{]}\mathbf{\&}\mathbf{,}$

$\mathbf{\text{Table}}\mathbf{[}\mathit{i}\mathbf{,}\mathbf{\{}\mathit{i}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{20}\mathbf{,}\mathbf{20}\mathbf{/}\mathbf{19}\mathbf{\}}\mathbf{]}\mathbf{]}\mathbf{;}$

$\mathbf{\text{solplot}}\mathbf{=}\mathbf{\text{Plot}}\mathbf{[}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{\text{/.}}\mathbf{\text{numsols}}\mathbf{,}\mathbf{\{}\mathit{t}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{10}\mathbf{\}}\mathbf{,}\mathbf{\text{PlotRange}}\mathbf{\to }\mathbf{\text{All}}\mathbf{,}$

$\mathbf{\text{PlotStyle}}\mathbf{\to }\mathbf{\{}\mathbf{\{}\mathbf{\text{CMYKColor}}\mathbf{[}\mathbf{0}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{.}\mathbf{5}\mathbf{]}\mathbf{,}\mathbf{\text{Thickness}}\mathbf{[}\mathbf{.}\mathbf{005}\mathbf{]}\mathbf{\}}\mathbf{\}}\mathbf{]}\mathbf{;}$

Last, we display the direction field along with the solution graphs in solplot using Show in Fig. 6.11 (b).

$\mathbf{\text{Show}}\mathbf{[}\mathbf{\text{GraphicsRow}}\mathbf{[}\mathbf{\{}\mathbf{\text{pvf20}}\mathbf{,}\mathbf{\text{Show}}\mathbf{[}\mathbf{\text{pvf20}}\mathbf{,}\mathbf{\text{solplot}}\mathbf{,}\mathbf{\text{Frame}}\mathbf{\to }\mathbf{\text{False}}\mathbf{]}\mathbf{\}}\mathbf{]}\mathbf{]}$

Notice that there are three nontrivial equilibrium solutions that are found by solving $y^{\prime }=0$ .

$\mathbf{\text{Solve}}\mathbf{[}\mathbf{\text{eqn}}\mathbf{[}\mathbf{20}\mathbf{.}\mathbf{]}\mathbf{[}\mathbf{[}\mathbf{2}\mathbf{]}\mathbf{]}\mathbf{\text{==}}\mathbf{0}\mathbf{,}\mathit{y}\mathbf{[}\mathit{t}\mathbf{]}\mathbf{,}\mathit{t}\mathbf{]}$

$\{\{y[t]\to 0.\},\{y[t]\to 0.926741\},\{y[t]\to 7.38645\},\{y[t]\to 11.6868\}\}$

In this case, $y\approx .926$ and $y\approx 11.687$ are stable, whereas $y\approx 7.386$ is unstable.   □

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128241639000112

Introduction to Differential Equations

Martha L. Abell , James P. Braselton , in Differential Equations with Mathematica (Fourth Edition), 2016

1.3 Initial and Boundary-Value Problems

In many applications, we are not only given a differential equation to solve but we are given one or more conditions that must be satisfied by the solution(s) as well. For example, suppose that we want to find an antiderivative of the function f(x) = 3x ² − 4x. Then, we solve the differential equation dy/dx = 3x ² − 4x by integrating.

$\begin{array}{l}\hfill \frac{dy}{dx}=3x^2-4x⟹y=\int \left(3x^2-4x\right)dx⟹y=x^3-2x^2+C.\end{array}$

step1=Integrate[3x ^∧ 2−4x, x]+c

c − 2x ² + x ³

Because the solution involves an arbitrary constant and all solutions to the equation can be obtained from it, we call this a general solution. On the other hand, if we want to find a solution that passes through the point (1, 4), we must find a solution that satisfies the auxiliary condition y(1) = 4. Substitution into y = x ³ − 2x ² + c yields y(1) = 1³ − 2 ⋅ 1² + c = 4⟹c = 5. Therefore, the member of the family of solutions y = x ³ − 2x ² + c that satisfies y(1) = 4 is y = x ³ − 2x ² + 5. The following commands illustrate how to graph some members of the family of solutions by substituting various values of c into the general solution. We also graph the solution to the problem

$\begin{array}{l}\hfill \left\{\begin{array}{l}{y}^{\prime }=3x^2-4x\\ y(1)=4\end{array}\right..\end{array}$

First, we use Table to generate a table of functions x ³ − 2x ² + c for c = −20, − 18, $\dots$ , 18, 20, naming the resulting set of functions toplot. Note that we use c to avoid conflict with the built-in symbol C. The set of functions toplot is not displayed (for length reasons) because a semi-colon (;) is included at the end of the command.

toplot=Table[step1, {c, −20, 20, 2}];

To graph the functions contained in toplot we use Table. Then, we graph toplot with Plot in Figure 1-6 (a).

Figure 1-6. (a) Plot of y = x ³ − 2x ² + C for various values of C. (b) Plot of the solution that satisfies y(1) = 4

q1=Plot[toplot, {x, −2, 3}, AxesLabel→{x, y}, PlotLabel→ "Various Plots", PlotRange→{−25, 25}];

The solution that satisfies the initial condition is graphed in Figure 1-6 (b).

q2=Plot[x ^∧ 3−2x ^∧ 2+5, {x, −2, 3}, AxesLabel→{x, y}, PlotLabel→"A Particular Solution", PlotRange→{−25, 25}, PlotStyle→Black];

Show[GraphicsRow[{q1, q2}]]

Notice that this first-order equation requires one auxiliary condition to eliminate the unknown coefficient in the general solution. Frequently, the independent variable in a problem is t, which usually represents time. Therefore, we call the auxiliary condition of a first-order equation an initial condition, because it indicates the initial-value (at t = t ₀) of the dependent variable. Problems that involve an initial condition are called initial-value problems.

If DSolve cannot find an exact solution to an initial-value problem or if numerical results are desired, the command NDSolve[{y'[x]==f[x,y[x]],y[x0]==y0},y[x],{x,a,b}]

attempts to find a numerical solution to the initial-value problem $\left\{y^{\prime }=f\left(x,y\right),y\right.$ $\left.\left(x_0\right)=y_0\right\}$ valid for a ≤ x ≤ b. We use the Help Browser to obtain information about NDSolve and its options as well as several examples illustrating its use.

Notice that the syntax of the NDSolve command is almost identical to that of the DSolve command, except what we must specify an interval [a, b] on which we want the numerical solution to be valid.

Example 1.3.1

Graph the solution to the initial-value problem $\left\{\begin{array}{l}{y}^{\prime }=sin\left(x^2\right)\\ y(0)=0\end{array}\right.$ on the interval [0, 10]. Evaluate y(5).

Solution

In this case, we see that DSolve is able to solve the initial-value problem although the result is given in terms of the FresnelS function.

exactsol=DSolve[{y′[x]==Sin[x ^∧ 2], y[0]==0}, y[x], x]

{{y[x]→π2FresnelS[2πx]}}

Here is Mathematica's description of the FresnelS function.

Using NDSolve, we obtain a numerical solution to the initial-value problem valid on the interval [0, 10]:

numsol=NDSolve[{y′[x]==Sin[x ^∧ 2], y[0]==0}, y[x], {x, 0, 10}]

{{y[x]→InterpolatingFunction[][x]}}

which we graph with Plot in Figure 1-7.

Figure 1-7. Plot of the solution of $y^{\prime }=sin\left(x^2\right)$ , y(0) = 0 for 0 ≤ x ≤ 10

Plot[y[x]/.numsol, {x, 0, 10}, PlotRange→All, PlotStyle→Black, AxesStyle→Black, AxesLabel→{x, y}]

The value of y(5) is found with ReplaceAll, \..

numsol/.x→5

{{y[5] → 0.527917}}

and indicates that y(5) ≈ 0.5279.

Because first-order equations involve a single auxiliary condition, which is usually referred to as an initial condition, we use the following examples to distinguish between initial-value and boundary-valueproblems which involve higher-order equations.

Example 1.3.2

Consider the second-order differential equation x″ + x = 0, which models the motion of a mass with m = 1 attached to the end of a spring with spring constant k = 1, where x(t) represents the displacement of the mass from the equilibrium position x = 0 at time t. A general solution of this differential equation is found to be $x(t)=Acost+Bsint$ , where A and B are arbitrary constants, with DSolve.

Clear[x] gensol=DSolve[x"[t]+x[t]==0, x[t], t]

{{x[t] → C[1]Cos[t] + C[2]Sin[t]}}

Because this is a second-order equation, we need two auxiliary conditions to determine the two unknown constants. Suppose that the initial displacement of the mass is x(0) = 0 and the initial velocity is x′(0) = 1. This is an initial-value problem because we have two auxiliary conditions given at the same value of t, namely t = 0. Use these initial conditions to determine the solution of this problem.

Solution

Because we need the first derivative of the general solution, we calculate $x^{\prime }(t)=Bcost-Asint$ . Substitution yields x(0) = A = 0 and x′(0) = B = 1. Hence, the solution is $x(t)=sint$ . DSolve can solve this initial-value problem as well.

Clear[x]toplot2=DSolve[{x"[t]+x[t] ==0 , x[0] ==0 , x′[0] ==1}, x[t], t]

{{x[t] →Sin[t]}}

We plot various solutions to the differential equation by generating a list of functions in toplot1 obtained by replacing C[1] and C[2] in gensol with specific numbers.

toplot1=Table[gensol[[1, 1, 2]]/.{C[1]→ i, C[2]→ j}, {i, −2, 2}, {j, −2, 2}];

We then plot the functions in toplot1 with Plot in q1

q1=Plot[toplot1, {t, 0, 4Pi}, AxesStyle→Black, AxesLabel→{x, y}, PlotLabel→"Plot 1"];

and the solution to the initial value problem in q2.

q2=Plot[toplot2[[1, 1, 2]], {t, 0, 4Pi}, AxesStyle→Black, AxesLabel→{x, y}, PlotLabel→"Plot 2"];

We then show the two graphs side-by-side using Show together with GraphicsRow in Figure 1-8

Figure 1-8. In Plot 1, various solutions of the second-order differential equation. In Plot 2, the specific function that satisfies the initial conditions x(0) = 0 and x′(0) = 1

Show[GraphicsRow[{q1, q2}]]

Example 1.3.3

The shape of a bendable beam of length 1 unit that is subjected to a compressive force at one end is described by the graph of the solution y(x) of the differential equation $\frac{d^{2}y}{d{x}^2}+\frac{{\pi }^2}{4}y=0$ , 0 < x < 1. If the height of the beam above the x-axis is known at the endpoints x = 0 and x = 1, then we have a boundary-value problem. Use the boundary conditions y(0) = 1 and y(1) = 2 to find the shape of the beam.

Solution

First, we use DSolve to find a general solution to the equation. The result indicates that a general solution is $y(x)=Acos\left(\frac{\pi }{2}x\right)+Bsin\left(\frac{\pi }{2}x\right)$ .

Later, we will see that under reasonable conditions, initial-value problems have unique solutions. On the other hand, boundary-value problems may have no solutions, infinitely many solutions or a unique solution.

DSolve[y"[x]+Pi ^∧ 2/4y[x]==0, y[x], x]

$\left\{\left\{y[x]\to C[1]\text{Cos}\left[\frac{\pi x}{2}\right]+C[2]\text{Sin}\left[\frac{\pi x}{2}\right]\right\}\right\}$

Applying the condition y(0) = 0 to the general solution yields

$\begin{array}{l}\hfill y(0)=Acos0+Bsin0=A=0.\end{array}$

Similarly, y(1) = 2 indicates that

$\begin{array}{l}\hfill y(1)=Bsin\frac{\pi }{2}=B=2,\end{array}$

so the solution to the boundary value problem is $y(x)=2sin\left(\frac{\pi }{2}x\right)$ , 0 < x < 1. DSolve is also able to solve this boundary value problem.

sol=DSolve[{y"[x]+Pi ^∧ 2/4y[x]==0, y[0]==0, y[1]==2}, y[x], x]

$\left\{\left\{y[x]\to 2\text{Sin}\left[\frac{\pi x}{2}\right]\right\}\right\}$

This function that describes the shape of the beam is graphed with Plot in Figure 1-9.

Figure 1-9. A plot of a solution to the beam equation

Plot[y[x]/.sol, {x, 0, 1}, PlotStyle→Black, AxesStyle→Black, AxesLabel→{x, y}]

We will see that it is usually impossible to find exact solutions of higher-order nonlinear initial-value problems. In those cases, we can often use NDSolve to generate an accurate approximation of the solution.

Example 1.3.4 (Rayleigh's Equation)

Rayleigh's equation is the nonlinear equation

Also see Example 1.4.4.

(1.7) $\begin{array}{l}\hfill \frac{d^{2}x}{d{t}^2}+\left[\frac{1}{3}{\left(\frac{dx}{dt}\right)}^2-1\right]\frac{dx}{dt}+x=0\end{array}$

and arises in the study of the motion of a violin string. Graph the solution to Rayleigh's equation on the interval [0, 15] if (a) x(0) = 1, x′(0) = 0; (b) x(0) = 0.1, x′(0) = 0; and (c) x(0) = 0, x′(0) = 1.9.

Solution

First, observer that DSolve cannot solve the differential equation.

DSolve[x"[t]+(1/3x′[t] ^∧ 2−1)x′[t]+x[t]==0, x[t], t]

DSolve [x[t]+x′[t](−1+13x′[t]²)+x″[t]==0, x[t], t]

In each case, we use NDSolve to approximate the solution to the initial-value problem, naming the results numsol1, numsol2, and numsol3, respectively.

numsol1=NDSolve[{x"[t]+(1/3x′[t] ^∧ 2−1)x′[t]+x[t]==0, x[0]==1, x′[0]==0}, x[t], {t, 0, 15}]

{{x[t]→InterpolatingFunction[][t]}}

numsol2=NDSolve[{x"[t]+(1/3x′[t] ^∧ 2−1)x′[t]+x[t]==0, x[0]==.1, x′[0]==0}, x[t], {t, 0, 15}]

{{x[t]→InterpolatingFunction[][t]}}

numsol3=NDSolve[{x"[t]+(1/3x′[t] ^∧ 2−1)x′[t]+x[t]==0, x[0]==0, x′[0]==1.9}, x[t], {t, 0, 15}]

{{x[t]→InterpolatingFunction[][t]}}

All three solutions are graphed together on the interval [0, 15] with Plot in Figure 1-10. Notice that the solution to (c) appears to be periodic.

Figure 1-10. Plot of three solutions of Rayleigh's equation (equation (1.7))

Plot[Evaluate[x[t]/.{numsol1, numsol2, numsol3}], {t, 0, 15}, PlotStyle→{GrayLevel[0], GrayLevel[.3], GrayLevel[.6]}, PlotRange→All, AxesLabel→{x, y}, AxesStyle→Black, PlotLabel→"Solutions to Rayleigh's Equation"]

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128047767000012

Systems of Ordinary Differential Equations

Martha L. Abell , James P. Braselton , in Differential Equations with Mathematica (Fourth Edition), 2016

6.3.3 Alternate Method for Solving Initial-Value Problems

An alternate method can be used to solve initial-value problems.

Let Φ(t) be a fundamental matrix for the system of equations X′(t) = A(t)X(t). Then, a general solution is X′(t) = Φ(t)C, where C is a constant vector. If the initial condition X(0) =X ₀ is given, then

$\begin{array}{l}\hfill \begin{array}{cc}\mathbf{X}(0)& =\mathbf{\Phi }(0)\mathbf{C}\\ {\mathbf{X}}_0& =\mathbf{\Phi }(0)\mathbf{C}\\ \mathbf{C}& ={\mathbf{\Phi }}^{-1}(0){\mathbf{X}}_0.\end{array}\end{array}$

Therefore, the solution to the initial-value problem $\left\{\begin{array}{l}{\mathbf{X}}^{\prime }(t)=\mathbf{A}(t)\mathbf{X}(t)\\ \mathbf{X}(0)={\mathbf{X}}_0\end{array}\right.$ is X(t) = Φ(t)Φ ⁻¹(0)X ₀.

Example 6.3.4

Use a fundamental matrix to solve the initial-value problem ${\mathbf{X}}^{\prime }=\left(\begin{array}{ll}\hfill 1\hfill & \hfill 1\hfill \\ \hfill 4\hfill & \hfill -2\hfill \end{array}\right)\mathbf{X}$ subject to $\mathbf{X}(0)=\left(\begin{array}{l}\hfill 1\hfill \\ \hfill -2\hfill \end{array}\right)$ .

Solution

We first remark that you can use DSolve to solve the initial-value problem directly with the command

Clear[x, y]DSolve[{x′[t]==x[t]+y[t], y′[t]==4x[t]−2y[t], x[0]==1, y[0]==−2}, {x[t], y[t]}, t]

$\{\{x[t]\to \frac{1}{5}{e}^{-3t}(3+2e^{5t}),y[t]\to \frac{2}{5}{e}^{-3t}(-6+e^{5t})\}\}$

The eigenvalues and corresponding eigenvectors of $\mathbf{A}=\left(\begin{array}{ll}\hfill 1\hfill & \hfill 1\hfill \\ \hfill 4\hfill & \hfill -2\hfill \end{array}\right)$ are found with Eigensystem.

a={{1, 1}, {4, −2}}; s1=Eigensystem[a]

{{−3, 2}, {{−1, 4}, {1, 1}}}

The results mean that the eigenvalues are λ ₁ = −3 and λ ₂ = 2 with corresponding eigenvectors ${\mathbf{v}}_1=\left(\begin{array}{l}\hfill -1\hfill \\ \hfill 4\hfill \end{array}\right)$ and ${\mathbf{v}}_1=\left(\begin{array}{l}\hfill 1\hfill \\ \hfill 1\hfill \end{array}\right)$ , respectively. A fundamental matrix is then given by $\mathbf{\Phi }(t)=\left(\begin{array}{ll}\hfill -e^{3t}\hfill & \hfill e^{2t}\hfill \\ \hfill 4e^{3t}\hfill & \hfill e^{2t}\hfill \end{array}\right)$ .

phi[t_]={s1[[2, 1]]Exp[−3t], s1[[2, 2]]Exp[2t]}//Transpose; MatrixForm[phi[t]]

$\left(\begin{array}{ll}-e^{-3t}& e^{2t}\\ 4e^{-3t}& e^{2t}\\ \end{array}\right)$

We calculate Φ ⁻¹(0) with Inverse.

Inverse[A] finds the inverse of the square matrix A, if A is invertible.

Inverse[phi[t]]//MatrixForm

$\left(\begin{array}{ll}-\frac{e^{3t}}{5}& \frac{e^{3t}}{5}\\ \frac{4e^{-2t}}{5}& \frac{e^{-2t}}{5}\\ \end{array}\right)$

Hence, the solution to the initial-value problem is X(t) = Φ(t)Φ ⁻¹(0)X ₀.

sol=phi[t].Inverse[phi[0]].{1, −2}// Simplify

$\{\frac{1}{5}{e}^{-3t}(3+2e^{5t}),\frac{2}{5}{e}^{-3t}(-6+e^{5t})\}$

As in the previous examples, we graph x(t) and y(t) together in Figure 6-17 (a) and parametrically in Figure 6-17 (b).

Figure 6-17. (a) x(t) and y(t). (b) Parametric plot of x(t) versus y(t)

p1=Plot[Evaluate[{x[t], y[t]}/.{x[t]-> sol[[1]], y[t]→sol[[2]]}], {t, −1, 3}, PlotRange→{{−1, 3}, {−2, 2}}, AspectRatio→1]

p2=ParametricPlot[Evaluate[{x[t], y[t]}/.{x[t]-> sol[[1]], y[t]→sol[[2]]}], {t, −1, 3}, PlotRange→{{0, 20}, {−10, 10}}, AspectRatio→1]

Show[GraphicsRow[{p1, p2}]]

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128047767000061

Flow reactors

Henry C. Foley , in Introduction to Chemical Engineering Analysis Using Mathematica (Second Edition), 2021

Solution of the steady state PFR

Firstly, we begin by solving for the concentrations with linear kinetics, and we so do in complete form. We will do this with DSolve as an exercise, even though the equations are trivial to solve:

$\left\{\{\text{Ca}[z]\to \text{Caf}{e}^{-\frac{\text{kab}z}{\text{vz}}},\text{Cb}[z]\to \text{Caf}-\text{Caf}{e}^{-\frac{\text{kab}z}{\text{vz}}}\}\right\}$

Recall that kab is just $\frac{1}{{\tau }_{\mathbf{\text{rxn}}}}$ , that is the reciprocal of the characteristic time for reaction and $\frac{z}{v_z}$ is the holding time τ for the PFR. Therefore these solutions become:

$$\begin{gathered} \text{Ca}[z]=\text{Caf}{e}^{-\frac{\text{kab}z}{\text{vz}}}⟹\text{undefined}\text{Ca}[\tau ]=\text{Caf}{e}^{-\text{kab}\tau } \\ \text{Ca}[\tau ]=\text{Caf}{e}^{-\frac{\tau }{{\tau }_{\text{rxn}}}} \\ \text{Cb}[z]=\text{Caf}(1-e^{-\frac{\text{kab}z}{\text{vz}}})⟹\text{undefined}\text{Cb}[\tau ]=\text{Caf}(1-e^{-\text{kab}\tau }) \\ \text{Cb}[\tau ]=\text{Caf}(1-e^{-\frac{\tau }{{\tau }_{\text{rxn}}}}) \end{gathered}$$

When we solved the transient, well-mixed batch reactor with linear kinetics, we obtained the same solution functionally, but instead of kab τ, we had kab t as the argument of the exponential, that is in terms of real time instead of holding time.

We return now to the Langmuir-Hinshelwood kinetics from the CSTR section to see how the PFR will behave and to compare the CSTR and the PFR. As in the case of the steady state CSTR, we will write a steady state PFR Module function. Recall that the rate law was:

$r_{A-}\text{ =}\frac{k{K}_A{C}_A}{1+K_A{C}_A}$

Therefore, once again invoking the pseudo-homogeneous approximation, the equations that we must solve are:

$$\begin{gathered} v_z\text{undefined}\frac{\partial C_a}{\partial z}\text{= -}\frac{(1-ϵ)}{ϵ}\frac{k{K}_A{C}_A}{1+K_A{C}_A} \\ {v}_z\text{undefined}\frac{\partial C_b}{\partial z}\text{= +}\frac{(1-ϵ)}{ϵ}\frac{k{K}_A{C}_A}{1+K_A{C}_A} \end{gathered}$$

The Module function for this PFR will be called "pfr1." The arguments in "pfr1" are the rate constant k, the adsorption equilibrium constant K1, the volume flow rate q, and the radius of the reactor cross-section, r. So that we can make comparisons to the CSTR, we have kept the total volume the same at 1000. Once we specify the radius, that fixes the circular cross-section. This divided into the volume gives the length of the reactor, zmax. Therefore instead of specifying the reactor length, we simply specify the reactor radius and at constant volume this fixes zmax. Everything else will be kept the same for comparison sake, especially ϵ and the holding time $\frac{V}{q}$ .

To see how the program works let's try it out for a radius of 10 units and a volume of 10,000. This makes zmax about 32 units.

Now, for comparison to the CSTR, we can build a new Module function that allows us to create a plot of the exit concentration from a CSTR under the same conditions and also as a "function" of z. Of course there is no functional z-dependence for the CSTR, and the plots will be simply horizontal lines. But, they will be graphed over the same range as the axial distance through the PFR. In this way we can put the two on one graph for comparison. Here is the steady state CSTR Module:

Putting the PFR and CSTR results on the same graph we find:

We see that conversion at the exit of the PFR is larger than that of the CSTR at the same condition. This is because there is no back mixing in the PFR.

We can vary the flow rate over a few orders of magnitude for both the PRF and the CSTR to see what will happen. This is done by combining the two functions into one Table and the plotting the results as a GraphicsArray.

For the PFR at shorter q values the cross-over in concentrations moves to shorter holding times. For the CSTR the trend is the same, as q rises the conversion drops. In fact, we can find that even at q = 1, the CSTR has still not achieved full conversion of the feedstock. Remember the only difference between the two cases is that the CSTR is well-mixed throughout its volume, but the PFR is well-mixed only radially, and not at all axially. So this leads to the oft quoted rule-of-thumb that the PFR is more efficient than the CSTR. At the same volume of reactor, one achieves higher conversion, or to achieve equivalent conversion the PFR volume can be smaller than the CSTR volume.

Now in the previous calculations, we assumed that K1 was small. This forced the rate toward first order dependence. Therefore how does the comparison between PFR and CSTR work out if we vary K1 over several orders of magnitude to make the kinetics range from first order to zeroth order? To do this we will fix the flow rate at q =10 and vary K1 from ${10}^{-3}$ to ${10}^0.$

The results indicate that even if the adsorption equilibrium constant is increased, the PFR shows better results than the CSTR until the very last value at K1 = 0.1. The only time this is violated for an isothermal reaction system is if the kinetics are negative order overall, then the CSTR will actually be more efficient than the PFR, otherwise the PFR is better. Ceteris paribus, the PFR will give the same conversion but with a smaller volume, or higher conversion at the same volume.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128200513000139

Applications of the Symbolic Toolbox

G.R. Lindfield , J.E.T. Penny , in Numerical Methods (Third Edition), 2012

9.11 Symbolic Solution of Ordinary Differential Equations

The Symbolic Toolbox can be used to solve, symbolically, first-order differential equations, systems of first-order differential equations, or higher-order differential equations, together with any initial conditions provided. This symbolic solution of differential equations is implemented in Matlab using the function dsolve , and its use is illustrated by a range of examples.

It is important to note that this approach only provides a symbolic solution if one exists. If no solution exists, then the user should apply one of the numerical techniques provided in Matlab, such as ode45.

The general form of a call of the function dsolve to solve a differential equation system is

sol = dsolve('de1, de2, de3, … , den, in1, in2, in3, … , inn');

The independent variable is assumed to be t unless given by an optional final parameter of dsolve. The parameters de1, de2, de3 up to den stand for the individual differential equations. These must be written in symbolic form using symbolic variables, standard Matlab operators, and the symbols D, D2, D3, and so on, which represent the first-, second-, third-, and higher-order differential operators, respectively. The parameters in1, in2, in3, in4, and so on, represent the initial conditions for the differential equations, if these conditions are required. An example of how these initial conditions should be written, assuming a dependent variable y, is

y(0) = 1,   Dy(0) = 0,   D2y(0) = 9.1

which means that the value of y is 1, dy/dt = 0, and d ² y/dt ² = 9.1 when t = 0. It is important to note that dsolve accepts up to a maximum of 12 input parameters. If initial conditions are required, this is a significant restriction!

The solution is returned to sol as a Matlab structure, and consequently the names of the dependent variables must be used to indicate the individual components. For example, if g and y are two dependent variables for the differential equation, sol.y gives the solution for the dependent variable y and sol.g gives the solution for the dependent variable g.

To illustrate these points we consider some examples. Consider the following first-order differential equation:

$\left(1+t^2\right)\frac{dy}{dt}+2ty=\cos t$

We may solve this without initial conditions using dsolve as follows:

>> s = dsolve('(1+t^2)*Dy+2*t*y=cos(t)')

s =

-(C3 - sin(t))/(t^2 + 1)

Notice that the solution contains the arbitrary constant C3. If we now solve the same equation using an initial condition, we proceed as follows:

>> s = dsolve('(1+t^2)*Dy+2*t*y=cos(t),y(0)=0')

s =

sin(t)/(t^2 + 1)

Note that in this case there is no arbitrary constant.

We now solve a second-order system

$\frac{d^{2}y}{d{x}^2}+y=\cos 2x$

with the initial conditions y = 0 and dy/dx = 1 at x = 0. To solve this differential equation dsolve has the form

>> dsolve('D2y+y=cos(2*x), Dy(0)=1, y(0)=0','x')

ans =

(2*cos(x))/3 + sin(x) + sin(x)*(sin(3*x)/6 + sin(x)/2) -

  (2*cos(x)*(6*tan(x/2)^2 - 3*tan(x/2)^4 + 1))/(3*(tan(x/2)^2 + 1)^3)

>> simplify(ans)

ans =

sin(x) + (2*sin(x)^2)/3 - (2*sin(x/2)^2)/3

Notice that since the independent variable is x, this is indicated in dsolve by a final parameter x in the list of parameters.

We now solve the fourth-order differential equation

$\frac{d^{4}y}{d{t}^4}=y$

with initial conditions

$\begin{array}{ccc}y=1,dy/dt=0,d^{2}y/d{t}^2=-1,d^{3}y/d{t}^3=0& \text{when}& t=\pi /2\end{array}$

We again use dsolve. In this example D4 stands for the fourth derivative operator, with respect to t, and so on.

>> dsolve('D4y=y, y(pi/2)=1, Dy(pi/2)=0, D2y(pi/2)=-1, D3y(pi/2)=0')

ans =

sin(t)

However, we note that if we try to solve the apparently simple problem

$\frac{dy}{dx}=\frac{e^{-x}}{x}$

with the initial condition y = 1 when x = 1, difficulties arise. Applying dsolve, we have

>> dsolve('Dy=exp(-x)/x, y(1)=1', 'x')

ans =

1 - Ei(1, x) - Ei(-1)

Note that Ei(-1) = -Ei(1,1). Clearly this result is not an explicit solution. The function Ei(1,x) is the exponential integral and can be found in mfunlist. Details of this mathematical function can be found in Abramowitz and Stegun (1965) and Olver et al. (2010). We can evaluate the function Ei using the mfun function for any parameters. For example

>> y = mfun('Ei',1,1)

y =

  0.2194

If we require more digits in the solution, we have

vpa('Ei(1,1)',20)

ans =

0.21938393439552027368

We will now attempt to solve the following apparently simple differential equation:

$\frac{dy}{dx}=\cos \left(\sin x\right)$

Applying dsolve we have

>> dsolve('Dy=cos(sin(x))','x')

ans =

C17 + int(cos(sin(x)), x, IgnoreAnalyticConstraints)

In this case dsolve fails to solve the equation.

The differential equation may also contain constants represented as symbols. For example, if we wish to solve the equation

$\frac{d^{2}x}{d{t}^2}+\frac{a}{b}\sin t=0$

with initial conditions x = 1 and dx/dt = 0 when t = 0, we enter the following:

>> syms x t a b

>> x = dsolve('D2x+(a/b)*sin(t)=0,x(0)=1,Dx(0)=0')

x =

(a*sin(t))/b - (a*t)/b + 1

Note how the variables a and b appear in the solution as expected.

As an example of solving two simultaneous differential equations, we note that this differential equation may be rewritten as

$\begin{array}{lll}\hfill \frac{du}{dt}=-\frac{a}{b}\sin t& & \hfill \\ \hfill \frac{dx}{dt}=u& & \hfill \end{array}$

Using the same initial conditions, dsolve may be applied to solve these equations by writing

>> syms u

>> [u x] = dsolve('Du+(a/b)*sin(t)=0,Dx=u,x(0)=1,u(0)=0')

u =

(a*cos(t))/b - a/b

x =

(a*sin(t))/b - (a*t)/b + 1

This gives the same solution as that obtained from dsolve applied directly to solve the second-order differential equation.

The following example provides an interesting comparison of the symbolic and numerical approach. It consists of a script and the output from the script. The script compares the use of dsolve for the symbolic solution of a differential equation and the use ode45 for the numerical solution of the same differential equation. Note that the symbolic solution is obtained in two ways: by solving the second-order equation directly and by separating it into two first-order simultaneous equations. Both approaches provide the same solution.

% e3s904.m   Simultaneous first order differential equations

% dx/dt = y, Dy = 3*t-4*x.

% Using dsolve this become

syms y t x

x = dsolve('D2x+4*x=3*t','x(0)=0', 'Dx(0)=1')

tt = 0:0.1:5; p = subs(x,t,tt); pp = double(p);

% Plot the symbolic solution to the differential equ'n

plot(tt,pp,'r')

hold on

xlabel('t'), ylabel('x')

sol = dsolve('Dx=y','Dy=3*t-4*x', 'x(0)=0', 'y(0)=1');

sol_x = sol.x, sol_y = sol.y

fv = @(t,x) [x(2); 3*t-4*x(1)];

options = odeset('reltol', 1e-5,'abstol',1e-5);

tspan = [0 5]; initx = [0 1];

[t,x] = ode45(fv,tspan,initx,options);

plot(t,x(:,1),'k+');

axis([0 5 0 4])

Executing the script gives the symbolic solution

x =

  (3*t)/4 + sin(2*t)/8

sol_x =

(3*t)/4 + sin(2*t)/8

sol_y =

cos(2*t)/4 + 3/4

This script also provides the graph of the symbolic solution with alternate numerical solution values also plotted (Figure 9.3). Note how the numerical solution is consistent with the symbolic one.

Figure 9.3. Symbolic solution and numerical solution indicated by +.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780123869425000096

Applications of the Symbolic Toolbox

George Lindfield , John Penny , in Numerical Methods (Fourth Edition), 2019

10.11 Symbolic Solution of Ordinary Differential Equations

The Symbolic Toolbox can be used to solve, symbolically, first-order differential equations, systems of first-order differential equations or higher-order differential equations, together with any initial conditions provided. This symbolic solution of differential equations is implemented in Matlab using the function dsolve , and its use is illustrated by a range of examples.

It is important to note that this approach only provides a symbolic solution if one exists. If no solution exists, then the user should apply one of the numerical techniques provided in Matlab, such as ode45.

The general form of a call of the function dsolve to solve a differential equation system is

sol = dsolve('de1, de2, de3, ... , den, in1, in2, in3, ... , inn');

The independent variable is assumed to be t unless given by an optional final parameter of dsolve. The parameters de1, de2, de3 up to den stand for the individual differential equations. These must be written in symbolic form using symbolic variables, standard Matlab operators and the symbols D, D2, D3, D4, ..., etc., which represent the first-, second-, third-, and higher-order differential operators respectively. The parameters in1, in2, in3, in4, and so on, represent the initial conditions for the differential equations, if these conditions are required. An example of how these initial conditions should be written, assuming an dependent variable y, is

y(0) = 1,   Dy(0) = 0,   D2y(0) = 9.1

which means that the value of y is 1, $dy/dt=0$ and $d^{2}y/d{t}^2=9.1$ when $t=0$ .

The solution is returned to sol as a Matlab structure, and consequently the names of the dependent variables must be used to indicate the individual components. For example, if g and y are two dependent variables for the differential equation, sol.y gives the solution for the dependent variable y and sol.g gives the solution for the dependent variable g.

To illustrate these points we consider some examples. Consider a first-order differential equation:

$(1+t^2)\frac{dy}{dt}+2ty=\cos t$

We may solve this without initial conditions using dsolve as follows:

>> s = dsolve('(1+t^2)*Dy+2*t*y=cos(t)')

s =

-(C4 - sin(t))/(t^2 + 1)

Notice that the solution contains the arbitrary constant C4. If we now solve the same equation using an initial condition, we proceed as follows:

>> s = dsolve('(1+t^2)*Dy+2*t*y=cos(t),y(0)=0')

s =

sin(t)/(t^2 + 1)

Note that in this case there is no arbitrary constant.

We now solve a second-order system

$\frac{d^{2}y}{d{x}^2}+y=\cos 2x$

with the initial conditions $y=0$ and $dy/dx=1$ at $x=0$ . To solve this differential equation dsolve has the form below, where x is the independent variable.

>> dsolve('D2y+y=cos(2*x), Dy(0)=1, y(0)=0','x')

ans =

(2*cos(x))/3 + sin(x) + sin(x)*(sin(3*x)/6 + sin(x)/2) -

(2*cos(x)*(6*tan(x/2)^2 - 3*tan(x/2)^4 + 1))/(3*(tan(x/2)^2 + 1)^3)

>> simplify(ans)

ans =

(5*cos(x))/3 + sin(x) - (8*cos(x/2)^4)/3 + 1

We now solve a fourth-order differential equation:

$\frac{d^{4}y}{d{t}^4}=y$

with initial conditions

$y=1,dy/dt=0,d^{2}y/d{t}^2=-1,d^{3}y/d{t}^3=0\text{ when}t=\pi /2.$

We again use dsolve. In this example D4 stands for the fourth derivative operator, with respect to t, and so on.

>> dsolve('D4y=y, y(pi/2)=1, Dy(pi/2)=0, D2y(pi/2)=-1, D3y(pi/2)=0')

ans =

sin(t)

However, we note that if we try to solve the apparently simple problem

$\frac{dy}{dx}=\frac{e^{-x}}{x}$

with the initial condition $y=1$ when $x=1$ , difficulties arise. Applying dsolve, we have

>> dsolve('Dy=exp(-x)/x, y(1)=1', 'x')

ans =

1 - expint(x) - ei(-1)

The function expint(z), where z is real or complex, is the exponential integral function defined by

$E_1(z)=-{\int }_z^{\infty }\frac{\exp (-t)}{t}dtz\ne 0$

$E_1(z)$ assumes its principal value. In fact $E_1(z)$ is a particular case of the more general $E_n(z)$ function.

Details of this mathematical function can be found in Abramowitz and Stegun (1965) and Olver et al. (2010). The function ei(z) is the one-argument exponential integral function and is defined, for $x>0$ by

$\text{Ei}(x)=-{\int }_{-x}^{\infty }\frac{\exp (-t)}{t}dt={\int }_{-\infty }^x\frac{\exp (t)}{t}dtx>0$

The functions $E_1$ and Ei are related. Thus

$E_1(x)={\int }_x^{\infty }\frac{\exp (-t)}{t}dt=-\text{Ei}(-x)x>0$

and

$E_1(-x)=-\text{Ei}(x)-ȷ\pi$

We can illustrate these relationships. For example, if $x=3$ ,

>> x = 3;

>> [expint(x) -ei(-x)]

ans =

0.0130   0.0130

>> [expint(-x) -ei(x)-j*pi]

ans =

-9.9338 - 3.1416i   -9.9338 - 3.1416i

Returning to the solution of the differential equation, we can evaluate the ei(-1) as follows:

>> ei(-1)

ans =

-0.2194

Thus, for example, when $x=2$

>> x = 2;

>> 1 - expint(x) - ei(-1)

ans =

1.1705

We will now attempt to solve the following apparently simple differential equation:

$\frac{dy}{dx}=\cos (\sin x)$

Applying dsolve we have

>> dsolve('Dy=cos(sin(x))','x')

ans =

C14 + int(cos(sin(x)), x, 'IgnoreSpecialCases', true,

'IgnoreAnalyticConstraints', true)

In this case dsolve fails to solve the equation.

The differential equation may also contain constants represented as symbols. For example, if we wish to solve the equation

$\frac{d^{2}x}{d{t}^2}+\frac{a}{b}\sin t=0$

with initial conditions $x=1$ and $dx/dt=0$ when $t=0$ . Thus

>> syms x t a b

>> x = dsolve('D2x+(a/b)*sin(t)=0,x(0)=1,Dx(0)=0')

x =

(a*sin(t))/b - (a*t)/b + 1

Note how the variables a and b appear in the solution as expected.

As an example of solving two simultaneous differential equations, we note that this differential equation may be rewritten as

$\begin{array}{c}\frac{du}{dt}=-\frac{a}{b}\sin t\hfill \\ \frac{dx}{dt}=u\hfill \end{array}$

Using the same initial conditions, dsolve may be applied to solve these equations by writing

>> syms u

>> [u x] = dsolve('Du+(a/b)*sin(t)=0,Dx=u,x(0)=1,u(0)=0')

u =

(a*cos(t))/b - a/b

x =

(a*sin(t))/b - (a*t)/b + 1

The same solution as that obtained from dsolve applied directly to solve the second-order differential equation.

The script e4sX04.m provides an interesting comparison between the symbolic and numerical approach. It consists of a script and the output from the script. The script compares the use of dsolve for the symbolic solution of a differential equation, and the use of ode45 for the numerical solution of the same differential equation. Note that the symbolic solution is obtained in two ways: by solving the second-order equation directly and by separating it into two first-order simultaneous equations. Both approaches provide the same solution.

% e4sX04.m   Simultaneous first order differential equations

% dx/dt = y, Dy = 3*t-4*x.

% Using dsolve this becomes

syms y t x

x = dsolve('D2x+4*x=3*t','x(0)=0', 'Dx(0)=1')

tt = 0:0.1:5; p = subs(x,t,tt); pp = double(p);

% Plot the symbolic solution to the differential equ'n

plot(tt,pp,'r')

hold on

xlabel('t'), ylabel('x')

sol = dsolve('Dx=y','Dy=3*t-4*x', 'x(0)=0', 'y(0)=1');

sol_x = sol.x, sol_y = sol.y

fv = @(t,x) [x(2); 3*t-4*x(1)];

options = odeset('reltol', 1e-5,'abstol',1e-5);

tspan = [0 5]; initx = [0 1];

[t,x] = ode45(fv,tspan,initx,options);

plot(t,x(:,1),'k+');

axis([0 5 0 4])

Executing script e4sX04.m gives the symbolic solution

x =

(3*t)/4 + sin(2*t)/8

sol_x =

(3*t)/4 + sin(2*t)/8

sol_y =

cos(2*t)/4 + 3/4

This script also provides the graph of the symbolic solution with alternate numerical solution values also plotted, Fig. 10.2. Note how the numerical solution is consistent with the symbolic one.

Figure 10.2. Symbolic solution and numerical solution indicated by +.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128122563000191

Applications of Systems of Ordinary Differential Equations

Martha L. Abell , James P. Braselton , in Differential Equations with Mathematica (Fourth Edition), 2016

7.1.4 Spring-Mass Systems

The displacement of a mass attached to the end of a spring was modeled with a second-order linear differential equation with constant coefficients in Chapter 5. This situation can then be expressed as a system of first-order ordinary differential equations as well. Recall that if there is no external forcing function, then the second-order differential equation that models this situation is mx″ + cx′ + kx = 0, where m is the mass attached to the end of the spring, c is the damping coefficient, and k is the spring constant found with Hooke's law. This equation is transformed into a system of equations by letting x′ = y so that $y^{\prime }=x^{″}=-\frac{k}{m}x-\frac{c}{m}{x}^{\prime }$ and then solving the differential equation for x″. After substitution, we have the system

(7.4) $\begin{array}{l}\hfill \left\{\begin{array}{c}\frac{dx}{dt}=y\\ \frac{dy}{dt}=-\frac{k}{m}x-\frac{c}{m}y\end{array}\right..\end{array}$

In previous chapters, the displacement of the spring was illustrated as a function of time. However, problems of this type may also be investigated using the phase plane.

Example 7.1.3

Solve the system of differential equations to find the displacement of the mass if m = 1, c = 0, and k = 1.

Solution

In this case, the system is $\left\{\begin{array}{c}dx/dt=y\\ dy/dt=-x\end{array}\right.$ , which in matrix form is ${\mathbf{X}}^{\prime }=\left(\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 0\hfill \end{array}\right)\mathbf{X}$ . The solution that satisfies the initial conditions x(0) = x ₀ and y(0) = y ₀ is found with DSolve and named gensol for later use.

Clear[x, y]gensol=DSolve[{D[x[t], t]==y[t], D[y[t], t]==−x[t], x[0]==x0, y[0]==y0}, {x[t], y[t]}, t]

{{x[t] →x0Cos[t] + y0Sin[t], y[t] →y0Cos[t] −x0Sin[t]}}

Note that this system is equivalent to the second-order differential equation x″ + x = 0, which we solved in Chapters 4 and 5. At that time, we found a general solution to be $x(t)=c_{1}cost+c_{2}sint$ which is equivalent to the first component of $\mathbf{X}=\left(\begin{array}{c}\hfill x(t)\hfill \\ \hfill y(t)\hfill \end{array}\right)$ , the result obtained with DSolve. Also notice that (0, 0) is the equilibrium point of the system. The eigenvalues of $\mathbf{A}=\left(\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 0\hfill \end{array}\right)$ are λ = ±i,

$\text{Eigenvalues}\left[\left(\begin{array}{cc}0& 1\\ -1& 0\\ \end{array}\right)\right]$

{i, −i}

so we classify the origin as a center.

We graph several members of the phase plane for this system with ParametricPlot in Figure 7-5.

Figure 7-5. The origin is a center

toplot=Table[{x[t], y[t]}/.gensol/.{x0→i, y0→0}, {i, 0, 4, 4/9}];

ParametricPlot[Evaluate[toplot], {t, 0, 2π}, PlotRange→ {{−5, 5}, {−5, 5}}, AspectRatio→1, AxesLabel→{x, y}]

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128047767000073

Ordinary Differential Equations

Frank E. Harris , in Mathematics for Physical Science and Engineering, 2014

Mathematica

It may be convenient to define a variable and assign to it the entire ODE to be solved. If, for example, our ODE is $y^{\prime }(x)+2y(x)=x^2$ , we might define

ode = y'[x] + 2*y[x] == xˆ2

Derivatives are indicated by primes (right single quotes); the second derivative of $y$ is y'', written using two single quotes, not a double quote. Derivatives can also (with identical results) be written using mathematica's D operator. When the above command is executed it is necessary that both $x$ and $y$ be undefined. Note also that the operator = causes everything to its right to be the object assigned to ode and that the equals sign of the ODE is represented by ==.

The command for solving an ODE is DSolve , which must contain as arguments the ODE to be solved, the dependent variable, and the independent variable. Thus, we might write

DSolve[ode, y[x], x]      $\left\{\left\{y[x]\to \frac{1}{4}(1-2x+2x^2)+e^{-2x}C[1]\right\}\right\}$

Note that both in ode and the DSolve command, $y$ is always written with its argu-ment: y[x]. Failure to include arguments in a consistent manner is regarded by mathematica as an error. mathematica does permit the argument to be omitted from y[x] and its derivatives (both in ode and DSolve), but doing so produces different output than we are showing here.

The output identifies the general solution to the ODE, containing one undetermined constant, given the name C[1]. If all we wanted was the solution we are done. But if we want to do further work with the solution, we need to retrieve it in a more usable form. Since each set of braces defines an ordered list, our solution is the first (and only) element of the first (and only) element of the output nested lists, and it therefore can be obtained as

%[[1]][[1]]        $y[x]\to \frac{1}{4}(1-2x+2x^2)+e^{-2x}C[1]$

To recover the solution from this mathematica construct we need to know that it is represented as a two-element list, so our solution can be isolated as

Y = %[[2]]        $\frac{1}{4}(1-2x+2x^2)+e^{-2x}C[1]$

This expression, which we named Y, can be processed by the usual mathematica commands (including those for differentiation and integration), but it has not been defined as a function and notations like Y' have no meaning. We can create a function $y(x)$ from Y by writing

y[x_] ≔ Evaluate[Y]

The command Evaluate is needed because mathematica does not substitute a value for Y before analyzing its dependence on $x$ .

The solution, in either of the forms $Y$ or $y(x)$ , can be specialized to meet an "initial condition" by making an appropriate choice of $C[1]$ . Alternatively, mathematica can find the $C[i]$ if the initial condition(s) are included in DSolve. To do so, we replace ode as the first argument of DSolve by a list containing the ODE and its initial conditions. To solve our present ODE subject to the condition $y(1)=0$ , we write

DSolve[ {ode, y[1]==0}, y[x], x];

$\left\{\left\{y[x]\to \frac{1}{4}{e}^{-2x}\left(-e^2+e^{2x}-2e^{2x}x+2e^{2x}{x}^2\right)\right\}\right\}$

Possible initial conditions for an ODE in $y$ can include setting the value of $y$ at a specific point to infinity. The number of initial conditions should not exceed the order of the ODE; if this condition is violated, mathematica may generate meaningless results.

For ODEs of second or higher order, initial conditions can include specifications of the values of derivatives. When more than one initial condition is to be specified, they can be included as additional elements in the list that contains ode. For example,

DSolve[{ode, y[1]==4, y'[1]==0}, y[x], x];

Input of the following kind is also accepted:

DSolve[{ode, y[0]==y'[0]==0}, y[x], x];

Example 10.2.2 Solving an ODE with MATHEMATICA

Let's solve the ODE

$y^{″}+5y^{\prime }+4y=e^{-2x}\text{,}$

subject to the initial conditions $y(0)=0\text{,}{y}^{\prime }(0)=1$ . After we get the solution let's verify that the initial conditions are satisfied, plot the solution, and characterize the extremum shown on the plot. mathematica code for these operations can be as follows:

Set up the ODE:

ode = y''[x] + 5*y'[x] + 4*y[x] = Eˆ(-2*x);    (y" would be wrong!)

Now insert the initial conditions and solve the ODE:

DSolve[ {ode, y[0]==0, y'[0]==1}, y[x], x]

                $\left\{\left\{y[x]\to \frac{1}{6}{e}^{-4x}(-1-3e^{2x}+4e^{3x})\right\}\right\}$

We need to save the solution for further use.

Y = %[[1]][[1]][[2]]    $\frac{1}{6}{e}^{-4x}(-1-3e^{2x}+4e^{3x})$

Here the [[1]][[1]] removes the two sets of braces, and the [[2]] selects the material to the right of the symbol $\to$ .

We continue by computing $\mathit{dY}/\mathit{dx}$ , placing the result in YY:

YY = Expand[D[Y, x]]    $\frac{2e^{-4x}}{3}+e^{-2x}-\frac{2e^{-x}}{3}$

The command Expand is needed to force the collection of similar terms.

To check the initial conditions we now evaluate Y and YY at $x=0$ :

x = 0; {Y, YY}	{0, 1}
Clear[x]	Undefine x

The next command will work even if x remains defined. However, potential errors later on may be avoided by undefining quantities as soon as their definitions are no longer needed.

Assuming that we are interested in $y(x)$ for $x>0$ we plot it for $0<x<6$ :

Plot[Y, {x, 0, 6}]

We characterize the maximum in $y(x)$ by finding the $x$ value at which $\mathit{dy}/\mathit{dx}=0$ (this is a zero of YY). Thus,

Solve[YY == 0., x]

$\{\{\mathrm{x}\to -0.29568-1.77822\mathrm{i}\}\text{,}\{\mathrm{x}\to -0.29568+1.77822\mathrm{i}\}\text{,}\{x\to 0.59136\}\}$

The coding of this mathematica statement deserves two comments. First, because we are solving an algebraic (and not a differential) equation, the command is Solve, not DSolve. Second, because we want a decimal result and not a complicated analytical expression, we set YY equal to a decimal zero so that the computation will be reduced to decimal numbers.

The positive real root is the one we want. To isolate it we need to select the third item of the outer list, the first (and only item) of the list that is that third item, and then take the second argument of the construct $x\to 0.59136$ . This can be done all at once by the command

x = %[[3]][[1]][[2]]      0.59136

The final step is to evaluate $y(x)$ for this x value:

Y                  0.200176

We have therefore established that the maximum is at $x=0.591$ , with $y(x)=0.200$ .

Sometimes we will want our ODE solution established as a function, thereby permitting us to invoke it with an arbitrary argument. We proceed as follows:

Clear[x]; Clear[y];        Undefine, in case it is necessary

y[x_] ≔ Evaluate[Y];        y is now a function

To check what we have done,

y[t]              $\frac{1}{6}{e}^{-4t}(-1-3e^{2t}+4e^{3t})$

$\blacksquare$

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128010006000109

Dsolve Unable to Find Explicit Solution

Source: https://www.sciencedirect.com/topics/mathematics/dsolve

Share :