NDSolve[.]
Accurate Numerical Solutions to Differential Equations

by
K. D.  Stroyan
University of Iowa

copyright 1997 by Academic Press, Inc. - All rights reserved.

Special Mathematica Functions used in this NoteBook

eulerApprox[f_,{t0_,x0_},dt_,final_] :=
Module[{steps,G},
steps = Floor[final/dt];
G[{tt_,xx_}] := {tt+dt,xx + f[xx] dt};
NestList[ G , {t0,x0} , steps ]]

Notebook Overview

This NoteBook shows how to use Mathematica's built-in accurate numerical differential equation solver  NDSolve[.].  Throughout the course we encounter the extremely important general problem: Given a starting value and a rule for how a quantity changes, find the quantity at later times.  This is called an "Initial Value Problem."

In 1 dimension a continuous initial value problem is:

GIVEN:
(1)  the start:  y[(Subscript[t, 0])] = (Subscript[y, 0])

(2)  the change:  (dy/dt) = f[t,y]
in terms of a change law function f[t,y]

FIND:
The function y[t]  for  t > (Subscript[t, 0]).

A 2 dimensional continuous initial value problem:

GIVEN:
(1) the start
	x[(Subscript[t, 0])] = (Subscript[x, 0])
	y[(Subscript[t, 0])] = (Subscript[y, 0])

(2)	
	(dx/dt) = f[t,x,y]
	(dy/dt)= g[t,x,y]

in terms of change laws f[t,x,y] and g[t,x,y]

FIND:
The vector function  X[t] = {x[t],y[t]}  for  t >(Subscript[t, 0]).

Sometimes we can find exact symbolic solutions of such problems, but often we need to approximate the solution numerically.  (Many such problems can not be solved symbolically.)

Euler's Approximation

Euler's method of approximation is simple conceptually; we just use the increment approximation recursively.  The Increment Approximation (or Microscope Equation) is

x[t + dt] = x[t] + x'[t] dt + ϵ dt,   with  ϵ  ≈  0  when  dt  ≈  0.

Since we know  f[t,x] and x[(Subscript[t, 0])], we can find x[(Subscript[t, 0])+dt] as follows:

x[(Subscript[t, 0]) + dt] ≈ x[(Subscript[t, 0])]+ f[(Subscript[t, 0]), x[(Subscript[t, 0])]] dt

since  x'[(Subscript[t, 0])] = f[(Subscript[t, 0]),x[(Subscript[t, 0])]].  Then we can find  x[(Subscript[t, 0])+2dt]  using  x[(Subscript[t, 0])+dt]:

x[(Subscript[t, 0]) + 2dt]  ≈ x[(Subscript[t, 0]) + dt]+ f[(Subscript[t, 0]) + dt, x[(Subscript[t, 0]) + dt]] dt

Continuing recursively,

x[(Subscript[t, 0]) + (n+1)dt]  ≈  x[(Subscript[t, 0]) + n dt]+ f[(Subscript[t, 0]) + n dt, x[(Subscript[t, 0]) + n dt]] dt
		
This gives us approximate vlaues for

x[(Subscript[t, 0])], x[(Subscript[t, 0])+dt], x[(Subscript[t, 0])+2dt], x[(Subscript[t, 0])+3dt], x[(Subscript[t, 0])+4dt], ...

Unfortunately, Euler's approximation is not very acurate.  Mathematica has a built-in solver  NDSolve[.] that gives a more accurate approximation.  (Roughly speaking it gains accuracy by taking the bend in the graph into account, rather than just the slope of the increment approximation.)

We studied Euler's method in more detail earlier in the course, see the NoteBook EulerApprox.

References to the Text

This NoteBook goes with Chapters 21and 22 of the text.

Step-by-step:Comparison of Euler and NDSolve[.]

The following specific instance has the simplest possible change rule in the dependent variable:

(1)	y[0] = 1
and
(2)	(dy/dt) = y

This initial value problem (or differential equation) has the analytical solution

y[t] = (e^t) = Exp[t]

The next program solves the problem numerically using Euler's approximation and NDSolve[.].  It then compares the numerical solutions to the exact one.

