"ExplicitRungeKutta" Method for NDSolve

Introduction

This loads packages containing some test problems and utility functions:

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

Euler's Method

One of the first and simplest methods for solving initial value problems was proposed by Euler:

(Subscript[y, n+1])==(Subscript[y, n])+hf((Subscript[t, n]),(Subscript[y, n])).

Euler's method is not very accurate.

Local accuracy is measured by how high terms are matched with the Taylor expansion of the solution. Euler's method is first-order accurate, so that errors occur one order higher starting at powers of (h^2).

Euler's method is implemented in (RefLink[NDSolve, paclet:ref/NDSolve]) as "ExplicitEuler":

NDSolve[{y'[t]==-y[t], y[0]==1},y[t],{t,0,1},
Method->"ExplicitEuler", "StartingStepSize"->1/10]

Generalizing Euler's Method

The idea of Runge–Kutta methods is to take successive (weighted) Euler steps to approximate a Taylor series. In this way function evaluations (and not derivatives) are used.

For example, consider the one-step formulation of the midpoint method.

({{(Subscript[k, 1]) , ==, f((Subscript[t, n]),(Subscript[y, n]))}, {(Subscript[k, 2]), ==, f((Subscript[t, n])+(1/2)h,(Subscript[y, n])+(1/2)h(Subscript[k, 1]))}, {(Subscript[y, n+1]), ==, (Subscript[y, n])+h(Subscript[k, 2])}})

The midpoint method can be shown to have a local error of O((h^3)), so it is second-order accurate.

The midpoint method is implemented in (RefLink[NDSolve, paclet:ref/NDSolve]) as "ExplicitMidpoint":

NDSolve[{y'[t]==-y[t], y[0]==1},y[t],{t,0,1},
Method->"ExplicitMidpoint", "StartingStepSize"->1/10]

Runge–Kutta Methods

Extending the approach in (), repeated function evaluation can be used to obtain higher-order methods.

Denote the Runge–Kutta method for the approximate solution to an initial value problem at (Subscript[t, n+1])==(Subscript[t, n])+h by

({{(Subscript[g, i]), ==, (Subscript[y, n])+h(∑_(j=1)^s)(Subscript[a, i,j])(Subscript[k, j]),, }, {(Subscript[k, i]), ==, f((Subscript[t, n])+(Subscript[c, i])h,(Subscript[g, i])),, i==1,2,…,s,}, {(Subscript[y, n+1]), ==, (Subscript[y, n])+h(∑_(i=1)^s)(Subscript[b, i])(Subscript[k, i]), }})

where s is the number of stages.

It is generally assumed that the row-sum conditions hold:

(Subscript[c, i])==(∑_(j=1)^s)(Subscript[a, i,j]), i == 1, …,s.

These conditions effectively determine the points in time at which the function is sampled and are a particularly useful device in the derivation of high-order Runge–Kutta methods.

The coefficients of the method are free parameters that are chosen to satisfy a Taylor series expansion through some order in the time step h. In practice other conditions such as stability can also constrain the coefficients.

Explicit Runge–Kutta methods are a special case where the matrix A is strictly lower triangular:

It has become customary to denote the method coefficients c==(([(Subscript[c, i])])^T), b==(([(Subscript[b, i])])^T), and A==[(Subscript[a, i,j])] using a Butcher table, which has the following form for explicit Runge–Kutta methods.

({{0, 0, 0, ⋯, 0, 0}, {(Subscript[c, 2]), (Subscript[a, 2,1]), 0, ⋯, 0, 0}, {⋮, ⋮, ⋮, ⋱, ⋮, ⋮}, {(Subscript[c, s]), (Subscript[a, s,1]), (Subscript[a, s,2]), ⋯, (Subscript[a, s,s-1]), 0}, {, (Subscript[b, 1]), (Subscript[b, 2]), ⋯, (Subscript[b, s-1]), (Subscript[b, s])}})

The row-sum conditions can be visualized as summing across the rows of the table.

Notice that a consequence of explicitness is (Subscript[c, 1])==0, so that the function is sampled at the beginning of the current integration step.

Example

The Butcher table for the explicit midpoint method () is given by:

({{0, 0, 0}, {(1/2), (1/2), 0}, {, 0, 1}}).

FSAL Schemes

A particularly interesting special class of explicit Runge–Kutta methods, used in most modern codes, is that for which the coefficients have a special structure known as First Same As Last (FSAL):

(Subscript[a, s,i])==(Subscript[b, i]),i==1,…,s-1and(Subscript[b, s])==0.

For consistent FSAL schemes the Butcher table () has the form:

({{0, 0, 0, ⋯, 0, 0}, {(Subscript[c, 2]), (Subscript[a, 2,1]), 0, ⋯, 0, 0}, {⋮, ⋮, ⋮, ⋱, ⋮, ⋮}, {(Subscript[c, s-1]), (Subscript[a, s-1,1]), (Subscript[a, s-1,2]), ⋯, 0, 0}, {1, (Subscript[b, 1]), (Subscript[b, 2]), ⋯, (Subscript[b, s-1]), 0}, {, (Subscript[b, 1]), (Subscript[b, 2]), ⋯, (Subscript[b, s-1]), 0}}).

The advantage of FSAL methods is that the function value (Subscript[k, s]) at the end of one integration step is the same as the first function value (Subscript[k, 1]) at the next integration step.

The function values at the beginning and end of each integration step are required anyway when constructing the (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) that is used for dense output in (RefLink[NDSolve, paclet:ref/NDSolve]).

Embedded Pairs and Local Error Estimation

An efficient means of obtaining local error estimates for adaptive step-size control is to consider two methods of different orders p and (OverHat[p]) that share the same coefficient matrix (and hence function values).

({{0, 0, 0, ⋯, 0, 0}, {(Subscript[c, 2]), (Subscript[a, 2,1]), 0, ⋯, 0, 0}, {⋮, ⋮, ⋮, ⋱, 0, ⋮}, {(Subscript[c, s-1]), (Subscript[a, s-1,1]), (Subscript[a, s-1,2]), ⋯, 0, 0}, {(Subscript[c, s]), (Subscript[a, s,1]), (Subscript[a, s,2]), ⋯, (Subscript[a, s,s-1]), 0}, {, (Subscript[b, 1]), (Subscript[b, 2]), ⋯, (Subscript[b, s-1]), (Subscript[b, s])}, {, (Subscript[(OverHat[b]), 1]), (Subscript[(OverHat[b]), 2]), ⋯, (Subscript[(OverHat[b]), s-1]), (Subscript[(OverHat[b]), s])}})

These give two solutions:

(Subscript[y, n+1])==(Subscript[y, n])+h(∑_(i=1)^s)(Subscript[b, i])(Subscript[k, i])

(Subscript[(OverHat[y]), n+1])==(Subscript[y, n])+h(∑_(i=1)^s)(Subscript[(OverHat[b]), i])(Subscript[k, i]).

A commonly used notation is p((OverHat[p])), typically with (OverHat[p])==p-1 or (OverHat[p])==p+1.

In most modern codes, including the default choice in (RefLink[NDSolve, paclet:ref/NDSolve]), the solution is advanced with the more accurate formula so that (OverHat[p])==p-1, which is known as local extrapolation.

The vector of coefficients e==(([(Subscript[b, 1])-(Subscript[(OverHat[b]), 1]),(Subscript[b, 2])-(Subscript[(OverHat[b]), 2]),…,(Subscript[b, s])-(Subscript[(OverHat[b]), s])])^T) gives an error estimator avoiding subtractive cancellation of (Subscript[y, n]) in floating-point arithmetic when forming the difference between () and ().

The quantity ‖(Subscript[err, n])‖ gives a scalar measure of the error that can be used for step size selection.

Step Control

The classical Integral (or I) step-size controller uses the formula

(Subscript[h, n+1])==(((Subscript[h, n])((Tol/(‖(Subscript[err, n])‖))))^(1/(Overscript[p, ∼])))

where (Overscript[p, ∼])==min((OverHat[p]),p)+1.

The error estimate is therefore used to determine the next step size to use from the current step size.

The notation Tol/‖(Subscript[err, n])‖ is explained within "Norms in NDSolve".

Overview

Explicit Runge–Kutta pairs of orders 2(1) through 9(8) have been implemented.

Formula pairs have the following properties:

Optimal formula pairs of orders 2(1), 3(2), and 4(3) subject to the already stated requirements have been derived using the Wolfram Language, and are described in [SS04].

The 5(4) pair selected is due to Bogacki and Shampine [BS89b, S94] and the 6(5), 7(6), 8(7), and 9(8) pairs are due to Verner [V10].

For the selection of higher-order pairs, issues such as local truncation error ratio and stability region compatibility should be considered (see [S94]). Various tools have been written to assess these qualitative features.

Methods are interchangeable so that, for example, it is possible to substitute the 5(4) method of Bogacki and Shampine with a method of Dormand and Prince.

Summation of the method stages is implemented using level 2 BLAS which is often highly optimized for particular processors and can also take advantage of multiple cores.

Example

Define the Brusselator ODE problem, which models a chemical reaction:

system=GetNDSolveProblem["BrusselatorODE"]

This solves the system using an explicit Runge–Kutta method:

sol = NDSolve[system,Method->"ExplicitRungeKutta"]

Extract the interpolating functions from the solution:

ifuns = system["DependentVariables"] /.First[sol];

Plot the solution components:

ParametricPlot[Evaluate[ifuns],Evaluate[system["TimeData"]]]

Method Comparison

Sometimes you may be interested to find out what methods are being used in (RefLink[NDSolve, paclet:ref/NDSolve]).

Here you can see the coefficients of the default 2(1) embedded pair:

NDSolve`EmbeddedExplicitRungeKuttaCoefficients[2,Infinity]

You also may want to compare some of the different methods to see how they perform for a specific problem.

Utilities

You will make use of a utility function CompareMethods for comparing various methods. Some useful (RefLink[NDSolve, paclet:ref/NDSolve]) features of this function for comparing methods are:

This displays the results of the method comparison using a (RefLink[GridBox, paclet:ref/GridBox]):

TabulateResults[labels_List,names_List,data_List]:=
DisplayForm[
FrameBox[
GridBox[
Apply[{labels,##}&,MapThread[Prepend,{data,names}]],
RowLines->True,ColumnLines->True
]
]
] /; SameQ[Length[names], Length[data]];

Reference Solution

A number of examples for comparing numerical methods in the literature rely on the fact that a closed-form solution is available, which is obviously quite limiting.

In (RefLink[NDSolve, paclet:ref/NDSolve]) it is possible to get very accurate approximations using arbitrary-precision adaptive step size; these are adaptive order methods based on "Extrapolation".

The following reference solution is computed with a method that switches between a pair of "Extrapolation" methods, depending on whether the problem appears to be stiff:

sol = NDSolve[system,Method->"StiffnessSwitching",WorkingPrecision->32];

refsol = First[FinalSolutions[system,sol]];

Automatic Order Selection

When you select "DifferenceOrder"->(RefLink[Automatic, paclet:ref/Automatic]), the code will automatically attempt to choose the optimal order method for the integration.

Two algorithms have been implemented for this purpose and are described within "SymplecticPartitionedRungeKutta Method for NDSolve".

Example 1

Here is an example that compares built-in methods of various orders, together with the method that is selected automatically.

This selects the order of the methods to choose between and makes a list of method options to pass to (RefLink[NDSolve, paclet:ref/NDSolve]):

orders = Join[Range[2,9],{Automatic}];

methods = Table[{"ExplicitRungeKutta","DifferenceOrder"->i, "StiffnessTest"->False},{i, orders}];

Compute the number of integration steps, function evaluations, and the endpoint global error:

data = CompareMethods[system, refsol, methods];

Display the results in a table:

labels = {"Method","Steps","Cost","Error"};

TabulateResults[labels, orders, data]

The default method has order nine, which is close to the optimal order of eight in this example. One function evaluation is needed during the initialization phase to determine the order.

Example 2

A limitation of the previous example is that it did not take into account the accuracy of the solution obtained by each method, so that it did not give a fair reflection of the cost.

Rather than taking a single tolerance to compare methods, it is preferable to use a range of tolerances.

The following example compares various "ExplicitRungeKutta" methods of different orders using a variety of tolerances.

This selects the order of the methods to choose between and makes a list of method options to pass to (RefLink[NDSolve, paclet:ref/NDSolve]):

orders=Join[Range[4,9],{Automatic}];

methods=Table[{"ExplicitRungeKutta","DifferenceOrder"->i,"StiffnessTest"->False},{i,orders}];

The data comparing accuracy and work is computed using CompareMethods for a range of tolerances:

tols = Range[3,14];

data=Table[Map[Rest,CompareMethods[system,refsol,methods,AccuracyGoal->i,PrecisionGoal->i]],{i,tols}];

The work-error comparison data for the various methods is displayed in the following logarithmic plot, where the global error is displayed on the vertical axis and the number of function evaluations on the horizontal axis. The default order selected (displayed in red) can be seen to increase as the tolerances are increased:

ListLogLogPlot[Transpose[data],Joined->True,Axes->False,Frame->True,PlotMarkers->Map[Style[#,Medium]&,Join[Drop[orders,-1],{"A"}]],PlotStyle->{{Black},{Black},{Black},{Black},{Black},{Black},{Red}}]

The order-selection algorithms are heuristic in that the optimal order may change through the integration but, as the examples illustrate, a reasonable default choice is usually made.

Ideally, a selection of different problems should be used for benchmarking.

Coefficient Plug-in

The implementation of "ExplicitRungeKutta" provides a default method pair at each order.

Sometimes, however, it is convenient to use a different method, for example:

The Classical Runge–Kutta Method

This shows how to define the coefficients of the classical explicit Runge–Kutta method of order four, approximated to precision p:

crkamat={{1/2},{0,1/2},{0,0,1}};
crkbvec={1/6,1/3,1/3,1/6};
crkcvec={1/2,1/2,1};ClassicalRungeKuttaCoefficients[4,p_]:=N[{crkamat,crkbvec,crkcvec},p];

The method has no embedded error estimate and hence there is no specification of the coefficient error vector. This means that the method is invoked with fixed step sizes.

Here is an example of the calling syntax:

NDSolve[system,Method->{"ExplicitRungeKutta","DifferenceOrder"->4,"Coefficients"->ClassicalRungeKuttaCoefficients},StartingStepSize->1/10]

ode23

This defines the coefficients for a 3(2) FSAL explicit Runge–Kutta pair.

The third-order formula is due to Ralston, and the embedded method was derived by Bogacki and Shampine [BS89a].

This defines a function for computing the coefficients to a desired precision:

BSamat={{1/2},{0,3/4},{2/9,1/3,4/9}};
BSbvec={2/9,1/3,4/9,0};
BScvec = {1/2,3/4,1};
BSevec={-5/72,1/12,1/9,-1/8};
BSCoefficients[4, p_]:=
N[{BSamat,BSbvec,BScvec,BSevec}, p];

The method is used in the Texas Instruments TI-85 pocket calculator, Matlab, and RKSUITE [S94]. Unfortunately it does not allow for the form of stiffness detection that has been chosen.

A Method of Fehlberg

This defines the coefficients for a 4(5) explicit Runge–Kutta pair of Fehlberg that was popular in the 1960s [F69].

The fourth-order formula is used to propagate the solution, and the fifth-order formula is used only for the purpose of error estimation.

This defines the function for computing the coefficients to a desired precision:

Fehlbergamat={
{1/4},
{3/32,9/32},
{1932/2197,-7200/2197,7296/2197},{439/216,-8,3680/513,-845/4104},{-8/27,2,-3544/2565,1859/4104,-11/40}};
Fehlbergbvec={25/216,0,1408/2565,2197/4104,-1/5, 0};
Fehlbergcvec = {1/4,3/8,12/13,1,1/2};
Fehlbergevec={-1/360,0,128/4275,2197/75240,-1/50,-2/55};
FehlbergCoefficients[4, p_]:=
N[{Fehlbergamat,Fehlbergbvec,Fehlbergcvec,Fehlbergevec}, p];

In contrast to the classical Runge–Kutta method of order four, the coefficients include an additional entry that is used for error estimation.

The Fehlberg method is not a FSAL scheme since the coefficient matrix is not of the form (); it is a six-stage scheme, but it requires six function evaluations per step because of the function evaluation that is required at the end of the step to construct the (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]).

A Dormand–Prince Method

Here is how to define a 5(4) pair of Dormand and Prince coefficients [DP80]. This is currently the method used by ode45 in Matlab.

This defines a function for computing the coefficients to a desired precision:

DOPRIamat={
{1/5},
{3/40,9/40},
{44/45,-56/15,32/9},
{19372/6561,-25360/2187,64448/6561,-212/729},
{9017/3168,-355/33,46732/5247,49/176,-5103/18656},
{35/384,0,500/1113,125/192,-2187/6784,11/84}};
DOPRIbvec={35/384,0,500/1113,125/192,-2187/6784,11/84, 0};
DOPRIcvec = {1/5,3/10,4/5,8/9,1, 1};
DOPRIevec={71/57600,0,-71/16695,71/1920,-17253/339200,22/525,-1/40};
DOPRICoefficients[5, p_]:=
N[{DOPRIamat,DOPRIbvec,DOPRIcvec,DOPRIevec}, p];

The Dormand–Prince method is a FSAL scheme since the coefficient matrix is of the form (); it is a seven-stage scheme, but effectively uses only six function evaluations.

Here is how the coefficients of Dormand and Prince can be used in place of the built-in choice. Since the structure of the coefficients includes an error vector, the implementation is able to ascertain that adaptive step sizes can be computed:

NDSolve[system,Method->{"ExplicitRungeKutta","DifferenceOrder"->5,"Coefficients"->DOPRICoefficients, "StiffnessTest"->False}]

Method Comparison

Here you solve a system using several explicit Runge–Kutta pairs.

For the Fehlberg 4(5) pair, the option "EmbeddedDifferenceOrder" is used to specify the order of the embedded method:

Fehlberg45 = {"ExplicitRungeKutta","Coefficients"->FehlbergCoefficients,"DifferenceOrder"->4,"EmbeddedDifferenceOrder"->5, "StiffnessTest"->False};

The Dormand and Prince 5(4) pair is defined as follows:

DOPRI54 ={"ExplicitRungeKutta","Coefficients"->DOPRICoefficients,"DifferenceOrder"->5, "StiffnessTest"->False};

The 5(4) pair of Bogacki and Shampine is the default order-five method:

BS54 ={"ExplicitRungeKutta","Coefficients"->"EmbeddedExplicitRungeKuttaCoefficients","DifferenceOrder"->5, "StiffnessTest"->False};

Put the methods and some descriptive names together in a list:

names = {"Fehlberg 4(5)", "Dormand-Prince 5(4)","Bogacki-Shampine 5(4)"};

methods = {Fehlberg45,DOPRI54,BS54};

Compute the number of integration steps, function evaluations, and the endpoint global error:

data = CompareMethods[system,refsol, methods];

Display the results in a table:

labels = {"Method","Steps","Cost","Error"};
TabulateResults[labels, names, data]

The default method was the least expensive and provided the most accurate solution.

Method Plug-in

This shows how to implement the classical explicit Runge–Kutta method of order four using the method plug-in environment.

This definition is optional since the method in fact has no data. However, any expression can be stored inside the data object. For example, the coefficients could be approximated here to avoid coercion from rational to floating-point numbers at each integration step:

ClassicalRungeKutta/:NDSolve`InitializeMethod[ClassicalRungeKutta,__]:=ClassicalRungeKutta[];

The actual method implementation is written using a stepping procedure:

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}
];

