NDSolve Method Extensibility

Writing your own methods for fun and science ...

Rob Knapp

Wolfram Research, Inc.

Initialization

Off[General::spell];
Off[General::spell1];

$Path  =  
        Join[$Path,  {"/Documents and Settings/rknapp/My \
Documents/Talks/2004 Tech Conference/Packages"}];

Needs["GUIKit`"]

Needs["NumericalMath`EquationTrekker`"]

Needs["NumericalMath`OrderStar`"]

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Needs["DifferentialEquations`NDSolveUtilities`"]

SetOptions[InvariantErrorPlot,  
        InvariantErrorFunction -> Subtract];

MaxTime[ifun_InterpolatingFunction]  :=  
      Max[First[InterpolatingFunctionDomain[ifun]]]

Needs["Impulse`"]

Heun[]@"Step"[rhs_,  t_,  h_,  y_,  yp_]  :=  
      Module[{k1  =  
            y  +  h yp}, 
{h,  (h/2) (yp  +  
                rhs[t  +  h,  k1])}]

penergy  =  y[t]^2 - 2 Cos[x[t]];

For Mathematica version 5, the differential equation solver, NDSolve, was completely redesigned.  One of the important features of the design is a framework which includes combination of methods and extensibility.  The NDSolve framework handles processing equations, calling methods to advance the solution, and construction of a solution.

A basic use of the framework is direct incorporation of particular methods.

NDSolve framework ⟶ ExplicitRungeKutta

pendulum  =  
        Apply[Sequence,  {{x'[t]  ==  
                y[t], y'[t]  +  Sin[x[t]]  ==  0,  
              x[0]  ==  1,  y[0]  ==  0},  {x, y}, {t, 0, 
              800}}];
penergy  =  y[t]^2 - 2 Cos[x[t]];

InvariantErrorPlot[penergy, {x[t], y[t]},t,
NDSolve[ pendulum, 
            Method -> "ExplicitRungeKutta"]];

The framework allows for nesting of methods inside of controller methods.

NDSolve framework  ⟶ Extrapolation ⟶ ExplicitMidpoint

InvariantErrorPlot[penergy, {x[t], y[t]},t,
NDSolve[ pendulum, 
            Method -> {"Extrapolation",  
                Method -> "ExplicitMidpoint"}]];

NDSolve framework  ⟶ Projection ⟶ ExplicitRungeKutta

InvariantErrorPlot[penergy, {x[t], y[t]},t,
NDSolve[ pendulum, 
            Method -> {"Projection",  Invariants -> penergy, 
                Method -> "ExplicitRungeKutta"}]];

The method invocation can be arranged as a tree.

NDSolve framework ({{⟶, Splitting f=(Subscript[f, 1])+(Subscript[f, 2])  ({{↗, LocallyExact  (Subscript[f, 1])}, {⟶, ({{ExplicitEuler (Subscript[f, 2]), }})}, {↘, ({{LocallyExact (Subscript[f, 1]), }})}})}})

StrangSplittingCoefficients[2,  p_] :=  
        N[{1/2, 1, 1/2},  p];

