"DoubleStep" Method for NDSolve

Introduction

The method "DoubleStep" performs a single application of Richardson's extrapolation for any one-step integration method.

Although it is not always optimal, it is a general scheme for equipping a method with an error estimate (hence adaptivity in the step size) and extrapolating to increase the order of local accuracy.

"DoubleStep" is a special case of extrapolation but has been implemented as a separate method for efficiency.

Given a method of order p:

e==(( (Subscript[y, 2])- (Subscript[y, 1]))/((2^p)- 1)).

(Subscript[(OverHat[y]), 2])==(Subscript[y, 2 ])+e.

Examples

Load some package with example problems and utility functions:

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

Select a nonstiff problem from the package:

nonstiffsystem = GetNDSolveProblem["BrusselatorODE"];

Select a stiff problem from the package:

stiffsystem = GetNDSolveProblem["Robertson"];

Extending Built-in Methods

The method "ExplicitEuler" carries out one integration step using Euler's method. It has no local error control and hence uses fixed step sizes.

This integrates a differential system using one application of Richardson's extrapolation (see ()) with the base method "ExplicitEuler".

The local error estimate () is used to dynamically adjust the step size throughout the integration:

eesol = NDSolve[nonstiffsystem,{T, 0, 1},Method->{"DoubleStep", Method->"ExplicitEuler"}]

This illustrates how the step size varies during the numerical integration:

StepDataPlot[eesol]

The stiffness detection device (described within "StiffnessTest Method Option for NDSolve") ascertains that the "ExplicitEuler" method is restricted by stability rather than local accuracy:

NDSolve[stiffsystem,Method->{"DoubleStep", Method->"ExplicitEuler"}]

An alternative base method is more appropriate for this problem:

liesol = NDSolve[stiffsystem,Method->{"DoubleStep", Method->"LinearlyImplicitEuler"}]

User-Defined Methods and Method Properties

Integration methods can be added to the (RefLink[NDSolve, paclet:ref/NDSolve]) framework.

In order for these to work like built-in methods it can be necessary to specify various method properties. These properties can then be used by other methods to build up compound integrators.

Here is how to define a top-level plug-in for the classical Runge–Kutta method (see "NDSolve Method Plug-in Framework: Classical Runge–Kutta" and "ExplicitRungeKutta Method for NDSolve" for more details):

ClassicalRungeKutta[___]["Step"[f_,t_,h_,y_,yp_]]:=
Block[{deltay,k1,k2,k3,k4},
k1=yp;
k2=f[t+1/2 h,y+1/2 h k1];
k3=f[t+1/2 h,y+1/2 h k2];
k4=f[t+h,y+h k3];
deltay=h (1/6 k1+1/3 k2+1/3 k3+1/6 k4);
{h,deltay}
];

Method properties used by "DoubleStep" are now described.

Order and Symmetry

This attempts to integrate a system using one application of Richardson's extrapolation based on the classical Runge–Kutta method:

NDSolve[nonstiffsystem,Method->{"DoubleStep",Method->ClassicalRungeKutta}];

Without knowing the order of the base method, "DoubleStep" is unable to carry out Richardson's extrapolation.

This defines a method property to communicate to the framework that the classical Runge–Kutta method has order four:

ClassicalRungeKutta[___]["DifferenceOrder"]:=4;

The method "DoubleStep" is now able to ascertain that ClassicalRungeKutta is of order four and can use this information when refining the solution and estimating the local error:

NDSolve[nonstiffsystem,Method->{"DoubleStep",Method->ClassicalRungeKutta}]

The order of the result of Richardson's extrapolation depends on whether the extrapolated method has a local error expansion in powers of h or (h^2) (the latter occurs if the base method is symmetric).

If no method property for symmetry is defined, the "DoubleStep" method assumes by default that the base integrator is not symmetric.

This explicitly specifies that the classical Runge–Kutta method is not symmetric using the "SymmetricMethodQ" property:

ClassicalRungeKutta[___]["SymmetricMethodQ"]:=False;

Stiffness Detection

Details of the scheme used for stiffness detection can be found within "StiffnessTest Method Option for NDSolve".

Stiffness detection relies on knowledge of the linear stability boundary of the method, which has not been defined.

Computing the exact linear stability boundary of a method under extrapolation can be quite complicated. Therefore a default value is selected which works for all methods. This corresponds to considering the p(^th) order power series approximation to the exponential at 0 and ignoring higher order terms.

This computes the linear stability boundary for a generic method of order 4:

Reduce[Abs[Sum[((z^i)/(i!)),{i,0,4}]]==1 &&z<0,z]

A default value for the "LinearStabilityBoundary" property is used:

NDSolve[stiffsystem,Method->{"DoubleStep",Method->ClassicalRungeKutta,"StiffnessTest"->True}];

This shows how to specify the linear stability boundary of the method for the framework. This value will only be used if "DoubleStep" is invoked with "LocalExtrapolation"->True :

ClassicalRungeKutta[___]["LinearStabilityBoundary"]:=Root[24+12 #1+4 ((#1)^2)+((#1)^3)&,1];

"DoubleStep" assumes by default that a method is not appropriate for stiff problems (and hence uses stiffness detection) when no "StiffMethodQ" property is specified. This shows how to define the property:

ClassicalRungeKutta[___]["StiffMethodQ"]:=False;

Higher Order

The following example extrapolates the classical Runge–Kutta method of order four using two applications of ().

The inner specification of "DoubleStep" constructs a method of order five.

A second application of "DoubleStep" is used to obtain a method of order six, which uses adaptive step sizes.

Nested applications of "DoubleStep" are used to raise the order and provide an adaptive step-size estimate:

NDSolve[nonstiffsystem,Method->{"DoubleStep", Method->{"DoubleStep",Method->ClassicalRungeKutta}}]

In general the method "Extrapolation" is more appropriate for constructing high-order integration schemes from low-order methods.

Option Summary

Options of the method "DoubleStep".

The default setting of (RefLink[Automatic, paclet:ref/Automatic]) for the option "StiffnessTest" indicates that the stiffness test is activated if a nonstiff base method is used.

The default setting of (RefLink[Automatic, paclet:ref/Automatic]) for the option "StepSizeSafetyFactors" uses the values {9/10,4/5} for a stiff base method and {9/10,13/20} for a nonstiff base method.