Notice that the implementation closely resembles the description that you might find in a textbook. There are no memory allocation/deallocation statements or type declarations, for example. In fact the implementation works for machine real numbers or machine complex numbers, and even using arbitrary-precision software arithmetic.

Here is an example of the calling syntax. For simplicity the method only uses fixed step sizes, so you need to specify what step sizes to take:

NDSolve[system,Method->ClassicalRungeKutta,StartingStepSize->1/10]

Many of the methods that have been built into (RefLink[NDSolve, paclet:ref/NDSolve]) were first prototyped using top-level code before being implemented in the kernel for efficiency.

Stiffness

Stiffness is a combination of problem, initial data, numerical method, and error tolerances.

Stiffness can arise, for example, in the translation of diffusion terms by divided differences into a large system of ODEs.

In order to understand more about the nature of stiffness it is useful to study how methods behave when applied to a simple problem.

Linear Stability

Consider applying a Runge–Kutta method to a linear scalar equation known as Dahlquist's equation:

y'(t)==λ y(t),λ∈ℂ,Re(λ)<0.

The result is a rational function of polynomials R(z) where z==h λ (see for example [L87]).

This utility function finds the linear stability function R(z) for Runge–Kutta methods. The form depends on the coefficients and is a polynomial if the Runge–Kutta method is explicit.

