Numerical Solving of Differential Equations
in Mathematica

In this talk, I plan to show you a brief overview of the possibilities available for solving differential equations in Mathematica.  I will start by showing a few examples with NDSolve.  Next I will use the Schrodinger equation as an example to demonstrate in detail how the error estimates used by NDSolve work and how you can use Mathematica's features to write a specialized code more efficient than NDSolve.  Finally, as time permits, I will show you some features planned for NDSolve in future versions.

You should be able to find the notebook and packages for this talk at http://wri.members.com/rknapp/DeveloperConference99

Here is an example which shows several of the basic features (and limitations) of NDSolve:

ndsol  =  
      NDSolve[{x'[t]^2  ==  x[t],  x[0]  ==  1}, 
        x, {t, 0, 100}]

First, note that the solution is returned in terms of an approximate function based on interpolation of the solution data.  The first argument gives the domain.

end  =  ndsol[[1, 1, 2, 1, 1, 2]]

The important point about this is that within the given domain, the function can be treated like any other function which evaluates numerically. For example, they can be plotted:

Plot[Evaluate[Abs[x[t]]  /.  ndsol], {t, 0, end}];

or differentiated:

Plot[
        Evaluate[Log[10, Abs[x'[t]^2  -  x[t]]]  /.  ndsol], {t, 0, 
          end},  PlotStyle -> {RGBColor[1, 0, 0],  
            RGBColor[0, 0, 1]}];

This is the residual function.  You should really look at this at the data point produced by the method

times  =  (x  /.  ndsol)[[All, 3, 1]]; 
    Show[Block[{$DisplayFunction  =  
            Identity}, {
ListPlot[
            Map[Function[
                Evaluate[{#,  
                    Log[10, 
                      Abs[x'[#]^2  -  x[#]  /.  
                          ndsol[[1]]]]}]],  times[[1]]],  
            PlotStyle -> {RGBColor[1, 0, 0],  PointSize[0.02]}], 
          ListPlot[
            Map[Function[
                Evaluate[{#,  
                    Log[10, 
                      Abs[x'[#]^2  -  x[#]  /.  
                          ndsol[[2]]]]}]],  times[[2]]],  
            PlotStyle -> {RGBColor[0, 0, 1],  
                PointSize[0.02]}]}],  %];

While the residual is realated to the error, you should not mistake it for a direct measure of the error.  An solution with small residual can have a large error.

For this simple example, the exact solution can be found with DSolve:

dsol  =  
      DSolve[{x'[t]^2  ==  x[t],  x[0]  ==  1}, x, 
        t]

The second solution in this case does have a small error:

Plot[
        Evaluate[(x[t]   /.  ndsol[[2]])  -  (x[t]  /.  
                dsol[[2]])], {t, 0, 100}];

The larger residual for the first indicates the possibility of a problem, which indeed there is:

Plot[
        Evaluate[
          Abs[(x[t]   /.  ndsol[[1]])  -  (x[t]  /.  
                  dsol[[1]])]], {t, 0, end}];

The exact solution is tangent to the axis here.

Plot[(x[t]  /.  dsol[[1]]), {t, 0, 4}];

The problem is that NDSolve has to work from the normal form x'[t] == -Sqrt[x[t]], and at x == 0, the function Sqrt[x] does not have Lipshitz continuity, which is the condition necessary to prove the basic existence-uniqueness theorm for initial vaue problems.

A stiff equation is one for which numerical stability requirements force some numerical methods to sue extremely small timesteps. Typically, a stiff problem has components of the solution varying on quite disparate time scales.

A good example of a stiff problem is the Van der Pol oscillator, which has a limit cycle

ϵ  =  50;

vsol  =  
      First[NDSolve[{x'[t]  ==  
              y[t], y'[
                t]  ==  ϵ (1  -  x[t]^2) y[t]  -  
                x[t], x[0]  ==  2,  y[0]  ==  0}, {x, 
            y}, {t, 0, 100}]]

A look at the graph of one component indiccates that the solution may be periodic.

Plot[x[t]  /.  vsol, {t, 0, 100}];

A phase plane view shows the limit cycle

ParametricPlot[{x[t], y[t]}  /.  vsol, {t, 0, 100},  
        PlotRange -> All,  PlotPoints -> 1000];

The RungeKutta method does not handle stiff equations and can't compute the solution nearly as well:

rvsol  =  
      First[NDSolve[{x'[t]  ==  
              y[t], y'[
                t]  ==  ϵ (1  -  x[t]^2) y[t]  -  
                x[t], x[0]  ==  2,  y[0]  ==  0}, {x, 
            y}, {t, 0, 100},  Method -> RungeKutta]]

Find (approximately) the period

half  =  
      y  /.  First[
          NDSolve[{x'[t]  ==  
                y[t], y'[
                  t]  ==  ϵ (1  -  x[t]^2) y[
                      t]  -  x[t], x[0]  ==  2,  
              y[0]  ==  0}, {x, y}, {t, 0, 100},  
            StoppingTest -> {y[t]  >  0}]]

end  =  half[[3, 1, -1]]; 
    Plot[half[t], {t, 0,  41.0302}];

FindRoot[half[t], {t, end,  half[[3, 1, -2]]}, 
      Compiled -> False]

This example comes from a notebook sent to me by Kevin McCann of John Hopkins University.  It provides a nice example not only to show how NDSolve can work with PDE's, but also to give you an idea of some of what NDSolve is doing to make error estimates.

Set up

Parameters

ℏ = 1.0;
m = 1.0;
σ = 5.0;
τ = 2  (σ^2);
x0 = -50.0;
k0 = 1.0;
L = 200.0;
T = 100.0;

Initial condition

Ψ[x_, 
        t_] := (1  -  
            UnitStep[
              Abs[x]  -  
                100*(1  -  $MachineEpsilon)]) ((((1/(2  π (σ^2)))))^(1/
              4)) (1/(Sqrt[1  +  I (t/τ)])) (E^(I (k0 (x - 
                        x0) - ((ℏ (k0^2))/(2  m)) 
                    t))) (E^(-((((x - (x0 + ((ℏ k0)/m) t)))^2)/(4  (σ^2) (1 + I t/τ)))))

Potential

V[x_]  =  100.  UnitStep[x];

Solution with NDSolve

sol  =  
      First[NDSolve[{D[u[x, t], 
                t]  ==  ((I ℏ)/(2  m)) 
                  D[u[x, t], x, x]  + (I/ℏ) 
                  V[x] u[x, t],  
            u[x, 0]  ==  Ψ[x, 0],  
            u[-L/2, t]  ==  0,  u[L/2, t]  ==  0}, 
          u, {x, -L/2, L/2}, {t, 0, T}]]

Plot3D[Abs[u[x, t]  /.  sol], {x, -100, 100}, {t, 0, 100},  
        PlotPoints -> 50,  PlotRange -> All];

Ouch!  This does not look so good.

The problem here is that NDSolve makes an error estimate based on a few spatial points, then increases the number of points until the Accuracy and Precision goal targets (default 4) are met.  Sometimes, the default number of points misses features of the solution, so the error estimate is initially wrong.  You can take care of this by requesting more points:

Timing[
      sol  =  First[
          NDSolve[{D[u[t, x], 
                  t]  ==  ((I ℏ)/(2  m)) 
                    D[u[t, x], x, x]  + (I/ℏ) 
                    V[x] u[t, x],  
              u[0, x]  ==  Ψ[x, 0],  
              u[t, -100]  ==  u[t, 100]  ==  0}, 
            u, {t, 0, 100}, {x, -100, 100},  
            MaxStepSize -> 2]]]

Now we get the expected reflection...

Plot3D[Abs[u[t, x]  /.  sol], {t, 0, 100}, {x, -100, 100},  
        PlotPoints -> 50,  PlotRange -> All];

So how good a solution is this?  It is even harder to guage the quality of the solution of a PDE than one of an ODE since more can go wrong and the behavior can be much more complicated.  One thing which is often quite useful is conserved quantities if there are any. The Schrodinger equation conserves Abs[u]^, or total probability.

ndsol  =  u  /.  sol;

Timing[
        ndc  =  
          ListPlot[
            Map[Function[{#, 
                  NIntegrate[Abs[ndsol[#, x]]^2,  
                    Evaluate[Join[{x},  ndsol[[3, 2]]]]]}], 
              Range[100]]]][[1]]

Can we do better?

Here is the notebook sent to me by Kevin McCann.  He had written a code to solve the equation using a Crank-Nicholson scheme, but was hoping it could do this problem in fewer than the 4-5 hundred seconds his code was taking.  In my previous talk on "Generally Efficient Numerical Programming in Mathematica " earlier today, I showed how this scheme could be speeded up tremendously.

Crank-Nicolson is an implicit second order scheme for 1 + n dimensional PDE's:

(Subscript[u, t]) ==  ℱ (u)

which discretizes symmetrically in time

((Subscript[u, i])^(n + 1)) -  
        ((Subscript[u, i])^n)  == (( Q (((Subscript[u, i])^(n + 1)))  +  Q (((Subscript[u, i])^n)))/2)

where Q is typically a finite difference operator:  The discretization used here is the standard centered finite difference operator:

Q ((Subscript[u, i]))  =  (((Subscript[u, i - 1]) -  2 (Subscript[u, i]) + (Subscript[u, i + 1]))/Δ)(x^2)

The differencing is implicit since ((Subscript[u, i])^(n + 1))appears on both sides of the equation.  For a linear equation like the Schordinger equation, this requires the solution of a system of linear equations (tridiagonal) at each step.  For a nonlinear equation, Newton-type iterations may be necessary.

Discretization Parameters

These are the parameters he was using:

Nx = 1001;
Nt = 401;

Δx  =  
        (L/(Nx - 1));

Δt  =  
        (T/(Nt - 1));

ρ  =  ((I ℏ Δt)/(2  m (((Δx))^2)));

grid   =  Δx Range[0, Nx - 1]  -  
          (L/2);

init  =  
        Developer`ToPackedArray[Ψ[grid, 0]];

vlist  =  
        Take[Developer`ToPackedArray[(I/ℏ) Δt Map[V, grid],  Complex], {2, -2}];

This sets up a sparse matrix  in the form {row, col}->value, with any matrix position not explicitly given taken to be 0.  In my previous talk today, I showed that this was the most efficient way of doing the matrix solution.

spam  =  {{
Join[Range[Nx - 2],  
              Range[Nx - 3],  Range[2, Nx - 2]], 
Join[
              Range[Nx - 2],  Range[2, Nx - 2],  
              Range[Nx - 3]]}, 
Developer`ToPackedArray[
            Join[2 + 2 ρ  +  vlist,  Table[-ρ, {Nx - 3}], 
              Table[-ρ, {Nx - 3}]]]};

This produces an LU decomposition of the sparse matrix...

time[tlu  =  Experimental`ExtendedLinearSolve[spam]]

... in a functional form which can then backsolve very quickly

(Note time is an adaptation I have made of the Timing command for handling timings which
vary over many different orders of magnitude.)

time[tlu[ Range[Nx - 2]]]

Crank

p is the finite difference stencil

onezero  =  
        Developer`ToPackedArray[{0.   +  
              0.  I}];
p  =  
        Developer`ToPackedArray[{ ρ,  2  -  2 ρ,  ρ},  
          Complex];

This takes one time step of the Crank-Nicholson method ...

crank  =  
        Function[
          PadLeft[ 
            tlu[ ListCorrelate[p, #]  -  vlist  Take[#, {2, -2}]], 
            Nx,  onezero,  1]];

time[crank[init]]

This computes the solution...

time[cnmat  =  NestList[crank,  init,  Nt - 1];]

... which can be put into an InterpolatingFunction as well.

cnsol  =  
      ListInterpolation[cnmat,  {{0, T}, {-L/2, L/2}}]

It looks quite similar to the one computed by NDSolve..

Plot3D[Abs[cnsol[t, x]], {t, 0, T}, {x, -L/2, L/2},  
        PlotPoints -> 50,  PlotRange -> All];

but there are differences

Plot[{Abs[ndsol[T, x]],  Abs[cnsol[T, x]]}, {x, -L/2, L/2}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]},  
        PlotRange -> All];

The smoothness and the conserved quantity indicates that the Crank-Nicholson solution may well be better:

Timing[
        cnc  =  
          ListPlot[
            Map[Function[{#, 
                  NIntegrate[Abs[cnsol[#, x]]^2,  
                    Evaluate[Join[{x},  cnsol[[3, 2]]]]]}], 
              Range[100]]]][[1]]

One big problem with the Crank-Nicholson method as used above is that you have to set up discretization parameters.  Since we have gone to all the trouble to get an efficient method it is worth spending a few extra minutes to make a function that I could use again sometime.

SchrodingerSolve[initf_,  {ℏ_,  m_,  v_},  time_,  
        space_]  :=  
      Catch[Module[{pack  =  
              Function[
                Developer`ToPackedArray[#,  
                  Complex]], 
t,  begin,  
            end,  Δt,  Nt, 
x,  left,  
            right,  Δx, Nx, 
tlu,  grid,  
            vlist,  init, 
ρ,  p,  spam,  
            onezero,  crank,  
            cnmat}, 
{t, begin,  
              end,  Δt,  Nt}  =  
            GetRange[time]; 
{x,  left,  
              right,  Δx, Nx}  =  
            GetRange[
              space]; 
ρ    = ((I ℏ Δt)/(2  m (((Δx))^2))); 
p  =  
            pack[{ρ,  
                2  -  
                  2 ρ,  ρ}]; 
spam  =  spam  =  {{
Join[Range[Nx - 2],  Range[Nx - 3],  
                    Range[2, Nx - 2]], 
Join[
                    Range[Nx - 2],  Range[2, Nx - 2],  
                    Range[Nx - 3]]}, 
pack[
                  Join[Table[2 + 2 ρ, {Nx  -  2}],  
                    Table[-ρ, {Nx - 3}], 
                    Table[-ρ, {Nx - 
                          3}]]]}; 
grid  =  
            left + Δx Range[0, 
                  Nx - 1]; 
onezero  =  
            pack[{0.   +  0.  I}]; 
If[
            v  ===  0, 
crank  =  
              Function[
                Join[onezero,  tlu[ListCorrelate[p, #]],  
                  onezero]], 
vlist  =  
              pack[Map[
                  Function[Evaluate[x], 
                    Evaluate[I Δt v/ℏ]],  
                  Take[grid, {2, -2}]]]; 
spam[[2, 
                  Range[Nx  -  2]]]  +=  vlist; 
            crank  =  
              Function[
                Join[onezero,  
                  tlu[ ListCorrelate[p, #]  -  
                      vlist  Take[#, {2, -2}] ],  
                  onezero]];
]; 
tlu  =  Experimental`ExtendedLinearSolve[spam]; 
init  =  
            Map[Function[Evaluate[x],  Evaluate[initf]],  
              grid]; 
cnmat  =  
            NestList[crank,  init,  
              Nt - 1]; 
ListInterpolation[
            cnmat, {{begin,  end}, {left,  right}}]]]

GetRange[{var_Symbol,  b_?NumberQ,  e_?NumberQ,  
            n_Integer  /;  n  >  1}]  :=  {var,  b,  e,  
            N[(e  -  b)/(n  -  1)],  n}  /;  
          e  >  b;

SchrodingerSolve::bad  =  "Variable range `1` should be in the \
form {var, begin, end, num}, where var is a symbol, end > begin, and num is \
an integer > 1.";

GetRange[r_]  :=  
        Throw[Message[SchrodingerSolve::bad, 
            r];  $Failed];

Timing[
      SchrodingerSolve[Ψ[x, 0], {ℏ,  m, 0}, {t, 0, T, 
          Nt}, {x, -L/2, L/2, Nx}]]

Timing[
      vsol  =  
        SchrodingerSolve[Ψ[x, 0], {ℏ,  m, V[x]}, {t, 0, T, 
            Nt}, {x, -L/2, L/2, Nx}]]

Adaptive time step

The solution above was computed very quickly and was reasonably good.  The problem is that you may have to try several time steps (and spatial discretizations) to determine convergence (showing how to do that is beyond the scope of this talk.).  One way to handle this is to implement an adaptive time-step.

Double Step

This is a crude, but effective way of setting up an adaptive time step.  It is not the most efficient possible because each step ends up being done doubly.

Definition

This is the module which handles automatic step control.  As arguments, it eccepts a method, a proposed step size, the current solution, the method order, and the local tolerance criteria.  The basic idea is to compute the solution using both one step of size h and two steps of size h/2 and to use Richardson extrapolation to make an estimate of the error. Suppose the method has error of order p.  Then, if (Subscript[A, h])refers to the approximation with step size h, and S the actual solution, we have asymptotically
(Subscript[A, h]) = S  +  c (h^(p + 1))  and (Subscript[A, h/2]) = S  + 2 c (((h/2))^(p + 1))
(p + 1 because this is the local error: combining steps gives one order less, namely p.  Solving these equations for c, the error estimate is
2 ((c(h/2))^p) ≅  ((Subscript[A, h]) - 
              (Subscript[A, h/2]))/((2^p) - 1).
In a similar way, an "optimal" next step size h can also be computed.  In practice we multiply this by a safety factor because it is relatively expensive to reject a step.

The scaled norm here is fairly typical for DE solvers.  The balance of absolute (AccuracyGoal) and relative (PrecisionGoal) tolerance is the same as for the built in NDSolve command.

DoubleStep[meth_, htry_, u_, mord_, {atol_, rtol_}] := 
        Block[{h = 2 htry, half, one, two, errest = Infinity, count = 0, 
            fac = 1}, 
          half = meth[h/2][u]; 
While[
            errest > 1, count++; 
h /= 
              2; 
one = 
              half; 
half = meth[h/2][
                u]; 
two = meth[h/2][
                half]; 
errest = (2^mord - 
                    1) WeightedMaxNorm[atol, rtol, u, two][one, 
                  two];]; 
  If[count < 2, 
            If[Developer`ZeroQ[errest], fac = 2, 
              fac = Min[
                  2, (95/
                        100)*(1/
                          errest)^(1/(mord + 
                              1))]]]; 
{h*fac, h, 
            two}];

This norm is worth a close look at since it is essentially identical to that used in the innards of NDSolve, and is the best way I know to explain how the balance between absolute and relative tolerance is made in NDSolve.

WeightedMaxNorm[atol_, rtol_, wvec__][v1_, v2_] := 
        With[{w = atol + rtol*Map[Max, Transpose[Abs[{wvec}]]]}, 
          Max[Abs[(v1 - v2)/w]]];

Example

SHEuler[h_][{x_,  y_}]  :=  
      With[{newx  =  x  +  h y}, {newx,  y  -  h newx}]

DoubleStep[SHEuler,   .1,  {1. , 0. },  
      1,  {10^-6,  10^-6}]

circle  =  
        Block[{newh  =   .001,  hold  =   .001,  v}, 
          NestWhileList[
Function[{newh,  hold, v}  =  
                DoubleStep[SHEuler,  hold,  #[[2]],  
                  1,  {10^-2,  
                    10^-2}]; 
{#[[1]]  +  hold,  
                v}], 
{0. , {1. , 
                0. }}, 
Function[hold  =  newh; 
              If[#[[1]]  ≥  2 Pi, 
                False, 
If[#[[1]]  +  
                      hold  ≥  2 Pi,  
                  hold  =  
                    2 Pi  -  #[[1]]]; 
True]], 1,  
            10000]];

ListPlot[circle[[All, 2]],  AspectRatio -> 1];

Euler[h_][x_]  :=  
        x  +  h {1, -1} Reverse[x];

spiral  =  
        Block[{newh  =   .001,  hold  =   .001,  v}, 
          NestWhileList[
Function[{newh,  hold, v}  =  
                DoubleStep[Euler,  hold,  #[[2]],  
                  1,  {10^-2,  
                    10^-2}]; 
{#[[1]]  +  hold,  
                v}], 
{0. , {1. , 
                0. }}, 
Function[hold  =  newh; 
              If[#[[1]]  ≥  2 Pi, 
                False, 
If[#[[1]]  +  
                      hold  ≥  2 Pi,  
                  hold  =  
                    2 Pi  -  #[[1]]]; 
True]], 1,  
            10000]];

ListPlot[spiral[[All, 2]],  AspectRatio -> 1];

Accuracy and Precision goals

NDSolve is convenient in that it can take AccuracyGoal and PrecisionGoal options and turn those into proper tolerances, with reasonable defaults.

Here is a function which handles Automatic values for the PrecisionGoal and AccuracyGoal options.  This is essentially identical to how built-in Mathematica functions handle these values.

\(\(AutomaticAGPG[wprec_, apgin_, pora_, \ head_] := 
        Module[{apg = apgin, \ mess}, 
If[
            apg === Automatic, 
apg = 
              If[wprec ≤ $MachinePrecision, 6, 
                wprec - 10]]; 
If[
            Not[NumberQ[apg] && Positive[apg]], 
mess\  = \ 
              If[pora === 
                  Precision, "precg", \ "accg"]; 
MessageEvaluate[head, \ mess, 
              apg]; 
Return[$Failed]]; 
N[
            10^\((\(-apg\))\), wprec]];\)\)

MessageEvaluate[head_,  mess_,  args___]  :=  
      Message[MessageName[head,  mess],  args]

AutomaticAGPG[$MachinePrecision,  Automatic,  Accuracy]

AutomaticAGPG[$MachinePrecision,  -3.1,  Accuracy,  
      nohead]

nohead::accg  =  NIntegrate::accg;
nohead::precg  =  NIntegrate::precg;

SchrodingerSolve with adaptive time step

Since the code is getting a bit more complicated, you can make your life much easier by putting together simple packages.  That way, you don't get errors if you forget to enter a cell or two.

Next, I put the commands together in a package file called Schrodinger.m.  This automatically loads packages for DoubleStep and AutomaticAGPG.

Restart the kernel to avoid any name conflicts.

Needs["Schrodinger`"]

Developer`SetSystemOptions["CatchMachineUnderflow" -> 
        False]

Parameters (just enter to set them for the new kernel)

ℏ = 1.0;
m = 1.0;
σ = 5.0;
τ = 2  (σ^2);
x0 = -50.0;
k0 = 1.0;
L = 200.0;
T = 100.0;

Ψ[x_, 
        t_] := (1  -  
            UnitStep[
              Abs[x]  -  
                100*(1  -  $MachineEpsilon)]) ((((1/(2  π (σ^2)))))^(1/
              4)) (1/(Sqrt[1  +  I (t/τ)])) (E^(I (k0 (x - 
                        x0) - ((ℏ (k0^2))/(2  m)) 
                    t))) (E^(-((((x - (x0 + ((ℏ k0)/m) t)))^2)/(4  (σ^2) (1 + I t/τ)))))

V[x_]  =  100.  UnitStep[x];

Nx = 1001;
Nt = 401;

Examples

Timing[
      csol  =  
        SchrodingerSolve[Ψ[x, 0], {ℏ,  m, V[x]}, {t, 0, T, 
            Nt}, {x, -L/2, L/2, Nx}]]

Timing[
      vsol  =  
        SchrodingerSolve[Ψ[x, 0], {ℏ,  m, V[x]}, {t, 0, T, 
            Automatic}, {x, -L/2, L/2, Nx},  AccuracyGoal -> 3,  
          PrecisionGoal -> 3]]

It takes longer, but at least you have some kind of estimate on the temporal error.  Comparing the two solutions shows a small difference:

Plot[Abs[csol[T, x]] - Abs[vsol[T, x]], {x, -L/2, L/2},  
        PlotRange -> All];

You can see that the time step was decreased when the potential had an effect:

ListPlot[
        vsol[[3, 1]]  -  RotateRight[vsol[[3, 1]]]];

A look at the package

A look at the comments in the  Schrodinger.m. package should highlight some of what I did to combine things together.

Multiple dimensions

An experimental commnd which handles multiple dimensions is also in the package.  Here are the results of a sample run

Timing[
      csol  =  
        SchrodingerSolveMD[
          Exp[I (x  +  y/2)] 
            Exp[-(x/20)^2  -  (y/20)^2], {ℏ,  m, 0}, {t, 
            0, 100, Automatic}, {x, -L/2, L/2, 51}, {y, -L/2, L/2, 
            51},  AccuracyGoal -> 2,  
          PrecisionGoal -> 2]]

Do[Plot3D[csol[t, x, y], {x, -L/2, L/2}, {y, -L/2, L/2},  
        PlotRange -> {0, 1},  PlotPoints -> 50,  
        Compiled -> False], {t, 0, T, 10}]