InvariantErrorPlot[penergy, {x[t], y[t]},t,
\(NDSolve[pendulum, 
            StartingStepSize -> 1/8, 
\ 
            Method -> {"Splitting", \ "Coefficients" -> 
                  StrangSplittingCoefficients, \ "DifferenceOrder" -> 2, \ "Equations" -> {{\(x'\)[t]\  == \ 
                        y[t], \ \(y'\)[t]\  == \ 
                        0}, {\(x'\)[t]\  == \ 
                        0, \(y'\)[t] == \(-Sin[x[t]]\)}, {\(x'\)[
                          t]\  == \ y[t], \ \(y'\)[t]\  == \ 
                        0}}, 
"Method" -> {"LocallyExact", \ "ExplicitEuler", \ "LocallyExact"}}]\)];

The NDSolve framework includes extensibility to methods which can be "plugged-in" through a particular interface.  Within the framework a method is initialized by name as a data object with that name as head.  Method invocation is handled by object "methods" (or subvalues associated with the method name).

mymethod[data]@Step[rhs, t, h, y, yp]] (invokes method to take a step) -> new step data

It is important to keep in mind that it is the solution increment, Δy, such that y(t + h) = y(t) + Δy not the solution y(t+h) itself which is returned.  Having the solution increment as the return allows NDSolve to better control roundoff error through compensated addition.

This is sort of like a "Hello World" example for differential equation methods.  It is simple and yet does enough to illustrate the principles. Note, this method is sometimes referred to as the extended Euler and sometimes as the midpoint method.  It is a two-stage Runge-Kutta method:

This defines a Step function which computes a step of Heuns's method.

Heun[]@"Step"[rhs_,  t_,  h_,  y_,  yp_]  :=  
      Module[{k1  =  
            y  +  h yp}, 
{h,  (h/2) (yp  +  
                rhs[t  +  h,  k1])}]

This computes the solution of the pendulum equation using Heun's method with fixed step size.

hsol  =  
      NDSolve[{x'[t]  ==  
            y[t],  y'[t]  +  Sin[x[t]]  ==  0,  
          x[0]  ==  3,  y[0]  ==  0}, {x, y},  {t, 0, 
          100},  
Method -> Heun,  
        StartingStepSize -> 1/8]

Here is a plot of the solution.

ParametricPlot[{x[t],  y[t]}  /.  First[hsol], {t, 0, 
          100}];

The change in orbits apparent near the separatrices is because of energy loss.  This is more apparent from the invariant plot.

InvariantErrorPlot[penergy,  {x[t], y[t]}, t, hsol];

Part of the problem above is that the problem changes its rate of variation tremendously depending on whether it is near the separatrices or not.  A fixed step size would have to be small enough to resolve the fastest variation, which might be prohibitive over the entire range.  An adaptive step size typically allows the whole integration to be done somewhat more efficiently.

NDSolve provides a simple extrapolative step controller called DoubleStep, which estimates error based on extrapolation of solution approximations found by using one step of size h versus two of size h/2.  The error estimate allows for adaptive step control.  The extrapolation has the further advantage of using the more accurate step approximations, effectively giving an extra order in asymptotic convergence.

For the extrapolation to work properly, it is necessary to know the difference order of the base method.  NDSolve gets this information by querying the method object for the "DifferenceOrder" property.

This sets the "DifferenceOrder" property for Heun's method to the correct value of 2.

Heun[]@"DifferenceOrder"[]  =  2;

The DoubleStep controller also checks for the possibility that the equations are too stiff to effectively integrate with the method.  To do this, it also needs to know the linear stability boundary of the base method along the negative real axis.

This loads a package which aids in plotting the stability regions for various methods.

Needs["NumericalMath`OrderStar`"]

This computes the polynomial in hλ for the linear test equation x' = λ x.

poly  =  
      Expand[1  + 
          Last[Heun[]@"Step"[Function[{t, x}, λ x],  t,  h,  
                  x,  λ x]]/x]

This plots the stability region, the region in the hλ plane where the method is contractive: i.e. |poly| < 1

OrderStar[poly  /. h -> 1,  1,  
        Ticks -> True];

This sets the "LinearStabilityBoundary" property to the value where the stability region crosses the negative real axis, -2.

Heun[]@"LinearStabilityBoundary"[]  =  -2;

This computes the solution of the pendulum equation with adaptive time steps.

dhsol  =  
      NDSolve[{x'[t]  ==  
            y[t],  y'[t]  +  Sin[x[t]]  ==  0,  
          x[0]  ==  3,  y[0]  ==  0}, {x, y},  {t, 0, 
          100},  
Method -> {DoubleStep, 
            Method -> Heun}]

Here is a plot of the solution as far as it got.

ParametricPlot[{x[t],  y[t]}  /.  First[dhsol], {t, 0, 
          72}];

Here is a plot of the error in the energy invariant.

InvariantErrorPlot[penergy,  {x[t], y[t]}, t, dhsol];

Here is a plot of the step sizes the adaptive step controller used combined with a plot of x (scaled and shown in red).

Show[
        Block[{$DisplayFunction  =  Identity}, {StepDataPlot[dhsol,  
              PlotJoined -> False], 
            Plot[Evaluate[ .001 {x[t], x'[t]}  /.  dhsol], {t, 0, 
                72},  PlotStyle -> {Red,  Blue}]}]];

The leapfrog method,

(Subscript[x, n + 1]) = 
        (Subscript[x, n - 1])  +  2 hf((Subscript[t, n]), (Subscript[x, n]))
 
 is one commonly used for strictly hyperbolic PDE's because it achieves second order with only one evaluation per step and has reasonable stability on the imaginary axis.  (Alas, it has no stability at all for dissipative terms.)

The leaping derivative approximation using (Subscript[x, n - 1]) requires an initialization step.  This can be done with any method second order or higher.  In the first implementation, this will be done by the Heun's method.  Because the method uses data ((Subscript[x, n - 1])) which is not directly provided in the arguments to the Step function, it will need to save data.

The initial step is handled by Heun's method.  The returned value needs to be extended to include the data object since it is modified by storing the value which, for the next step is (Subscript[x, n - 1])

This defines the "Step" function for the initial step.  It returns a new Leapfrog method object which

Leapfrog[]@"Step"[rhs_,  t_,  h_,  y_,  yp_]  :=  
        Module[{yinc}, 
yinc  = 
            Part[Heun[]@"Step"[rhs,  t,  h,  y,  yp], 
              2]; 
{h,  yinc,  None,  
            Leapfrog[yinc]}];

This defines the "Step" function for all subsequent steps.

Leapfrog[yincprev_?VectorQ]@"Step"[rhs_,  t_,  h_,  y_,  
            yp_]  :=  
Module[{yinc  =  
              2 h yp  -  yincprev}, {h, yinc,  None,  
            Leapfrog[yinc]}];

This computes the solution to the pendulum equation using the Leapfrog method.  (This is stable where the Jacobian has purely imaginary eigenvalues as long as h is small enough.

lsol  =  
      NDSolve[{x'[t]  ==  
            y[t],  y'[t]  +  Sin[x[t]]  ==  0,  
          x[0]  ==  3,  y[0]  ==  0}, {x, y},  {t, 0, 
          100},  
Method -> Leapfrog,  
        StartingStepSize -> 1/8]

Here is a plot of the solution.

ParametricPlot[{x[t],  y[t]}  /.  First[lsol], {t, 0, 
          100}];

Here is a plot of the energy invariant error.

InvariantErrorPlot[penergy,  {x[t], y[t]}, t, lsol];

Of course, since the method is usually used for PDEs, a PDE example is in order...
Keep in mind that the second order method is really only cheap if you are looking for low overall accuracy in the solution (which is typical for PDE solutions).  Usually, it will make sense to have the spatial discretization order match the second order of the leapfrog method.

The sine-Gordon equation is hyperbolic and so it can be integrated with the leapfrog method.

This computes the solution of the sine-Gordon equation with periodic boundary conditions using the leapfrog method with a second order spatial discretization.  The solution accuracy is determined by a combination of the spatial grid spacing (0.1) and the temporal grid spacing (0.05).  The temporal spacing cannot be increased for this particular spatial grid because of stability restrictions.  One can think of the stability restriction in terms of the CFL criteria.

\(sgsol\  = \ 
      NDSolve[{D[u[t, x], t, t]\  == \ 
            D[u[t, x], x, x]\  - \ Sin[u[t, x]], \ 
          u[0, x]\  == \ 
            Exp[\(-x^2\)], \ \((D[u[t, x], t]\  /. \ 
                t -> 0\ )\) == \ 0, \ 
          u[t, \(-5\)]\  == \ u[t, 5]}, u, {t, 0, 5}, {x, \(-5\), 5}, \ 
        Method -> {"MethodOfLines", \ 
"SpatialDiscretization" -> {"TensorProductGrid", \ "DifferenceOrder" -> 
                  2, \ "MaxStepSize" ->  .1, \ "MinStepSize" ->  .1}, 
Method -> 
              Leapfrog}, 
StartingStepSize -> 
          0.05]\)

Here is a surface plot of the solution.

Plot3D[
        First[u[t, x]  /.  sgsol], {t, 0, 5}, {x, -5, 5}];

To demonstrate handling method options along with method initialization, the implementation of the Leapfrog method will be extended to include an option specifying the starting method.

Method initialization is done when NDSolve actually needs to use the method by evaluating NDSolve`IntializeMethod with arguments according to a specific template:

This sets the default option for the Leapfrog method.

\(\(Options[
          Leapfrog]\  = \ {StartingMethod -> "ExplicitMidpoint"};\)\)

TagSet is used to attach the initialization rule to Leapfrog, thus not slowing down initialization of other methods.  NDSolve`InitializeSubmethod is a utility function which is convenient for initializing a submethod.

\(\(Leapfrog\  /: 
NDSolve`InitializeMethod[Leapfrog, \ 
          stepmode_, \ rhs_, \ 
          state_, \ {opts___}]\  := \ 
Module[{startmethod, \ startmethodobject}, 
Print["Initializing Leapfrog with method options ", \ {opts}]; 
startmethod\  = \ 
            StartingMethod\  /. \ 
              Flatten[{opts, \ 
                  Options[
                    Leapfrog]}]; 
startmethodobject\  = \ 
            NDSolve`InitializeSubmethod[startmethod, \ stepmode, \ rhs, \ 
              state]; 
Leapfrog[
            startmethodobject]];\)\)

Now the initial step is computed simply by calling the "Step" function of the submethod.  Then the Leapfrog method object with the increment is returned so that the leapfrog step will be used for all subsequent steps.

Leapfrog[startmethodobject_]@"Step"[rhs_,  t_,  h_,  y_,  
            yp_]  :=  Module[{yinc}, 
yinc  = 
            Part[startmethodobject@"Step"[rhs,  t,  h,  y,  yp], 
              2]; 
{h,  yinc,  None,  
            Leapfrog[yinc]}];

You give method options to your method just like I have shown above for built-in methods.

This integrates the sine-Gordon equation using Leapfrog initialized by a third order explicit Runge-Kutta method.

\(sgsol\  = \ 
      NDSolve[{D[u[t, x], t, t]\  == \ 
            D[u[t, x], x, x]\  - \ Sin[u[t, x]], \ 
          u[0, x]\  == \ 
            Exp[\(-x^2\)], \ \((D[u[t, x], t]\  /. \ 
                t -> 0\ )\) == \ 0, \ 
          u[t, \(-5\)]\  == \ u[t, 5]}, u, {t, 0, 5}, {x, \(-5\), 5}, \ 
        Method -> {"MethodOfLines", \ 
"SpatialDiscretization" -> {"TensorProductGrid", \ "DifferenceOrder" -> 
                  2, \ "MaxStepSize" ->  .05, \ "MinStepSize" ->  .05}, 
Method -> {Leapfrog, \ 
                StartingMethod -> {ExplicitRungeKutta, \ "DifferenceOrder" -> 
                      3}}}, 
StartingStepSize -> 
          0.025]\)

As a final example, I'll show you a method that works as a controller method which is all defined in a package.

The impulse method does some symbolic preprocessing of the equations to find any  δ (DiracDelta) functions in them and sets up the zeroes of the δ functions as events for the EventLocator method.  When a zero event is detected in a step, the Impulse controller modifies the step size to integrate almost, but not quite all the way, up to the zero.  Then the impulse is computed, based on the derivative of the argument to the δ function, and it is provided with a step of small width sufficient to move just past the zero.  Then integration continues as before.

Before showing you the code, I'll demonstrate how is works by showing an example.

This loads the Impulse method package.

Needs["Impulse`"]

Here is a system which was submitted by user Mats Jirstrand.

system  =  
        Apply[Sequence, {
{x'[t] + 
                  x[t] - (1/2) y[t] == 
                (1/5), 
y'[t] + 
                  y[t] - (3/10) x[t] == 
                0 + (1/2) DiracDelta[t - 5] - y[t] DiracDelta[t - 7] + 
                  DiracDelta[t - 11], 
x[0] == 
                1, 
y[0] == 
                0}, 
{x, y}, {t, 0, 12}}];

This integrates the system with the default method and plots the solutions.

sol = NDSolve[system];
Plot[Evaluate[{x[t], y[t]}  /.  First[sol]], {t, 0, 12},  
        PlotRange -> All, 
        PlotStyle -> {Red,  Blue}];

This integrates the system with the impulse method and plots the solutions.

sol = 
        NDSolve[system,  
          Method -> Impulse];
Plot[Evaluate[{x[t], y[t]}  /.  First[sol]], {t, 0, 12},  
        PlotRange -> All, 
        PlotStyle -> {Red,  Blue}];

The code for the impulse method is in a package file Impulse.m.  I show most of it here in an annotated form for illustration.

The only method option is the submethod to use for integrating the system.

Options[Impulse]={"Method"->"StiffnessSwitching"};

Method initialization.  This ultimately returns an Impulse method object with needed data saved.

Impulse/:NDSolve`InitializeMethod[Impulse,stepmode_,nf_,state_,rawmopts_]:=
Module[{submethod,submethoddata,tvar,dvars,fun,fatnf,cnf,hasdep,deltaevents,deltatime,prec},

This gets the rhs function.  It is generally a Function expression, so to get something that can be scanned symbolically, it is applied to the independent and dependent variables.

      fun=nf@"FunctionExpression"[];
      tvar=state@"TemporalVariable"[];
      dvars=nf@"ArgumentNames"[];
      fun=fun[tvar,Apply[Sequence,dvars]];

The function hasdep will be used to check for the presence of dependent variables in the δ functions, which this method cannot handle.
      
      Map[Function[hasdep[#]:=Throw[Message[NDSolve::ddvar,#];$Failed]],dvars];

 The function is scanned for the presence of any δ functions.  For each one found, the argument is added to the list of events.
      
      deltaevents = Cases[
		fun,
		Rule[
			Condition[DiracDelta[z_],And[MemberQ[z,tvar,Infinity],
			(Scan[hasdep,{z},-2];True)]],z],
		Infinity];

Get the submethod from the option and set up an EventLocator method using that method for integration.  
The events are the arguments of the δ functions.  The action to take upon an event is to set a variable deltatime to the event time.  Note: this is only set if it has not previously been set by another event detected in the same step so that it catches the first event in the integration direction.  The Function substitution is used to set up the assignment expression without evaluating it. 
      
      submethod = Method /. Flatten[{rawmopts, Options[Impulse]}];
      submethod={"EventLocator",
      	"Event"->deltaevents,
      	"EventAction":>(#1[MatchQ[deltatime,_Symbol],#2[deltatime,#3]]),
      	"Method"->submethod
      	}&[If,Set,tvar];
      submethoddata=NDSolve`InitializeSubmethod[submethod,stepmode,nf,state];

The size of the interval used around the δ functions is chosen to be half of the precision.  This is consistent with the default tolerances for EventLocator.  Really, this should be an option and/or ties to AccuracyGoal, PrecisionGoal.
      
      prec=nf@"WorkingPrecision"[];
      width=SetPrecision[10.^(-prec/2),prec];

Set up numerical functions 
cnf which gives the rhs of the ODE without the δ functions
fatnf which gives the contribution of the delta functions which have zeros within the interval determined above.

The "ProcessExpression" method of NDSolve`StateData gives a NumericalFunction object with arguments based on the variables given to NDSolve.
      
      size=nf@"ResultDimensions"[];
      fatnf=fun/.DiracDelta[z_]:>(If[Abs[z]≤width,rdz,0]/.rdz->1/D[z,tvar]);
      cnf=fun/.DiracDelta[z_]->0;
      fatnf=Expand[fatnf-cnf];
      fatnf=state@"ProcessExpression"[All,fatnf,size];
      cnf=state@"ProcessExpression"[All,cnf,size];

Finally (whew!) the Impulse method object is returned.
      
      Impulse[{"Step",submethoddata,deltatime,cnf,fatnf,width}]];

The Impulse method has three modes.  The usual mode is the "Step" mode where it simply takes a step of the integration method and checks to see if an event (δ zero) was present.

Impulse[{"Step", submethod_, deltatimesymbol_, cnf_, fatnf_, eps_}]["Step"[rhs_, t_, h_, y_, yp_]] := 
Module[{tnext, newh, dy, newyp, newsub, dt, deltatime},

This invokes the EventLocator method (which invokes the integration method).  Note that the rhs with no δ functions is used to avoid the unpleasant consequences of evaluating at a zero of one the δ functions.
      
    invoke = NDSolve`InvokeMethod[submethod, cnf, t, h, y, yp];
    If[!ListQ[invoke], Throw[$Failed]];
    {newh, yinc, newyp, newsub} = invoke;

Unset the deltatimesymbol for the next step.  
      
    deltatime = deltatimesymbol;
    Unset[deltatimesymbol];
    If[NumberQ[deltatime],

If deltatime was set to a number, this was done by the action for the first event which was found in the step.   We return a step rejection ($Failed in second arg) with a new requested step to bring us just short of the event (δ zero).  The Impulse object returned with deltatime as first argument instead of "Step" to change modes.
      
        {deltatime - t - Sign[h] eps, $Failed, None, 
        	Impulse[{deltatime, submethod, deltatimesymbol, cnf, fatnf, eps}]},

Else, return the step data as given by the integration method.
      
        {newh, yinc, newyp, 
        	Impulse[{"Step", newsub, deltatimesymbol, cnf, fatnf, eps}]}]
];

The next mode is used to step just short of the δ zero found and to set the stepsize to leap over the δ zero.

Impulse[{deltatime_?NumberQ, submethod_, deltatimesymbol_, cnf_, fatnf_, eps_}]["Step"[rhs_, t_, h_, y_, yp_]] := 
Module[{tnext, newh, dy, newyp, newsub, dt, impulse, invoke},

Invoke the integration method.  
      
    invoke = NDSolve`InvokeMethod[submethod, cnf, t, h, y, yp];
    If[!ListQ[invoke], Throw[$Failed]];
    {newh, yinc, newyp, newsub} = invoke;
    Unset[deltatimesymbol];
    If[Not[VectorQ[yinc]],

Step rejection: go back to "Step" mode with the shorter step.  This is unlikely to happen because the method already took a larger step.  However, it happens occasionally when the smaller step allows the method to resolve variations that were missed with a larger step.   
      
        {newh, yinc, newyp, Impulse[{impulse, newsub, deltatimesymbol, cnf, fatnf, eps}]},
 
Else, compute the impulse by evaluating the version of the rhs with the delta functions having zero near this point give their proper contribution.  The stepsize is very smal: just enough to jump as far past the δ zero as we stopped short of it.  The first argument is the computed impulse which is not used until the next step (done by the third mode).r
      
       impulse = fatnf[deltatime, y];
        {2*eps, yinc, newyp, Impulse[{impulse, newh, newsub, deltatimesymbol, cnf, fatnf, eps}]}
    ]
];

The third mode for the impulse method takes the step past the δ zero, giving the inpulse as the solution increment.  The next step size was saved from the step size being used by the integration method before the δ zero was encounted.  The Impulse object returned is reset to normal "Step" mode.

Impulse[{impulse_?VectorQ, nexth_, submethod_, deltatimesymbol_, cnf_, fatnf_, eps_}]["Step"[rhs_, t_, h_, y_, yp_]] := 
    {nexth, impulse, None, Impulse[{"Step", submethod, deltatimesymbol, cnf, fatnf, eps}]}

And to conclude a little fun with an impulsive use of EquationTrekker.

This loads the EquationTrekker standard package.

Needs["NumericalMath`EquationTrekker`"]

This starts EquationTrekker to display solutions of a damped pendulum equation which gets a kick every so often.

\(EquationTrekker[{\(x'\)[t]\  == \ 
          y[t], \ \(y'\)[t]\  + \  .1\ y[t]\  + \ 
            Sin[x[t]]\  == \ \ ϵ\/\((x[t]^2\  + \ y[t]^2)\)\ \ DiracDelta[Cos[ω\ t^2]]}, \ {x, y}, {t, 0, 10}, \ 
      TrekGenerator -> {DifferentialEquationTrek, \ 
          Method -> {Impulse, \ Method -> "ExplicitRungeKutta"}}, \ 
      TrekParameters -> {ϵ -> {0, {0, 
                1}}, \ ω -> 1}, \ 
      PlotRange -> {{\(-5\), 5}, {\(-3\), 3}}]\)

I three examples of methods:

These should demonstrate that by following a particular template you can write particular methods or controllers for handling integration of ODEs or PDEs within the NDSolve framework.  Using the framework gives you the advantage of equation processing, step control, and solution formation all within an easy package.