Numerical ODE and PDE
Solvers in Mathematica.

Introduction and more to NDSolve

NDSolve is Mathematica's general-purpose ODE, and to some extent PDE solver.

ODE's

NDSolve takes equations which specify an ODE or a system of ODE's

A simple first-order ODE

sol  =  
      NDSolve[{y'[x]  ==  y[x], y[0]  ==  1},  y,  {x,  0,  1}]

Plot[Evaluate[y[x]  /.  %],  {x, 0, 1}];

Having the solution returned as an InterpolatingFunction object is very convienient.  You can use the InterpolatingFunction essentially anywhere in Mathematica where you would use any function evaluating to a number.

Plot[Evaluate[(y'[x]  -  y[x])  /.  sol], {x, 0, 1}];

NIntegrate[y[x]  /.  First[sol], {x, 0, 1}]

People often wonder why you have to use the First in the command above -- why wrap the solution in an extra level of List?

NDSolve[{y'[x]^2 ==  y[x], y[0]  ==  1}, y, {x, 0, 1}]

Plot[Evaluate[y[x]  /.  %],  {x, 0, 1}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}];

There are some constraints to keep the method from taking longer than you may have intended for solving a problem.

NDSolve[{y'[x]  ==  y[x], y[0]  ==  1}, y, {x, 0, 250}]

Generally, you can override these constraints by using options

NDSolve[{y'[x]  ==  y[x], y[0]  ==  1}, y, {x, 0, 250}, 
      MaxSteps -> 10000]

StoppingTest

This option allows you to put a conditional on the solution which is checked after every time step.

tab  =  
      Table[First[
          x  /.  NDSolve[{x'[t]  ==  x[t]^2 -  t ,  
                x[0]  ==  z}, x, {t, 0, 20}, 
              StoppingTest -> (Abs[x[t]]  > 2)]], {z, -2, 2,  .1}]

Show[
      Table[Plot[tab[[j]][t], {t, 0, tab[[j, 1, 1, 2]]}, 
          DisplayFunction -> Identity], {j, Length[tab]}], 
      DisplayFunction -> $DisplayFunction, PlotRange -> All];

Accuracy and Precision

Lets see how good the solution is

sol  =  y  /.  First[%]

Lets see how good the solution is at the endpoint

Exp[250. ]  -  sol[250. ]

Plot[
      Evaluate[Log[Abs[Exp[x]  -  y[x]]] /.  First[sol]], {x, 1, 250}];

Plot[Log[10. , Abs[Exp[x]  -  sol[x]]/Exp[x]], {x, 0, 250}];

Controlling the Accuracy and Precision

\(ShowPA[pg_, \ ag_, \ wp_:  16]\  := \ 
      Module[{sol}, 
		
        sol\  = \ 
          First[NDSolve[{\(y'\)[x]\  == \ \(-y[x]\), y[0]\  == \ 1}, 
              y, {x, 0, 5}, \ PrecisionGoal -> pg, \ AccuracyGoal -> ag, 
              WorkingPrecision -> wp]]; 
		sol\  = \ y\  /. \ sol; 
		
        pa\  = \ \(-
            Log[10. , Abs[Exp[\(-x\)]\  - \ sol[x]]/{Exp[\(-x\)], 1}]\); 

        		Plot[Evaluate[pa], \ {x, 0, 5}, PlotRange -> All, \ 
          PlotStyle -> {RGBColor[1, 0, 0], \ RGBColor[0, 0, 1]}]; 
	
        pa\  /. \ x -> 5]\)

ShowPA[Automatic,  Automatic]

ShowPA[10,  Infinity]

ShowPA[Infinity, 10]

ShowPA[Infinity,  20]

ShowPA[Infinity,  20,  32]

This plot does not come out well because Plot by default uses machine precision numbers, which results in catastrophic cancellation for numbers this precise.  The solution is to make sure that the arithmetic needed to produce the data for the plot gets done in higher precision.

\(p[xx_]\  := Block[{x\  = \ SetPrecision[xx, 32]}, \ First[pa]]; 

    a[xx_]\  := \ Block[{x\  = \ SetPrecision[xx, 32]}, \ Last[pa]]; \)

This may take some while.

Plot[{p[x], a[x]}, {x, 0, 5}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}, PlotRange -> All];

Higher Order ODE's

NDSolve handles higher order ODE's naturally

sol  =  
      NDSolve[{x''[t]  +   .1 x'[t]  +  Sin[x[t]]  ==  0, 
          x[0]  == 3,  x'[0]  ==  0}, x, {t, 0, 25}]

ParametricPlot[
      Evaluate[{x[t], x'[t]}  /.  First[sol]], {t, 0, 25}];

sol  =  
      NDSolve[{Derivative[4][x][t] - x[t]  ==  4 Sin[t], 
          x[0]  ==  0,  x'[0]  ==  1,  x''[0]  ==  0,  
          x'''[0]  ==  -3}, x, {t, 0, 25}]

Plot[{t Cos[t],  Evaluate[x[t]  /.  sol]}, {t, 0, 25},  
      PlotStyle -> {RGBColor[1, 0, 0],  RGBColor[0, 0, 1]}];

Systems of ODE's

NDSolve also handles systems of equations, either first order or higher.

\(sol\  = \ 
      NDSolve[{x[t] \(\(θ'\)'\)[t]\  + \ Sin[θ[t]]\  == \ 0, 
          θ[0]\  == \ 1, \(θ'\)[0]\  == \ 0, 
		
          \(\(x'\)'\)[t]\  + \ x[t]\  - \ Cos[θ[t]] - \ 1\  == \ 0, 
          x[0]\  == \ 1.5, \ \(x'\)[0]\  == \ 0}, {x, θ}, {t, 0, 50}]
        \)

Plot[Evaluate[{x[t], θ[t]}  /.  First[sol]], {t, 0, 50},  
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}];

ParametricPlot[
      Evaluate[{x[t], x'[t]}  /.  First[sol]], {t, 0, 50}];

ParametricPlot[
      Evaluate[{θ[t], θ'[t]}  /.  First[sol]], {t, 0, 
        50}];

Boundary Value Problems

NDSolve can also handle linear boundary value problems

\(s\  = \ 
      NDSolve[{\(\(\(y'\)'\)'\)[x] + \ 4\ y[x] == \ 
            \(-8\)*x^3*Cos[x^2]\  + \ 4*Sin[x^2]\  - \ 12*x*Sin[x^2], 

          \ \ \ \ \ \ \(\(y'\)'\)[0] + \(y'\)[0] == 2, \ 
          \(\(y'\)'\)[Pi/2] - y[Pi/2] == 
            2*Cos[Pi^2/4]\  - \ Sin[Pi^2/4]\  - 
\ \ \ \ \ \ 
              Pi^2*Sin[Pi^2/4], \ 
          \(\(y'\)'\)[3\ Pi/4] - \(y'\)[3\ Pi/4] == \ 
            2*Cos[\((9*Pi^2)\)/16]\  - 
\ \ \ \ \ \ 
              \((3*Pi*Cos[\((9*Pi^2)\)/16])\)/2\  - \ 
              \((9*Pi^2*Sin[\((9*Pi^2)\)/16])\)/4}, y, \ {x, 0, Pi}]\)

Stiffness

A problem is stiff when two components vary on extremely different time scales. For example

sol  =  
      NDSolve[{x'[t]  ==  -1000 x[t]  +  Cos[t],  
          y'[t]  ==  -y[t]  +  (t^2),  x[0]  == 1,  
          y[0]  ==  0}, {x, y}, {t, 0, 2}]

Plot[Evaluate[y[t]  /.  sol], {t, 0, 2}];

Plot[Evaluate[x[t]  /.  sol], {t, 0,  .005}];

Advanced Topic: Local versus global error

sol  =  
      NDSolve[{Derivative[4][x][t] - x[t]  ==  4 Sin[t], 
          x[0]  ==  0,  x'[0]  ==  1,  x''[0]  ==  0,  
          x'''[0]  ==  -3}, x, {t, 0, 20}]

Plot[{t Cos[t],  Evaluate[x[t]  /.  sol]}, {t, 0, 20},  
      PlotStyle -> {RGBColor[1, 0, 0],  RGBColor[0, 0, 1]}, 
      PlotRange -> All];

Here we could see that the solution was beginning to diverge from the actual trajectory.  If you don't a-priori know the correct solution, how can you determine if your global error is under control?

Time Reversal

\(end\  = \ 20; 
if\  = \ First[x\  /. \ sol]; 

    init\  = \ 
      Table[\(\(Derivative[j]\)[x]\)[end]\  == \ 
          \(\(Derivative[j]\)[if]\)[end], {j, 0, 3}]; 

    rsol\  = \ 
      First[NDSolve[{\(\(Derivative[4]\)[x]\)[t] - x[t]\  == \ 4\ Sin[t], 
            init}, x, {t, end, 0}]]\)

Plot[Evaluate[{x[t]  /.  sol, x[t]  /.  rsol}], {t, 0, end},  
      PlotStyle -> {RGBColor[1, 0, 0],  RGBColor[0, 0, 1]}, 
      PlotRange -> All];

\(\(err\  = \ 
      Function[{end}, if\  = \ First[x\  /. \ sol]; 

        init\  = \ 
          Table[\(\(Derivative[j]\)[x]\)[end]\  == \ 
              \(\(Derivative[j]\)[if]\)[end], {j, 0, 3}]; 

        rsol\  = \ 
          x\  /. \ First[
              NDSolve[{\(\(Derivative[4]\)[x]\)[t] - x[t]\  == \ 4\ Sin[t], 
                  init}, x, {t, end, 0}]]; {end, if[0]\  - \ rsol[0]}]; 
    \)\)

errlist  =  err /@ Drop[if[[3, 1]], 1];

ListPlot[errlist, PlotRange -> All];

Needs["Graphics`Graphics`"]

LinearLogListPlot[Abs[errlist]]

Perturb the initial condition

psol  =  
      NDSolve[{Derivative[4][x][t] - x[t]  ==  4 Sin[t], 
          x[0]  ==  (10^(-6)),  x'[0]  ==  1,  
          x''[0]  ==  0,  x'''[0]  ==  -3}, 
        x, {t, 0, 20}]

Plot[
      Evaluate[(x[t]  /.  sol)  -  (x[t]  /.  psol)], {t, 0, 
        20}];

LinearLogPlot[
      Evaluate[Abs[(x[t]  /.  sol)  -  (x[t]  /.  psol)]], {t, 
        0, 20}];

Conserved Quantities

The above example does not have any clearly conserved quantities, so I'll use another one that does...the nonlinear Schrodinger equation (NLSE)

NLSE[u_[x_, t_]]  :=  
      (I D[u[x, t], t]  ==  
          D[u[x, t], {x, 2}]  +  
            2 \[LeftBracketingBar]u[x, t](\[RightBracketingBar]^2) u[x, t])

We will start with a nice innocent looking periodic initial condition...

({{({{(Subscript[f, n])_[x_]  :=  
                    a (1  +  ϵ Cos[(Subscript[μ, n]) x])}, {ϵ  =  1/10}, {L  =  4 (Sqrt[2]) π}})}, {a  =   .35}, {(Subscript[μ, n])_ :=  2  π n/L}, {T  =  250}});

Plot[(Subscript[f, 1])[x], {x, 0, L}];

u  =   
      NDSolve[{NLSE[ϕ[x, t]], ϕ[x, 0]  ==  (Subscript[f, 1])[x],  
            ϕ[0, t]  ==  ϕ[L, t]}, ϕ, {x, 0, L}, {t, 0, T}, 
          StartingStepSize -> {L/32, Automatic}, MaxSteps -> 20000][
        [1, 1, 2]];

Take a look at the value of the evolved solution at x = 0.

Plot[Abs[u[0, t]], {t, 0, T}, PlotRange -> All, PlotPoints -> 100];

The NLSE has infinitely many conserved quantities to check.  One of the simplest is \(TraditionalForm
      \`M = ∫\_0\%L
            \(|u|\^2\) 
            dx\) which we can approximate by

Length[spaces  =  u[[3, 1]]]Length[times  =  u[[3, 2]]]

M[u_][t_]  :=  
      NIntegrate[Abs[u[x, t]]^2, Evaluate[{x, Sequence @@ u[[3, 1]]}]]

Timing[M[u][0]]

Timing[M[u][T]]

\(Clear[M]; 
    \(M[u_]\)[t_]\  := \ 
      Module[{grid\  = \ Drop[u[\([3, 1]\)], \(-1\)], \ h}, 
		
        h\  = \ grid[\([2]\)]\  - \ grid[\([1]\)]; 
		
        \((Plus\ @@\ \(Abs[u[grid, t]]\^2\))\)*h]\)

Timing[M[u][0]]

ListPlot[M[u] /@ times];

Try higher precision

sol16  =  First[x  /.  sol]

hsol[prec_] :=  
      First[x  /.  
          NDSolve[{Derivative[4][x][t] - x[t]  ==  4 Sin[t], 
              x[0]  ==  0,  x'[0]  ==  1,  x''[0]  ==  0, 
               x'''[0]  ==  -3}, x, {t, 0, 20},  
            WorkingPrecision -> prec,  MaxSteps -> 10000]]

Timing[sol20  =  hsol[20]; ]

sol16[20]  -  sol20[20]

Timing[sol25  =  hsol[25]; ]

sol20[20]

sol25[20]

Timing[sol30  =  hsol[30]; ]

sol30[20]

sol20[20]  -  sol30[20]

PDE's

Basic

\(sol = 
      u /. First[
          NDSolve[{
              ∂\_t u[x, t] == 
                \((2\ ∂\_{x, 2}u[x, t])\)/\((9\ π^2)\), 
              u[x, 0] == Sin[\ 5  π\ x/2], u[0, t] == Sin[t], 
              \(\(Derivative[1, 0]\)[u]\)[1, t] == 0}, 
            u, {x, 0, 1}, {t, 0, 6  π}]]\)

Plot3D[sol[x, t], {x, 0, 1}, {t, 0, 6  π}, PlotRange -> All, 
      PlotPoints -> 25];

Or systems of PDE's

\({usol, vsol} = {u, v} /. 
        First[NDSolve[{
              ∂\_t u[t, x] == 
                ∂\_x\((\((v[t, x] - 1)\)\ ∂\_x u[t, x])
                      \) + \((16  x\ t - 2  t - 16\ \((v[t, x] - 1)\))\)\ 
                    \((u[t, x] - 1)\) + 10  x\ Exp[\(-4\)\ x], 
              ∂\_t v[t, x] == 
                ∂\_{x, 2}v[t, x] + ∂\_x u[t, x] + 
                  4\ u[t, x] - 4 + x^2 - 2  t - 10  t\ Exp[\(-4\)\ x], 
              u[0, x] == 1, v[0, x] == 1, u[t, 0] == 1, v[t, 0] == 1, 
              3\ u[t, 1] + \(\(Derivative[0, 1]\)[u]\)[t, 1] == 3, 
              5\ \ \(\(Derivative[0, 1]\)[v]\)[t, 1] == 
                Exp[4]\ \((u[t, 1] - 1)\)}, {u, v}, {x, 0, 1}, {t, 0, 2}]]
          \)

\(ContourPlot[usol[t, x], {x, 0, 1}, {t, 0, 2}]; 

    ContourPlot[vsol[t, x], {x, 0, 1}, {t, 0, 2}]; \)

Limitations

What the method cannot handle is equations that cannot be properly posed as a Cauchy problem

sol = 
      u /. First[
          NDSolve[{
              (∂_(x, x))u[x, y] + (∂_(y, y))u[x, y] == 
                0, u[x, 0] == Cos[2  π x], 
              Derivative[0, 1][u][x, 0] == 0, u[0, y] == 1, 
              u[1, y] == 1}, u, {x, 0, 1}, {y, 0, 1}]]

Plot3D[sol[x, y], {x, 0, 1}, {y, 0, 1}];

At the current time, it also is not able to handle PDEs which form discontinuities.

\(sol = 
      u /. First[
          NDSolve[{
              ∂\_t u[x, t] == \(-u[x, t]\)\ ∂\_x u[x, t], 
              u[x, 0] == Sin[2  π\ x], u[0, t] == u[1, t]}, {u}, {x, 0, 
              1}, {t, 0,  .5}]]; 
Plot3D[sol[x, t], {x, 0, 1}, {t, 0,  .5}]; 
    \)

DifferenceOrder

Airy's equation contains a third derivative term, which is hard to resolve accurately with lower order finite differences.

\(sol = 
      u /. First[
          NDSolve[{∂\_t u[x, t] == \(-∂\_{x, 3}u[x, t]\), 
              u[x, 0] == Exp[\(-x^2\)], u[\(-5\), t] == u[5, t]}, {u}, {t, 
              0,  .4}, {x, \(-5\), 5}]]; 

    Do[Plot[sol[x, t], {x, \(-5\), 5}, PlotRange -> {\(-1\), 1}], 

      	{t, 0,  .4,  .025}]; \)

Try specifying sixth order differences.  Note how you can specify different values for the spatial and temporal parts.

\(sol = 
      u /. First[
          NDSolve[{∂\_t u[x, t] == \(-∂\_{x, 3}u[x, t]\), 
              u[x, 0] == Exp[\(-x^2\)], u[\(-5\), t] == u[5, t]}, {u}, {x, 
              \(-5\), 5}, {t, 0,  .4}, 
            DifferenceOrder -> {6, \ Automatic}]]; 

    Do[Plot[sol[x, t], {x, \(-5\), 5}, PlotRange -> {\(-1\), 1}], 

      	{t, 0,  .4,  .025}]; \)

Balancing Spatial and Temporal error tolerances

u1 =  
      First[ϕ  /.   
          NDSolve[{NLSE[ϕ[x, t]], ϕ[x, 0]  ==  (Subscript[f, 1])[x],  
              ϕ[0, t]  ==  ϕ[L, t]}, 
            ϕ, {x, 0, L}, {t, 0, T}, 
            StartingStepSize -> {L/32, Automatic}, MaxSteps -> 20000, 
            AccuracyGoal -> {4,  6},  PrecisionGoal -> {4, 6}]];

Plot[Abs[u1[0, t]], {t, 0, T}, PlotRange -> All, PlotPoints -> 100];

ListPlot[M[u1] /@ u1[[3, 2]]];

User Defined grid

\(\(Timing[
        sol = u /. 
            \((\(NDSolve[{
                    ∂\_t u[x, t] == 
                      \(-u[x, t]\)\ ∂\_x u[x, t] +  .01\ 
                          ∂\_{x, 2}u[x, t], 
                    u[x, 0] == Sin[2  π\ x], u[0, t] == u[1, t]}, {u}, {
                    x, 0, 1}, {t, 0,  .5}]\)[[1
                ]])\)]\)[[1
      ]]\)

Plot[sol[x,  .5], {x, 0, 1}];

\(grid\  = \ sol[\([3, 1]\)]; 

    ListPlot[Transpose[{grid, \ sol[grid, \  .5]}]]; \)

Clear[xgrid]; \(xgrid[n_, bunch_] := 
      Table[N[k/n + \((bunch\ Sin[\((2  π\ k)\)/n])\)/\((2  π)\)], {
          k, 0, n}]; 
ListPlot[xgrid[Length[grid], \  .9]]; \)

\(\(Timing[
        sol = u /. 
            \((\(NDSolve[{
                    ∂\_t u[x, t] == 
                      \(-u[x, t]\)\ ∂\_x u[x, t] +  .01\ 
                          ∂\_{x, 2}u[x, t], 
                    u[x, 0] == Sin[2  π\ x], u[0, t] == u[1, t]}, {u}, 
                  Join[{x}, xgrid[Length[grid],  .9]], {t, 0,  .5}]
                  \)[[1]])\)]
        \)[[1]]\)

\(grid\  = \ sol[\([3, 1]\)]; 

    ListPlot[Transpose[{grid, \ sol[grid, \  .5]}]]; \)

\(\(Timing[
        sol = u /. 
            \((\(NDSolve[{
                    ∂\_t u[x, t] == 
                      \(-u[x, t]\)\ ∂\_x u[x, t] +  .01\ 
                          ∂\_{x, 2}u[x, t], 
                    u[x, 0] == Sin[2  π\ x], u[0, t] == u[1, t]}, {u}, 
                  Join[{x}, xgrid[30,  .9]], {t, 0,  .5}]
                  \)[[1]])\)]
        \)[[1]]\)

\(grid\  = \ sol[\([3, 1]\)]; 

    ListPlot[Transpose[{grid, \ sol[grid, \  .5]}]]; \)

How does it work

As I have mentioned before, it is based on the numerical method of lines.  The way this works is by choosing one of the variables to be "spatial" and the other to be "temporal".  Sometimes this choice is arbitrary, but typically it is based on given intial conditions and the differential order with respect to the two variables.

Once the types of the two variables are chosen, the equations are discretized w.r.t. the "spatial" variable.

I have loaded a special version of the PDE routines which show you the discretizations that are used (sorry, but this is not available with the downloaded version.)

NDSolve[{D[u[x, t], t]  ==  D[u[x, t], {x, 2}], 
        u[x, 0]  ==  Exp[-(x^2)],  u[-10, t]  ==  Exp[-100], 
         u[10, t]  ==  Exp[-100]}, u, {x, -10, 10}, {t, 0, 1}]

You may be more familar with the second order difference formula:

NDSolve[{D[u[x, t], t]  ==  D[u[x, t], {x, 2}], 
        u[x, 0]  ==  Exp[-(x^2)],  u[-10, t]  ==  Exp[-100], 
         u[10, t]  ==  Exp[-100]}, u, {x, -10, 10}, {t, 0, 1}, 
      DifferenceOrder -> {2,  Automatic}]

Notice that in this case the boundaries are much simpler because of the narrower stencil of the second order method.

Nonlinearities are currently handled in the obvious way. (Currently no attempt is made to use conservative differencing, though in the future we will be working on this.)

NDSolve[{
        D[u[x, t], t]  ==  u[x, t] D[u[x, t], x] D[u[x, t], {x, 2}], 
        u[x, 0]  ==  Exp[-(x^2)],  u[-10, t]  ==  Exp[-100], 
         u[10, t]  ==  Exp[-100]}, u, {x, -10, 10}, {t, 0, 1}, 
      DifferenceOrder -> {2,  Automatic}]