Here is the stability function for the fifth-order scheme in the Dormand–Prince 5(4) pair:

DOPRIsf = RungeKuttaLinearStabilityFunction[DOPRIamat, DOPRIbvec,z]

This function finds the linear stability function R(z) for Runge–Kutta methods. The form depends on the coefficients and is a polynomial if the Runge–Kutta method is explicit.

The following package is useful for visualizing linear stability regions for numerical methods for differential equations:

Needs["FunctionApproximations`"];

You can now visualize the absolute stability region |R(z)|==1:

OrderStarPlot[DOPRIsf, 1,z]

Depending on the magnitude of λ in (), if you choose the step size h such that |R(h λ)|<1, then errors in successive steps will be damped, and the method is said to be absolutely stable.

If \[LeftBracketingBar]R(h λ)\[RightBracketingBar]>1, then step-size selection will be restricted by stability and not by local accuracy.

Stiffness Detection

The device for stiffness detection that is used with the option "StiffnessTest" is described within "StiffnessTest Method Option for NDSolve".

Recast in terms of explicit Runge–Kutta methods, the condition for stiffness detection can be formulated as

(Overscript[λ, ∼])==((‖(Subscript[k, s])-(Subscript[k, s-1])‖)/(‖(Subscript[g, s])-(Subscript[g, s-1])‖))

