**3.9.7 Numerical Solution of Differential Equations**

The function NDSolve discussed in Section 1.6.4 allows you to find numerical solutions to differential equations. NDSolve handles both single differential equations, and sets of simultaneous differential equations. It can handle a wide range of ordinary differential equations as well as some partial differential equations. In a system of ordinary differential equations there can be any number of unknown functions , but all of these functions must depend on a single "independent variable" x

, which is the same for each function. Partial differential equations involve two or more independent variables.

Finding numerical solutions to ordinary differential equations.

NDSolve represents solutions for the functions as InterpolatingFunction objects. TheInterpolatingFunction objects provide approximations to the over the range of values xmin to xmax for the independent variable x.

NDSolve finds solutions iteratively. It starts at a particular value of x, then takes a sequence of steps, trying eventually to cover the whole range xmin to xmax.

In order to get started, NDSolve has to be given appropriate initial or boundary conditions for the and their derivatives. These conditions specify values for

[x], and perhaps derivatives

'[x], at particular points x. In general, at least for ordinary differential equations, the conditions you give can be at any x: NDSolve will automatically cover the range xmin to xmax.

This finds a solution for y with x in the range 0 to 2, using an initial condition for y[0].
In[1]:= **NDSolve[{y'[x] == y[x], y[0] == 1}, y, {x, 0, 2}]**

Out[1]=

This still finds a solution with x in the range 0 to 2, but now the initial condition is for y[3].
In[2]:= **NDSolve[{y'[x] == y[x], y[3] == 1}, y, {x, 0, 2}]**

Out[2]=

Here is a simple boundary value problem.
In[3]:= **NDSolve[{y''[x] + x y[x] == 0, y[0] == 1, y[1] == -1},**

y, {x, 0, 1}]

Out[3]=

When you use NDSolve, the initial or boundary conditions you give must be sufficient to determine the solutions for the completely. When you use DSolve to find symbolic solutions to differential equations, you can get away with specifying fewer initial conditions. The reason is that DSolve automatically inserts arbitrary constants C[

i] to represent degrees of freedom associated with initial conditions that you have not specified explicitly. Since NDSolve must give a numerical solution, it cannot represent these kinds of additional degrees of freedom. As a result, you must explicitly give all the initial or boundary conditions that are needed to determine the solution.

In a typical case, if you have differential equations with up to

derivatives, then you need to give initial conditions for up to

derivatives, or give boundary conditions at

points.

With a third-order equation, you need to give initial conditions for up to second derivatives.
In[4]:= **NDSolve[**

{ y'''[x] + 8 y''[x] + 17 y'[x] + 10 y[x] == 0,

y[0] == 6, y'[0] == -20, y''[0] == 84},

y, {x, 0, 1} ]

Out[4]=

This plots the solution obtained.
In[5]:= **Plot[Evaluate[ y[x] /. % ], {x, 0, 1}]**

With a third-order equation, you can also give boundary conditions at three points.
In[6]:= **NDSolve[**

{ y'''[x] + Sin[x] == 0,

y[0] == 4, y[1] == 7, y[2] == 0 }, y, {x, 0, 2}]

Out[6]=

Mathematica allows you to use any appropriate linear combination of function values and derivatives as boundary conditions.
In[7]:= **NDSolve[{ y''[x] + y[x] == 12 x,**

2 y[0] - y'[0] == -1, 2 y[1] + y'[1] == 9},

y, {x, 0, 1}]

Out[7]=

In most cases, all the initial conditions you give must involve the same value of x, say . As a result, you can avoid giving both xmin and xmax explicitly. If you specify your range of x as

x,

, then Mathematica will automatically generate a solution over the range to

.

This generates a solution over the range 0 to 2.
In[8]:= **NDSolve[{y'[x] == y[x], y[0] == 1}, y, {x, 2}]**

Out[8]=

You can give initial conditions as equations of any kind. In some cases, these equations may have multiple solutions. In such cases, NDSolve will correspondingly generate multiple solutions.

The initial conditions in this case lead to multiple solutions.
In[9]:= **NDSolve[{y'[x]^2 - y[x]^2 == 0, y[0]^2 == 4},**

y[x], {x, 1}]

Out[9]=

Here is a plot of all the solutions.
In[10]:= **Plot[Evaluate[ y[x] /. % ], {x, 0, 1}]**

You can use NDSolve to solve systems of coupled differential equations.

This finds a numerical solution to a pair of coupled equations.
In[11]:= **sol = NDSolve[**

{x'[t] == -y[t] - x[t]^2, y'[t] == 2 x[t] - y[t],

x[0] == y[0] == 1}, {x, y}, {t, 10}]

Out[11]=

This plots the solution for y from these equations.
In[12]:= **Plot[Evaluate[y[t] /. sol], {t, 0, 10}]**

This generates a parametric plot using both x and y.
In[13]:= **ParametricPlot[Evaluate[{x[t], y[t]} /. sol],**

{t, 0, 10}, PlotRange -> All]

Unknown functions in differential equations do not necessarily have to be represented by single symbols. If you have a large number of unknown functions, you will often find it more convenient, for example, to give the functions names like y[i].

This constructs a set of five coupled differential equations and initial conditions.
In[14]:= **eqns = Join[**

Table[ y[i]'[x] == y[i-1][x] - y[i][x], {i, 2, 4} ],

{y[1]'[x] == -y[1][x], y[5]'[x] == y[4][x],

y[1][0] == 1},

Table[ y[i][0] == 0, {i, 2, 5}]

]

Out[14]=

This solves the equations.
In[15]:= **NDSolve[eqns, Table[y[i], {i, 5}], {x, 10}]**

Out[15]=

Here is a plot of the solutions.
In[16]:= **Plot[ Evaluate[Table[y[i][x], {i, 5}] /. %],**

{x, 0, 10} ]

Options for NDSolve.

NDSolve allows you to specify the precision or accuracy of result you want. In general, NDSolve makes the steps it takes smaller and smaller until the solution it gets satisfies either the AccuracyGoalor the PrecisionGoal you give. The setting for AccuracyGoal effectively determines the absolute error to allow in the solution, while the setting for PrecisionGoal determines the relative error. If you need to track a solution whose value comes close to zero, then you will typically need to increase the setting for AccuracyGoal. By setting AccuracyGoal->Infinity, you tell NDSolve to use PrecisionGoal only.

NDSolve uses the setting you give for WorkingPrecision to determine the total number of digits to use in its internal computations. If you specify large values for AccuracyGoal or PrecisionGoal, then you typically need to give a somewhat larger value for WorkingPrecision. With the default setting of Automatic, AccuracyGoal and PrecisionGoal are both equal to the setting for WorkingPrecision minus 10 digits.

This generates a high-precision solution to a complex differential equation.
In[17]:= **NDSolve[{y'[x] == I/4 y[x], y[0] == 1}, y, {x, 1},**

AccuracyGoal -> 20, PrecisionGoal -> 20,

WorkingPrecision -> 25]

Out[17]=

Here is an approximation to

found from the solution.
In[18]:= **y[1] /. %**

Out[18]=

As mentioned above, NDSolve works by taking a sequence of steps in the independent variable x. NDSolve uses an adaptive procedure to determine the size of these steps. In general, if the solution appears to be varying rapidly in a particular region, then NDSolve will reduce the step size so as to be able to track the solution better.

This solves a differential equation in which the derivative has a discontinuity.
In[19]:= **NDSolve[**

{y'[x] == If[x < 0, 1/(x-1), 1/(x+1)],

y[-5] == 5},

y, {x, -5, 5}]

Out[19]=

NDSolve reduced the step size around

so as to reproduce the kink accurately.
In[20]:= **Plot[Evaluate[y[x] /. %], {x, -5, 5}]**

Through its adaptive procedure, NDSolve is able to solve "stiff" differential equations in which there are several components which vary with x at very different rates.

In these equations, y varies much more rapidly than z.
In[21]:= **sol = NDSolve[**

{y'[x] == -40 y[x], z'[x] == -z[x]/10,

y[0] == z[0] == 1},

{y, z}, {x, 0, 1}]

Out[21]=

NDSolve nevertheless tracks both components successfully.
In[22]:= **Plot[Evaluate[{y[x], z[x]} /. sol], {x, 0, 1},**

PlotRange -> All]

NDSolve follows the general procedure of reducing step size until it tracks solutions accurately. There is a problem, however, when the true solution has a singularity. In this case, NDSolve might go on reducing the step size forever, and never terminate. To avoid this problem, the option MaxSteps specifies the maximum number of steps that NDSolve will ever take in attempting to find a solution. For ordinary differential equations the default setting is MaxSteps->1000.

NDSolve stops after taking 1000 steps.
In[23]:= **NDSolve[{y'[x] == -1/x^2, y[-1] == -1}, y[x], {x, -1, 0}]**

-22

NDSolve::mxst: Maximum number of 1000 steps reached at the point x == -7.57086 10 .

Out[23]=

There is in fact a singularity in the solution at

.
In[24]:= **Plot[Evaluate[y[x] /. %], {x, -1, 0}]**

The default setting MaxSteps->1000 should be sufficient for most equations with smooth solutions. When solutions have a complicated structure, however, you may sometimes have to choose larger settings for MaxSteps. With the setting MaxSteps->Infinity there is no upper limit on the number of steps used.

To reproduce the full structure of the solution to the Lorenz equations, you need to give a larger setting for MaxSteps.
In[25]:= **NDSolve[ {x'[t] == -3 (x[t] - y[t]),**

y'[t] == -x[t] z[t] + 26.5 x[t] - y[t],

z'[t] == x[t] y[t] - z[t],

x[0] == z[0] == 0, y[0] == 1},

{x, y, z}, {t, 0, 20}, MaxSteps->3000 ]

Out[25]=

Here is a parametric plot of the solution in three dimensions.
In[26]:= **ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. %],**

{t, 0, 20}, PlotPoints -> 1000]

When NDSolve solves a particular set of differential equations, it always tries to choose a step size appropriate for those equations. In some cases, the very first step that NDSolve makes may be too large, and it may miss an important feature in the solution. To avoid this problem, you can explicitly set the option StartingStepSize to specify the size to use for the first step.

Finding numerical solutions to partial differential equations.

This finds a numerical solution to the wave equation. The result is a two-dimensional interpolation function.
In[27]:= **NDSolve[{D[y[x, t], t, t] == D[y[x, t], x, x],**

y[x, 0] == Exp[-x^2], Derivative[0,1][y][x, 0] == 0,

y[-5, t] == y[5, t]}, y, {x, -5, 5}, {t, 0, 5}]

Out[27]=

This generates a plot of the result.
In[28]:= **Plot3D[Evaluate[y[x, t] /. First[%]],**

{x, -5, 5}, {t, 0, 5}, PlotPoints->30]

This finds a numerical solution to the nonlinear sine-Gordon equation.
In[29]:= **NDSolve[{D[y[x, t], t, t] == D[y[x, t], x, x] + Sin[y[x, t]],**

y[x, 0] == Exp[-x^2], Derivative[0,1][y][x, 0] == 0,

y[-5, t] == y[5, t]}, y, {x, -5, 5}, {t, 0, 5}]

Out[29]=

Here is a plot of the result.
In[30]:= **Plot3D[Evaluate[y[x, t] /. First[%]],**

{x, -5, 5}, {t, 0, 5}, PlotPoints->30]