ExplicitRungeKutta Method for NDSolve
Introduction
Extrapolation methods are a class of arbitrary-order methods with automatic order and step size control. The error estimate comes from computing a solution over an interval using the same method with a varying number of steps and using extrapolation on the polynomial that fits through the computed solutions, giving a composite higher-order method [
BS64]. At the same time, the polynomials give a means of error estimation.
Typically, for low precision, the extrapolation methods have not been competitive with Runge-Kutta-type methods. For high precision, however, the arbitrary order means that they can be arbitrarily faster than fixed-order methods for very precise tolerances.
The order and step size control is based on the codes
odex.f and
seulex.f described in [
HNW93] and [
HW96].
This loads packages that contain some utility functions for plotting step sequences and some predefined problems. |
Extrapolation
The method
DoubleStep performs a single application of Richardson's extrapolation for any one-step integration method and is described within
"DoubleStep Method for NDSolve".
Extrapolation generalizes the idea of Richardson's extrapolation to a sequence of refinements. A solution is computed over an interval using the same method with a varying number of steps and using extrapolation on the polynomial that fits through the computed solutions, giving a composite higher-order method. The polynomials at the same time provide a means of error estimation.
Consider a differential system
Let
H>0 be a basic step size and choose a monotonically increasing sequence of positive integers
and define the corresponding step sizes
Choose a numerical method of order
p and compute the solution of the initial value problem by carrying out
ni steps with step size
hi to obtain:
Extrapolation is performed using the Aitken-Neville algorithm by building up a table of values:
A dependency graph of the values illustrates the relationship:
For
k=2,
n1=1,
n2=2 extrapolation is identical to Richardson's extrapolation.
Extrapolation sequences
Any extrapolation sequence can be specified in the implementation. Some common choices are as follows.
This is the Romberg sequence.
| Out[5]= |  |
|
This is the Bulirsch sequence.
| Out[6]= |  |
|
This is the Harmonic sequence.
| Out[7]= |  |
|
A sequence that satisfies
(ni/ni-j+1)p≥2 has the effect of minimizing the roundoff errors for an order-
p base integration method.
For a base method of order two the first entries in the sequence are given by the following.
| Out[8]= |  |
|
Here is an example of adding a function to define the Harmonic sequence is where the method order is an optional pattern. |
Rounding error accumulation
For high-order extrapolation an important consideration is the accumulation of rounding errors in the Aitken-Neville algorithm (
1).
As an example consider Exercise 5 of Section II.9 in [
HNW93].
Suppose that the entries
T11, T21, T31, ... are disturbed with rounding errors
, -, , ... and compute the propagation of these errors into the extrapolation table.
Due to the linearity of the extrapolation process (
2), suppose that the
Ti, j are equal to zero and take
=1.
This shows the evolution of the Aitken-Neville algorithm (
3) on the initial data using the Harmonic sequence and a symmetric order-two base integration method.
Hence, for an order-sixteen method approximately two decimal digits are lost due to rounding error accumulation.
This model is somewhat crude because, as you will see later, it is more likely that rounding errors are made in
Ti+1, 1 than in
Ti, 1 for
i≥1.
Rounding error reduction
It seems worthwhile to look for approaches that can reduce the effect of rounding errors in high-order extrapolation.
Selecting a different step sequence to diminish rounding errors is one approach, although the drawback is that the number of integration steps needed to form the
Ti, 1 in the first column of the extrapolation table requires more work.
An alternative strategy, which does not appear to have received a great deal of attention in the context of extrapolation, is to modify the base-integration method in order to reduce the magnitude of the rounding errors in floating-point operations. This approach, based on ideas that dated back to [
G51] and used to good effect for the two body problem in [
F96b] (for background see also [
K65], [
M65a], [
M65b], [
V79]), are explained next.
Base methods
The following methods are the most common choices for base integrators in extrapolation.
- ExplicitModifiedMidpoint (Gragg smoothing step)
- LinearlyImplicitMidpoint (Bader-Deuflhard formulation without smoothing step)
- LinearlyImplicitModifiedMidpoint (Bader-Deuflhard formulation with smoothing step)
For efficiency, these have been built into
NDSolve and are callable via the
Method option.
The implementation of these methods has a special interpretation for multiple substeps within
DoubleStep and
Extrapolation.
The methods are now described together with an increment reformulation that is used to reduce rounding error accumulation throughout
NDSolve.
Multiple Euler steps
Given
t0,
y0 and
H, consider a succession of
n integration steps with step size
h=H/n carried out using Euler's method:
Correspondence with explicit Runge-Kutta methods
It is well known that, for certain base integration schemes, the entries
Ti, j in the extrapolation table produced from (
4) correspond to explicit Runge-Kutta methods (see Exercise 1, Section II.9 in [
HNW93]).
For example, (
5) is equivalent to an
n-stage explicit Runge-Kutta method:
where the coefficients are represented by the Butcher table:
Reformulation
Let
yn=yn+1-yn. Then the integration (
6) can be rewritten to reflect the correspondence with an explicit Runge-Kutta method (
7,
8) as:
where terms in the right-hand side of (
9) are now considered as departures from the
same value y0.
The
yi in (
10) correspond to the
ki in (
11).
Mathematically the formulations (
12) and (
13,
14) are equivalent. For
n>1, however, the computations in (
15) have the advantage of accumulating a sum of smaller
O (h) quantities, or increments, which reduces rounding error accumulation in finite-precision floating-point arithmetic.
Multiple linearly implicit Euler steps
Next consider a succession of integration steps carried out using the linearly implicit Euler method. Increments arise naturally in the description of many semi-implicit and implicit methods.
Here
I denotes the Identity matrix and
J denotes the Jacobian of
f:
The solution of the equations for the increments in (
16) is accomplished using a single LU decomposition of the matrix
I-h J followed by the solution of triangular linear systems for each right-hand side.
The desired result is obtained from (
17) as:
Reformulation
Reformulation in terms of increments as departures from
y0 can be accomplished as follows:
The result for
yn using (
18) is obtained from (
19).
Notice that (
20) and (
21) are equivalent when
J=0.
Multiple explicit midpoint steps
"Expansions in even powers of
h are extremely important for an efficient implementation of Richardson's extrapolation" [
S70].
Consider a succession of integration steps carried out using one Euler step followed by multiple explicit midpoint steps:
If (
22) is computed with
2m-1 midpoint steps, then the method has a symmetric error expansion ([
G65], [
S70]).
Reformulation
Reformulation of (
23) can be accomplished in terms of increments as:
Gragg's smoothing step
The
smoothing step of Gragg has its historical origins in the weak stability of the explicit midpoint rule:
In order to make use of (
24), the formulation (
25) is computed with
2m steps. This has the advantage of increasing the stability domain and evaluating the function at the end of the basic step [
HNW93].
Notice that because of the construction, a sum of increments is available at the end of the algorithm together with two consecutive increments. This leads to the following formulation:
Moreover (
26) has an advantage over (
27) in finite-precision arithmetic because the values
yi, which typically have a larger magnitude than the increments
yi, do not contribute to the computation.
Gragg's smoothing step is not of great importance if the method is followed by extrapolation, and Shampine proposes an alternative smoothing procedure that is slightly more efficient [
SB83].
The method
ExplicitMidpoint uses
2m-1 steps and
ExplicitModifiedMidpoint uses
2m steps followed by the smoothing step (
28).
Stability regions
This illustrates the effect of the smoothing step on the linear stability domain (carried out using the package
FunctionApproximations.m).
Multiple linearly implicit midpoint steps
Consider one step of the linearly implicit Euler method followed by multiple linearly implicit midpoint steps, using the formulation of Bader and Deuflhard [
BD83]:
If (
29) is computed for
2m-1 linearly implicit midpoint steps, then the method has a symmetric error expansion [
BD83].
Reformulation
Reformulation of (
30) in terms of increments can be accomplished as follows:
Smoothing step
An appropriate smoothing step for the linearly implicit midpoint rule is [
BD83]:
Bader's smoothing step (
31) rewritten in terms of increments becomes:
The required quantities are obtained when (
32) is run with
2m steps.
The smoothing step for the linearly implicit midpoint rule has a different role to Gragg's smoothing for the explicit midpoint rule (see [
BD83] and [
SB83]). Since there is no weakly stable term to eliminate, the aim is to improve the asymptotic stability.
The method
LinearlyImplicitMidpoint uses
2m-1 steps and
LinearlyImplicitModifiedMidpoint uses
2m steps followed by the smoothing step (
33).
Polynomial extrapolation in terms of increments
You have seen how to modify
Ti, 1, the entries in the first column of the extrapolation table, in terms of increments.
However, for certain base integration methods,
each of the
Ti, j corresponds to an explicit Runge-Kutta method.
Therefore, it appears that the correspondence has not yet been fully exploited and further refinement is possible.
Since the Aitken-Neville algorithm (
34) involves linear differences, the entire extrapolation process can be carried out using increments.
This leads to the following modification of the Aitken-Neville algorithm:
The quantities
Ti, j=Ti, j-y0 in (
35) can be computed iteratively, starting from the initial quantities
Ti, 1 that are obtained from the modified base integration schemes without adding the contribution from
y0.
The final desired value
Tk, k can be recovered as
Tk, k+y0.
The advantage is that the extrapolation table is built up using smaller quantities, and so the effect of rounding errors from subtractive cancellation is reduced.
Implementation issues
There are a number of important implementation issues that should be considered, some of which are mentioned here.
Jacobian reuse
The Jacobian is evaluated only once for all entries
Ti, 1 at each time step by storing it together with the associated time that it is evaluated. This also has the advantage that the Jacobian does not need to be recomputed for rejected steps.
Dense linear algebra
For dense systems, the LAPACK routines xyytrf can be used for the LU decomposition and the routines xyytrs for solving the resulting triangular systems [
LAPACK99].
Adaptive order and work estimation
In order to adaptively change the order of the extrapolation throughout the integration, it is important to have a measure of the amount of work required by the base scheme and extrapolation sequence.
A measure of the relative cost of function evaluations is advantageous.
The dimension of the system, preferably with a weighting according to structure, needs to be incorporated for linearly implicit schemes in order to take account of the expense of solving each linear system.
Stability check
Extrapolation methods use a large basic step size that can give rise to some difficulties.
"Neither code can solve the van der Pol equation problem in a straightforward way because of overflow..." [
S87].
Two forms of stability check are used for the linearly implicit base schemes (for further discussion, see [
HW96]).
One check is performed during the extrapolation process. Let
errj=
Tj, j-1-Tj, j
.
In order to interrupt computations in the computation of
T1, 1, Deuflhard suggests checking if the Newton iteration applied to a fully implicit scheme would converge.
For the implicit Euler method this leads to consideration of:
Notice that (
36) differs from (
37) only in the second equation. It requires finding the solution for a different right-hand side but no extra function evaluation.
Increments are a more accurate formulation for the implementation of both forms of stability check.
Examples
Stiff systems
One of the simplest nonlinear equations describing a circuit is van der Pol's equation. |
This solves the equations using Extrapolation with the LinearlyImplicitEuler base method with the default sub-Harmonic sequence 2, 3, 4, ....
| Out[14]= |  |
|
Notice that the Jacobian matrix is computed automatically (user specifiable by using either numerical differences or symbolic derivatives) and appropriate linear algebra routines are selected and invoked at run time.
This plots the first solution component over time.
| Out[15]= |  |
|
This plots the step sizes taken in computing the solution.
| Out[16]= |  |
|
High-precision comparison
Select the Lorenz equations. |
This invokes a bigfloat, or software floating-point number, embedded explicit Runge-Kutta method of order 9(8) [ V78].
| Out[18]= |  |
|
This invokes the Adams method using a bigfloat version of LSODA. The maximum order of these methods is twelve.
| Out[19]= |  |
|
This invokes the Extrapolation method with ExplicitModifiedMidpoint as the base integration scheme.
| Out[20]= |  |
|
As many as 14 entries in the first column of the extrapolation table are used, which results in a method of order 28. This means that a much larger step size can be taken.
Here are the step sizes taken by the various methods.
| Out[23]= |  |
|
Here are the solutions obtained by the various methods at the end of the integration.
Out[24]//DisplayForm= |
| |  |
|
For comparison, here is a higher-precision reference solution obtained using Extrapolation.
| Out[26]= |  |
|
Work-error comparison
For comparing different extrapolation schemes consider an example from [ HW96]. |
The exact solution is given by
y (t)=1/cos (t).
Increment formulation
This example involves an eighth-order extrapolation of
ExplicitEuler with the Harmonic sequence. Approximately two digits of accuracy are gained by using the increment-based formulation throughout the extrapolation process.
- The results for the standard formulation (38) are depicted in green.
- The results for the increment formulation (39) followed by standard extrapolation (40) is depicted in blue.
- The results for the increment formulation (41) with extrapolation carried out on the increments using (42) is depicted in red.
Approximately two decimal digits of accuracy are gained by using the increment-based formulation throughout the extrapolation process.
This compares the relative error in the integration data that forms the initial column of the extrapolation table for the previous example.
Reference values were computed using software arithmetic with 32 decimal digits and converted to the nearest IEEE double-precision floating-point numbers, where an ULP signifies a Unit in the Last Place or Unit in the Last Position.
Notice that the rounding error model that was used to motivate the study of rounding error growth is limited because in practice errors in
Ti, 1 can exceed 1 ULP.
The increment formulation used throughout the extrapolation process produces rounding errors in
Ti, 1 that are smaller than 1 ULP.
Method comparison
This compares the work required for extrapolation based on
ExplicitEuler (red), the
ExplicitMidpoint (blue), and
ExplicitModifiedMidpoint (green).
All computations are carried out using software arithmetic with 32 decimal digits.
Fine-tuning
StepSizeSafetyFactors
As with most methods, you want a balance between taking too small a step and trying to take too big a step that will be frequently rejected. The option
StepSizeSafetyFactors->{s1, s2} constrains the choice of step size as follows. The step size chosen by the method for order
p satisfies:
This includes both an order-dependent factor and an order-independent factor.
StepSizeRatioBounds
The option
StepSizeRatioBounds->{srmin, srmax}specifies bounds on the next step size to take such that:
OrderSafetyFactors
An important aspect in
Extrapolation is the choice of order.
Each extrapolation of order
k has an associated work estimate
k which is based on, for example, the number of function evaluations of the base method and the step sequence used.
The work estimate for explicit base methods is based on the number of function evaluations and the step sequence used.
The work estimate for linearly implicit base methods also includes an estimate of the cost of evaluating the Jacobian, the cost of an LU decomposition, and the cost of back solving the linear equations.
Estimates for the work per unit step are formed from the work estimate
k and the expected new step size to take for a method of order
k (computed from (
43)):
.
Comparing consecutive estimates,
k allows a decision about when a different order method will be more efficient.
The option
OrderSafetyFactors->{f1, f2} specifies safety factors to be included in the comparison of estimates
k.
An order decrease is made when
k-1<f1k.
An order increase is made when
k+1<f2k.
There are some additional restrictions, such as when the maximal order increase per step is one (two for symmetric methods), and when an increase in order is prevented immediately after a rejected step.
For a nonstiff base method the default values are {4/5, 9/10} whereas for a stiff method they are {7/10, 9/10}.
Option summary
| | |
| ExtrapolationSequence | Automatic | specify the sequence to use in extrapolation |
| MaxDifferenceOrder | Automatic | specify the maximum order to use |
| Method | ExplicitModifiedMidpoint | specify the base integration method to use |
| MinDifferenceOrder | Automatic | specify the minimum order to use |
| OrderSafetyFactors | Automatic | specify the safety factors to use in the estimates for adaptive order selection |
| StartingDifferenceOrder | Automatic | specify the initial order to use |
| StepSizeRatioBounds | Automatic | specify the bounds on a relative change in the new step size hn+1 from the current step size hn as low ≤hn+1/hn≤ high |
| StepSizeSafetyFactors | Automatic | specify the safety factors to incorporate into the error estimate used for adaptive step sizes |
| StiffnessTest | Automatic | specify whether to use the stiffness detection capability |
Options of the method Extrapolation.
The default setting of
Automatic for the option
ExtrapolationSequence selects a sequence based on the stiffness and symmetry of the base method.
The default setting of
Automatic for the option
MaxDifferenceOrder bounds the maximum order by two times the decimal working precision.
The default setting of
Automatic for the option
MinDifferenceOrder selects the minimum number of two extrapolations starting from the order of the base method. This also depends on whether the base method is symmetric.
The default setting of
Automatic for the option
OrderSafetyFactors uses the values
{7/10, 9/10} for a stiff base method and
{4/5, 9/10} for a nonstiff base method.
The default setting of
Automatic for the option
StartingDifferenceOrder depends on the setting of
MinDifferenceOrder pmin. It is set to
pmin+1 or
pmin+2 depending on whether the base method is symmetric.
The default setting of
Automatic for the option
StepSizeRatioBounds uses the values
{1/10, 4} for a stiff base method and
{1/50, 4} for a nonstiff base method.
The default setting of
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.
The default setting of
Automatic for the option
StiffnessTest indicates that the stiffness test is activate if a nonstiff base method is used.
| | |
| StabilityCheck | True | specify whether to carry out a stability check on consecutive implicit solutions (see e.g. (44)) |
Options of the method LinearlyImplicitEuler, LinearlyImplicitMidpoint, and LinearlyImplicitModifiedMidpoint.