Clear[y, t, f, soln]; f[y_] := y; t0 = 0.0; y0 = 1.0; finalTime = 1.0; dt = 0.1; eulerPlot = ListPlot[eulerApprox[f, {t0, y0}, dt, finalTime]]; iVprob={y[t0] == y0,y'[t]==f[y[t]]};soln = NDSolve[iVprob, {y[t]}, {t, 0, finalTime}]; accPlot = Plot[y[t] /. soln, {t, 0, finalTime}]; exactPlot = Plot[Exp[t], {t, 0, finalTime}]; Show[accPlot, eulerPlot, exactPlot]; Table[{t, y[t], N[Exp[t]]} /. soln, {t, 0, finalTime, 0.1}]

The computer NDSolve[.] solution "soln" to the above program is a list of "Interpolating Functions."  This means that Mathematica solves the differential equation at certain times and can interpolate to other times.  The program above shows how to extract graphs and values from this using "/.".

Mathematica's Description of NDSolve[.]

Enter the following command to see Mathematica's explaination of NDSolve[.]

??NDSolve

Practice Exercises

1)  Modify the program above to solve the initial value problem:
	y[0] = 1
	dy = y dt,  for 0 < t < 3

2)  Modify the program above to solve the initial value problem:
	y[0] = 2
	dy = 5 y dt,  for 0 < t < 7	
In this case, you have to change the exact solution as well as the differential equation.

Combined Computations

Sample  NDSolve[.] Solution to  (dx/dt) = f[x, t]

This example shows the accurate solution to a one equation or a one dimensional initial vlaue problem.

Clear[t, x, f]; f[x_, t_] := x (1 - (x/10) - (t/5))equn=x'[t]==f[x[t], t];tInitial = 0; xInitial = 4; initCond = x[tInitial] == xInitial; iVprob = {equn, initCond}tFinal = 10; sol = NDSolve[iVprob, x[t], {t, tInitial, tFinal}]

Plot[x[t] /. sol, {t, tInitial, tFinal}, PlotRange -> {0, 10}];

Practice Exercises

Modify the program above to solve the initial value problem:
	x[0] = 1
	dx = x (1 - (x/3)) dt,  for 0 < t < 15

Sample  NDSolve[.] Solution to
	dx/dt = f[x,y]
	dy/dt = g[x,y]

This example shows the accurate solution to the two equation or two dimensional differential equation.

Clear[t, x, y, f, g]; f[x_, y_] := x (1 - (x/10) - (y/5)); g[x_, y_] := y (1 - (x/7) - (y/7)); tInitial = 0; xInitial = 4; yInitial = 6; iVprob={x[tInitial] == xInitial,y[tInitial] == yInitial,x'[t]==f[x[t], y[t]],y'[t]==g[x[t], y[t]]};tFinal = 10; sol = NDSolve[iVprob, {x[t], y[t]}, {t, tInitial, tFinal}]

\(Plot[x[t] /. sol, {t, tInitial, tFinal}, PlotRange -> {0, 10}, 
      AxesLabel -> {"t", "x"}]; 

    Plot[y[t] /. sol, {t, tInitial, tFinal}, PlotRange -> {0, 10}, 
      AxesLabel -> {"t", "y"}]; 

    ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, tInitial, tFinal}, 
      AspectRatio -> 1, PlotRange -> {{0, 10}, {0, 10}}, 
      AxesLabel -> {"x", "y"}]; \)

Practice Exercises

1)  Modify the example above to solve these equations with the initial condition
	x[0] = 5
	y[0] = 5

2)  Modify the example above to solve the equations
	dx = x(1 - x/7 - y/7) dt
	dy = y(1 -x/10 -y/5) dt
			
for the initial conditions
	
	x[0] = 4		
	y[0] = 6

for  0 < t < 20

Exercises

Exercise 21.2.2

Solve the following differential equations using a modification of the NDSolve[.] program above and plot the solution.

x[0] = 2;
(dx/dt) =  3 x;

x[0] = 2;
(dx/dt) = 3 x (1 - (x/5));

x[0] = 7;
(dx/dt) = 3 x (1 - (x/5));

x[0] = 5;
(dx/dt) = 3 x (1 - (x/5));

x[0] = 5;
(dx/dt) = 3 x (1 - x Cos[t]);

y[0] = 2;
(dy/dx) =  3 x;

y[0] = 2;
(dy/dx) =  3 y;

y[0] = 2;
(dy/dx) =  3 x y;

Exercise 21.4.3

Use this NoteBook to solve one of the differential equations in Exercise 21.35 where Euler's method seemed rather inaccurate.  This can be done by modifying the comparison of Euler's method and NDSolve[.] for the sample of Exp[t] above.