# Preface

This section is about application of standard Mathematica command NDSolve to the systems of ordinary differential equations.

Introduction to Linear Algebra with Mathematica

# Numerical solutions using NDSolve

The Wolfram Language function NDSolve is a general numerical differential equation solver. It can handle a wide range of ordinary differential equations (ODEs) as well as some partial differential equations (PDEs) and some differential-algebraic equations (DAEs). 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, which is the same for each function. NDSolve provides approximations to the solutions as InterpolatingFunction objects over the specified range of independent variable.

In order to get started, NDSolve has to be given appropriate initial or boundary conditions for the which the problem has a unique solution. These conditions specify values for the unknown function/solution, and perhaps its derivatives, if needed. A boundary value occurs when the unknown function is specified at multiple points. NDSolve can solve nearly all initial value problems that can symbolically be put in normal form (i.e. are solvable for the highest derivative order), but only linear boundary value problems.

Example: Consider the driven pendulum that is modeled by the following initial value problem:

$\ddot{x} (t) + \sin \left( x (t) \right) = 4\,\sin (2t) , \qquad x(0) =1, \quad \dot{x} (0) =0.$
This pendulum is set into motion with no initial velocity, denoted by $$\dot{x} ,$$ and it is subject by the periodic external force, 4 sin(2t). As usual, we denote the independent variable by t associated with time, and denote derivatives of x (which is angle of inclination with respect to downward position of a bob of unit mass) by dots. Here is the code to determine the approximate solution in the interval t ∈ [0, 100].
solution = NDSolve[{x''[t] + Sin[x[t]] == 4*Sin[2*t], x[0] == 1, x'[0] == 0}, x, {t, 0, 100}]
$$\left\{ \left\{ x \to \mbox{InterpolatingFunction}\left[ \begin{array}{c} \mbox{Domain:}\{\{0., \ 100.\}\} \\ \mbox{Output:\ scalar} \end{array} \right] \right\} \right\}$$
Next we plot separately the solution and its derivative/velocity and the using ParametricPlot we plot x vs its derivative, denoted by v.
Plot[Evaluate[{x[t], x'[t]} /. solution], {t, 0, 15}, PlotRange -> All, PlotStyle -> {{Thick, Blue}, {Thick, Red}}, PlotLabel -> "Solution (blue) and its derivative (red)" ]
ParametricPlot[{x[t], x'[t]} /. solution, {t, 0, 14}, AxesLabel -> {"x", "v"}, PlotLabel -> "Driven pendulum equation"]
■

As mentioned earlier, NDSolve works by taking a sequence of steps in the independent variable t. 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 to be able to better track the solution. 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 reached satisfies either the AccuracyGoal or 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. Generally, AccuracyGoal and PrecisionGoal are used to control the error local to a particular time step. For some differential equations, this error can accumulate, so it is possible that the precision or accuracy of the result at the end of the time interval may be much less than what you might expect from the settings of AccuracyGoal and PrecisionGoal.

NDSolve uses the setting you give for WorkingPrecision to determine the precision 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, both AccuracyGoal and PrecisionGoal are equal to half of the setting for WorkingPrecision. NDSolve uses the setting you give for WorkingPrecision to determine the precision to use in its internal computations.

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->10000. The default setting MaxSteps->10000 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.

Example: The pendulum equation contains sine function depending on the angle of inclination. If we approximate this sine function by its two term Maclaurin series $$\sin x \approx x - \frac{x^3}{6} + \cdots ,$$ we come to the so called the Duffing equation (undamped):

$\ddot{y} (t) + y - \frac{y^3}{6} = 4\,\sin (2t) , \qquad y(0) =1, \quad \dot{y} (0) =0.$
We expect that the solutions to the corresponding initial value problem of the above Duffing equation will be closed to our pendulum equation. However, it is not the case, as the following graphs show.
duffing = NDSolve[{y''[t] + y[t] - (1/6)*(y[t])^3 == 4*Sin[2*t], y[0] == 1, y'[0] == 0}, y, {t, 0, 100}]
... NDSolve: At t == 2.8049247505261308, step size is effectively zero; singularity or stiff system suspected.
So we use WorkingPrecision option:
duf100 = NDSolve[{y''[t] + y[t] - (1/6)*(y[t])^3 == 4*Sin[2*t], y[0] == 1, y'[0] == 0}, y, {t, 0, 100}, WorkingPrecision -> 30]
... NDSolve: At t == 2.804924871868654265896030557579754267854949671645843957253730., step size is effectively zero; singularity or stiff system suspected.
When we plot two solutions, for the pendulum equation and for Duffing equations, we see that the coincide only within short time interval
Plot[Evaluate[{x[t] /. solution , y[t] /. duffing} ], {t, 0, 2.4}, PlotRange -> All, PlotStyle -> {{Thick, Blue}, {Thick, Red}}, PlotLabel -> "Pendulum (blue) and Duffing (red)" ]
When we increase the interval, we come to a similar problem:
sol150 = NDSolve[{x''[t] + Sin[x[t]] == 4*Sin[2*t], x[0] == 1, x'[0] == 0}, x, {t, 0, 150}, WorkingPrecision -> 30]
... NDSolve: Maximum number of 10000 steps reached at the point t == 102.1744813327041831765056200018136296715330..
This message actually indicates that the solution of this driven pendulum equation is unbounded. We try to reduce the amplitude of the oscillating force:
sol150 = NDSolve[{x''[t] + Sin[x[t]] == 1*Sin[2*t], x[0] == 1, x'[0] == 0}, x, {t, 0, 150}, WorkingPrecision -> 30]
ParametricPlot[{x[t], x'[t]} /. sol150, {t, 0, 100}, AxesLabel -> {"x", "v"}, PlotLabel -> "Driven pendulum equation"]
Plot[Evaluate[{x[t], x'[t]} /. sol150], {t, 0, 15}, PlotRange -> All, PlotStyle -> {{Thick, Blue}, {Thick, Red}}, PlotLabel -> "Solution (blue) and its derivative (red)" ]
Then we repeat calculations with the same driven force for Duffing equation.
duf150 = NDSolve[{y''[t] + y[t] - (1/6)*(y[t])^3 == 1*Sin[2*t], y[0] == 1, y'[0] == 0}, y, {t, 0, 150}, WorkingPrecision -> 30] Plot[Evaluate[{x[t] /. sol150 , y[t] /. duf150} ], {t, 0, 24}, PlotRange -> All, PlotStyle -> {{Thick, Blue}, {Thick, Red}}, PlotLabel -> "Pendulum (blue) and Duffing (red)" ]
You can also specify the solution method:
duf150 = NDSolve[{y''[t] + y[t] - (1/6)*(y[t])^3 == 1*Sin[2*t], y[0] == 1, y'[0] == 0}, y, {t, 0, 150}, WorkingPrecision -> 30, Method -> "ExplicitRungeKutta"]
However, upon plotting, we will not recognize the difference.
As we see, for such amplitude (equals 1), both equations have bounded solutions that we very closed.    ■

NDSolve has several different methods built in for computing solutions as well as a mechanism for adding additional methods. With the default setting Method->Automatic, NDSolve will choose a method which should be appropriate for the differential equations. For example, if the equations have stiffness, implicit methods will be used as needed, or if the equations make a DAE, a special DAE method will be used. In general, it is not possible to determine the nature of solutions to differential equations without actually solving them: thus, the default Automatic methods are good for solving as wide variety of problems, but the one chosen may not be the best one available for your particular problem. Also, you may want to choose methods, such as symplectic integrators, which preserve certain properties of the solution. \right) \)
$$y_{n+1} = y_n + h_n \,f\left( x_{n+1/2}, y_{n+1/2} \right)$$
Order Method Formula
1 Explicit Euler $$y_{n+1} = y_n + h_n f\left( x_n, y_n \right)$$
2 Explicit Midpoint $$y_{n+1/2} = y_n + h_n \frac{1}{2} \,f\left( x_n, y_n \right)$$
1 Backward or implicit Euler $$y_{n+1} = y_n + h_n f\left( x_{n+1}, y_{n+1} \right)$$
2 Implicit Midpoint $$y_{n+1} = y_n + h_n f \left( x_{n+1/2} , \frac{1}{2} \left( y_n + y_{n+1} \right) \right)$$
2 Trapezoidal $$y_{n+1} = y_n + h_n \frac{1}{2}\, f (x_n , y_n ) + h_n \frac{1}{2}\, f \left( x_{n+1} , y_{n+1} \right)$$
1 Linearly Implicit Euler $$\left( I - h_n J \right) \left( y_{n+1} - y_n \right) = h_n f \left( x_n , y_n ) \right)$$
2 Linearly Implicit Midpoint $$\left( I - \frac{h_n}{2}\, J \right) \left( y_{n+1/2} - y_n \right) = \frac{h_n}{2} \, f\left( x_n , y_n \right)$$
$$\left( I - \frac{h_n}{2}\, J \right) \frac{\Delta y_n - \Delta y_{n-1/2}}{2} = \frac{h_n}{2} \,f\left( x_{n+1/2}, y_{n+1/2} \right) - \Delta y_{n-1/2}$$
Here Δyn = yn+1 - yn, I denotes the ideentity matrix, and J denoted the Jacobian matrix.

When NDSolve computes a solution, there are typically three phases. First, the equations are processed, usually into a function that represents the right-hand side of the equations in normal form. Next, the function is used to iterate the solution from the initial conditions. Finally, data saved during the iteration procedure is processed into one or more InterpolatingFunction objects. Using functions in the NDSolve context, you can run these steps separately and, more importantly, have more control over the iteration process. The steps are tied by an NDSolve`StateData object, which keeps all of the data necessary for solving the differential equations.