with (Subscript[g, i]) and (Subscript[k, i]) defined in ().

The difference (Subscript[g, s])-(Subscript[g, s-1]) can be shown to correspond to a number of applications of the power method applied to hJ.

The difference is therefore a good approximation of the eigenvector corresponding to the leading eigenvalue.

The product |h(Overscript[λ, ∼])| gives an estimate that can be compared to the stability boundary in order to detect stiffness.

An s-stage explicit Runge–Kutta has a form suitable for () if (Subscript[c, s-1])==(Subscript[c, s])==1.

({{0, 0, 0, ⋯, 0, 0}, {(Subscript[c, 2]), (Subscript[a, 2,1]), 0, ⋯, 0, 0}, {⋮, ⋮, ⋮, ⋱, ⋮, ⋮}, {1, (Subscript[a, s-1,1]), (Subscript[a, s-1,2]), ⋯, 0, 0}, {1, (Subscript[a, s,1]), (Subscript[a, s,2]), ⋯, (Subscript[a, s,s-1]), 0}, {, (Subscript[b, 1]), (Subscript[b, 2]), ⋯, (Subscript[b, s-1]), (Subscript[b, s])}})

The default embedded pairs used in "ExplicitRungeKutta" all have the form ().

An important point is that () is very cheap and convenient; it uses already available information from the integration and requires no additional function evaluations.

