IMSE Chapter 6

IMSE-06(NDSolve).nb
Numerical Solution of Ordinary Differential Equations
RCH 8/6/01
revised 6/11/02 - rjerwin

This notebook illustrates how Mathematica solves differential equations numerically.  For linear differential equations, this is sometimes a useful alternative to analytic solutions.  For nonlinear differential equations, it is a necessity since nonlinear differential equations usually do not have solutions that can be expressed in analytic form.

The differential equation for a simple harmonic oscillator is set up in the next Input cell using the command NDSolve[] (numerical differential equation solver).  There are several things to note about the format.  First, the differential equation is specified with the double equals (==) form.  In addition to specifying the differential equation, you must specify the appropriate initial conditions.  For a second-order differential equation, you need to specify two initial conditions.  Here we specify the initial value of the first-derivative (in physics terms, the initial velocity) and the initial value of the dependent function itself (in physics terms, the initial position).  The initial conditions also make use of the == symbol.  The equation and the initial conditions form a list (and hence are enclosed in curly brackets).  We then specify the dependent variable, here y.  Finally comes a list (more curly brackets) that specifies the independent variable (here x) and the range of values for that variable.

Activate the following cell.

solution  =  
      NDSolve[{y''[x] + y[x] == 0, 
          y[0] == 0, y'[0] == 1}, y, {x, 0, 12}]

The response from Mathematica is a message indicating that the differential equation has been solved numerically and is represented by what Mathematica calls an "Interpolating Function."  This means that Mathematica has developed rules for giving the dependent variable as a function of the independent variable over the range specified.  The interpolating function is in the form of a list assigned to the name solution. 

To see what is in the list solution, activate the following three cells.

solution[[1]]

This level indicates that solution is a list.  Note the curly brackets.

solution[[1, 1]]

This level indicates that the solution for y is the first sublist.  (Later we shall see that we can solve two or more differential equations simultaneously.  Each solution will have its own sublist.)

solution[[1, 1, 1]]

The first element of the third level sublist is the name associated with the Interpolating Function.

solution[[1, 1, 2]]

Finally, we get down to the Interpolating Function itself.  We see that the actual Interpolating Function is located at the second element of the third sublist level "down" into the list.

If you want to plot the solution or use it in further calculations, you need to evaluate the solution or assign it to a function.  There are several ways to do this.  One way is to assign the relevant part of the list to a function, which, with characteristic lack of imagination, we call f. 

Activate the following cell.

f = solution[[1, 1, 2]]

When we now use f[x] in expressions, Mathematica uses the Interpolating Function assigned to f[x] to evaluate the function.  For example, in a Plot[] command, Mathematica will use the Interpolating Function to calculate the appropriate numerical values for the graph.

Plot[f[x], {x, 0, 12}, 
        AxesLabel -> {"x", "y[x]"}];

An Aside:  You can, with consumer warnings attached, use the Interpolating Function outside the range of definition.  Usually the result is nonsense.  Activate the following cell to see what happens.

Plot[f[x], {x, 0, 15}, 
        AxesLabel -> {"x", "y[x]"}];

Yet another way of doing the plot is to use the Evaluate[] command.  This tells Mathematica to evaluate the Interpolating Function at the points required to generate the graph.  The /. symbol combination (called the "replacement operator")  tells Mathematica to use the interpolating function associated with solution to evaluate the function y[x].  The replacement operator is also used in various symbolic transformations in Mathematica.

Plot[Evaluate[y[x] /. solution], {x, 0, 12}, 
        AxesLabel -> {"x", "y[x]"}];

We can use the solution in the form of f[x] for further computations.  For example, we might want to find the area between the curve and the x axis in the range 0 to π.

NIntegrate[f[x], {x, 0, π}]

Is the result consistent with what you expect?

NOTE:  If you try to do the numerical integration with y[x], Mathematica gives you an error message because you still have not assigned the interpolating function to the name y[x].  (The Evaluate[] command simply worked out the numerical values without making any name assignment.)   Activate the next cell to see what happens.

NIntegrate[y[x], {x, 0, π}]

The next cell assigns the solution to y and then does the numerical integration.

y = solution[[1, 1, 2]]
NIntegrate[y[x], {x, 0, π}]

Coupled Differential Equations

We can also find numerical solutions for coupled differential equations.  Here is a simple example with two dependent variables y and z and the independent variable t.

Two NOTES:

1.  If we want to use y[t] in the definition of the differential equations, we must use Clear[y] because we have defined y previously.
2.  Pay attention to how the assignment statements are used to associate the functions y1 and z1 with the list sol.  The Interpolating Function for y1 is associated with the first sublist while the Interpolating Function for z1 is associated with the second sublist.

Clear[
        y];  (*y was used above .   Need to Clear it.*) 
\(sol = 
      NDSolve[{\(y'\)[t] - z[t] == 0, \(z'\)[t] + y[t] == 0, 
          y[0] == 0, z[0] == 1}, {y, z}, {t, 0, 
          12}]

 (*Assign\ the\ \(\(solutions\)\(.\)\)*) \)
y1 = sol[[1, 1, 2]];
z1 = sol[[1, 2, 2]];
Plot[{y1[t], z1[t]}, {t, 0, 12}, 
        AxesLabel -> {"time", "y and z"}];

A cleaner way to set up the solution of multiple equations is to put the equations and the initial conditions into a list.  This makes the reading and editing much easier.  Let's repackage the previous example.

Clear[y, 
      z]  (*Clear y and z .   y was used above.*) 
 (*First construct the list of equations and initial conditions. *) 
\(\(eqList = {\(y'\)[t] - z[t] == 0, \(z'\)[t] + y[t] == 0, 
          y[0] == 0, z[0] == 1};\)
 (*\ 
      Now\ solve\ the\ equations\ \(\(numerically\)\(.\)\)\ *) \)
\(sol = 
      NDSolve[eqList, {y, z}, {t, 0, 
          12}]

 (*Assign\ the\ \(\(solutions\)\(.\)\)*) \)
y = sol[[1, 1, 2]];
z = sol[[1, 2, 2]];
Plot[{y[t], z[t]}, {t, 0, 12}, 
        AxesLabel -> {"time", "y and z"}];

Numerical Accuracy

How do you know that the solutions you obtain numerically are accurate?  This is a subtle and important issue.  If you have some idea of what the solution should look like, as we did in the examples given above, then you can make a reasonable evaluation of the results.  If you don't know in advance what the solution should be, then there are several strategies for judging accuracy:

1.  With your solutions, calculate some quantity, such as the total energy, which is expected to remain constant.  Since the numerical solutions are not perfectly accurate, you will expect the energy to change by some amount.  By setting some limits on the acceptable energy variations, you can judge the accuracy of the results.  However, to use this method you do need some "conserved" quantity and that may not always be available.
2.  Change the numerical procedures and see how the the solutions change.  If the solutions do not change (within some limits you deem acceptable), then you have some justification for assuming that the solutions are accurate to within that limit.  However, you still have no guarantee that they are accurate.

Mathematica's NDSolve[] uses an adaptive stepsize procedures, that is, it changes the integration step size, with large step size used for regions in which the solution is changing slowly and small step sizes for those regions in which the solution is varying rapidly.  You can tell Mathematica when to stop adjusting the step size by invoking various options for NDSolve[].  For example, you can change the integration step size and/or the accuracy that Mathematica imposes on the solution.  Mathematica also can switch procedures for differential equations that have rapidly varying parts.  The details of these methods are described in S. Wolfram's book Mathematica (section 3.9.7, Numerical Solution of Differential Equations).  Here is a brief summary of the options:

WorkingPrecision->n sets the number of digits n used in internal calculations.  It can be made as large as you want, but the execution time increases rapidly.  The default value is $MachinePrecision, typically 16 digits.  As mentioned before, you can find the value of $MachinePrecision for your computer by activating the following command.

$MachinePrecision

AccuracyGoal->n sets the desired accuracy n of the solution (basically the goal for the change in solution with decreasing step size) given in absolute terms -- for example, 0.0001.  The default value is the value of WorkingPrecision minus 10 digits.  By setting AccuracyGoal->Infinity you force Mathematica to use the PrecisionGoal.

PrecisionGoal->n sets the desired accuracy n of the solution relative to the actual numerical value of the solution.  The default value is WorkingPrecision minus 10 digits.

StartingStepSize->n sets the initial step size n that Mathematica uses for its integration procedures.  Sometimes the default value can be too large, particularly for solutions that change rapidly, and the procedure may appear to have high accuracy when it actually doesn't.

MaxSteps->n sets the maximum number of integration steps n that Mathematica will use to integrate the differential equation.  Rapidly varying functions require smaller (and hence more) steps to achieve a specified accuracy in the solution.  The default setting is 1000.

We will now illustrate these options by using our first example of the simple harmonic oscillator.  We have added explicit statements about WorkingPrecision and the maximum number of integration steps MaxSteps.  With the values given below, the integration proceeds well until about x = 200, but then blows up as Mathematica exhausts the default number of steps.

Clear[y]

xmax  = 
      250;    
solution2  =  
      NDSolve[{y''[x] + y[x] == 0, 
          y[0] == 0, y'[0] == 1}, y, {x, 0, xmax}, 
        AccuracyGoal -> Automatic, 
        MaxSteps -> 1000];
g = solution2[[1, 1, 2]];
Plot[g[x], {x, 0, xmax}];

Change the maximum number of steps to 5000 and reactivate the cell.  Everything should run smoothly now.

Change the maximum number of steps back to 1000 and decrease the AccuracyGoal from Automatic to 4.  You should see that the solution doesn't blow up, but the amplitude of the sine wave is slowly changing.

To make the effect more dramatic, change the AccuracyGoal to 3.  Now the sine wave amplitude seems to be increasing -- a clear sign that the solution is not accurate.

The lesson here is that the parameters of the numerical routines can affect your results, sometimes in unexpected ways. As a general rule, you should always check your numerical results by modifying those parameters to see how the results change.