Another advantage of () is that it is straightforward to make use of consistent FSAL methods ().

Examples

Select a stiff system modeling a chemical reaction:

system = GetNDSolveProblem["Robertson"];

This applies a built-in explicit Runge–Kutta method to the stiff system.

By default stiffness detection is enabled, since it only has a small impact on the running time:

NDSolve[system, Method->"ExplicitRungeKutta"];

The coefficients of the Dormand–Prince 5(4) pair are of the form () so stiffness detection is enabled:

NDSolve[system, Method->{"ExplicitRungeKutta","DifferenceOrder"->5,"Coefficients"->DOPRICoefficients}];

Since no "LinearStabilityBoundary" property has been specified, a default value is chosen. In this case the value corresponds to a generic method of order 5:

genlsb = NDSolve`LinearStabilityBoundary[5]

You can set up an equation in terms of the linear stability function and solve it exactly to find the point where the contour crosses the negative real axis:

DOPRIlsb = Reduce[Abs[DOPRIsf] == 1 && z < 0,z]

The default generic value is very slightly smaller in magnitude than the computed value:

N[{genlsb, DOPRIlsb[[2]]}]

In general, there may be more than one point of intersection, and it may be necessary to choose the appropriate solution.

The following definition sets the value of the linear stability boundary:

DOPRICoefficients[5]["LinearStabilityBoundary"]=Root[600+300*#1+100*#1^2+25*#1^3+5*#1^4+#1^5&,1,0];

Using the new value for this example does not affect the time at which stiffness is detected:

NDSolve[system, Method->{"ExplicitRungeKutta","DifferenceOrder"->5,"Coefficients"->DOPRICoefficients}];

The Fehlberg 4(5) method does not have the correct coefficient structure () required for stiffness detection, since (Subscript[c, s])==1/2≠1.

The default value "StiffnessTest"->(RefLink[Automatic, paclet:ref/Automatic]) checks to see if the method coefficients provide a stiffness detection capability; if they do, then stiffness detection is enabled.

Step Control Revisited

There are some reasons to look at alternatives to the standard Integral step controller () when considering mildly stiff problems.

This system models a chemical reaction:

system = GetNDSolveProblem["Robertson"];

This defines an explicit Runge–Kutta method based on the Dormand–Prince coefficients that does not use stiffness detection:

IERK = {"ExplicitRungeKutta","Coefficients"->DOPRICoefficients,"DifferenceOrder"->5,"StiffnessTest"->False};

This solves the system and plots the step sizes that are taken using the utility function StepDataPlot:

isol = NDSolve[system, Method->IERK];
StepDataPlot[isol]

Solving a stiff or mildly stiff problem with the standard step-size controller leads to large oscillations, sometimes leading to a number of undesirable step-size rejections.

The study of this issue is known as step-control stability.

It can be studied by matching the linear stability regions for the high- and low-order methods in an embedded pair.

One approach to addressing the oscillation is to derive special methods, but this compromises the local accuracy.

PI Step Control

An appealing alternative to Integral step control () is Proportional-Integral or PI step control.

In this case the step size is selected using the local error in two successive integration steps according to the formula:

(Subscript[h, n+1])==(((Subscript[h, n])((Tol/(‖(Subscript[err, n])‖))))^((Subscript[k, 1])/(Overscript[p, ∼])))(((((‖(Subscript[err, n-1])‖)/(‖(Subscript[err, n])‖))))^((Subscript[k, 2])/(Overscript[p, ∼]))).

This has the effect of damping and hence gives a smoother step-size sequence.

Note that Integral step control () is a special case of () and is used if a step is rejected:

The option "StepSizeControlParameters"->{(Subscript[k, 1]),(Subscript[k, 2])} can be used to specify the values of (Subscript[k, 1]) and (Subscript[k, 2]).

The scaled error estimate in () is taken to be ‖(Subscript[err, n-1])‖==‖(Subscript[err, n])‖ for the first integration step.

Examples

Stiff Problem

This defines a method similar to IERK that uses the option "StepSizeControlParameters" to specify a PI controller.

Here you use generic control parameters suggested by Gustafsson:

This specifies the step-control parameters:

PIERK ={"ExplicitRungeKutta","Coefficients"->DOPRICoefficients,"DifferenceOrder"->5,"StiffnessTest"->False,"StepSizeControlParameters"->{3/10, 2/5}};

Solving the system again, it can be observed that the step-size sequence is now much smoother:

pisol = NDSolve[system, Method->PIERK];
StepDataPlot[pisol]

Nonstiff Problem

In general the I step controller () is able to take larger steps for a nonstiff problem than the PI step controller () as the following example illustrates.

Select and solve a nonstiff system using the I step controller:

system = GetNDSolveProblem["BrusselatorODE"];

isol = NDSolve[system, Method->IERK];
StepDataPlot[isol]

Using the PI step controller the step sizes are slightly smaller:

pisol = NDSolve[system, Method->PIERK];
StepDataPlot[pisol]

For this reason, the default setting for "StepSizeControlParameters" is (RefLink[Automatic, paclet:ref/Automatic]) , which is interpreted as:

Fine-Tuning

Instead of using () directly, it is common practice to use safety factors to ensure that the error is acceptable at the next step with high probability, thereby preventing unwanted step rejections.

The option "StepSizeSafetyFactors"->{(Subscript[s, 1]),(Subscript[s, 2])} specifies the safety factors to use in the step-size estimate so that () becomes:

(Subscript[h, n+1])==(Subscript[h, n])(Subscript[s, 1])((((((Subscript[s, 2])Tol)/(‖(Subscript[err, n])‖))))^((Subscript[k, 1])/(Overscript[p, ∼])))(((((‖(Subscript[err, n-1])‖)/(‖(Subscript[err, n])‖))))^((Subscript[k, 2])/(Overscript[p, ∼]))).

Here (Subscript[s, 1]) is an absolute factor and (Subscript[s, 2]) typically scales with the order of the method.

The option "StepSizeRatioBounds"->{(Subscript[sr, min]),(Subscript[sr, max])} specifies bounds on the next step size to take such that:

(Subscript[sr, min])≤|((Subscript[h, n+1])/(Subscript[h, n]))|≤(Subscript[sr, max]).

Option Summary

Options of the method "ExplicitRungeKutta".

The default setting of (RefLink[Automatic, paclet:ref/Automatic]) for the option "DifferenceOrder" selects the default coefficient order based on the problem, initial values, and local error tolerances, balanced against the work of the method for each coefficient set.

The default setting of (RefLink[Automatic, paclet:ref/Automatic]) for the option "EmbeddedDifferenceOrder" specifies that the default order of the embedded method is one lower than the method order. This depends on the value of the "DifferenceOrder" option.

The default setting of (RefLink[Automatic, paclet:ref/Automatic]) for the option "StepSizeControlParameters" uses the values {1,0} if stiffness detection is active and {3/10,2/5} otherwise.

The default setting of (RefLink[Automatic, paclet:ref/Automatic]) for the option "StepSizeSafetyFactors" uses the values {17/20,9/10} if the I step controller () is used and {9/10,9/10} if the PI step controller () is used. The step controller used depends on the values of the options "StepSizeControlParameters" and "StiffnessTest".

The default setting of (RefLink[Automatic, paclet:ref/Automatic]) for the option "StiffnessTest" will activate the stiffness test if the coefficients have the form ().