This is documentation for Mathematica 5.2, which was
based on an earlier version of the Wolfram Language.
View current documentation (Version 11.2)



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.



The method DoubleStep performs a single application of Richardson's extrapolation for any one step integration method and is described within DoubleStep.

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  and compute the solution of the initial value problem by carrying out  steps with step size  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,  = 1,  = 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.



This is the Bulirsch sequence.



This is the Harmonic sequence.



A sequence that satisfies  has the effect of minimizing the roundoff errors for an order- base integration method.

For a base method of order two the first entries in the sequence are given by the following.



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  are disturbed with rounding errors  ,  ,  , ... and compute the propagation of these errors into the extrapolation table.

Due to the linearity of the extrapolation process (-1), suppose that the  are equal to zero and take  .

This shows the evolution of the Aitken-Neville algorithm (-1) 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 we will see later, it is more likely that rounding errors are made in  than in  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  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] (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  ,  and H, consider a succession of  integration steps with step size  carried out using Euler's method:

where  .

Correspondence with explicit Runge-Kutta methods

It is well known that, for certain base integration schemes, the entries  in the extrapolation table produced from (-1) correspond to explicit Runge-Kutta methods (see Exercise 1, Section II.9 in [HNW93]).

For example, (-1) is equivalent to an  -stage explicit Runge-Kutta method:

where the coefficients are represented by the Butcher table:


Let  =  . Then the integration (-1) can be rewritten to reflect the correspondence with an explicit Runge-Kutta method (-1, -1) as:

where terms in the right-hand side of (-1) are now considered as departures from the same value  .

The  in (-1) correspond to the  in (-1).

Let  , then the required result can be recovered as:

Mathematically the formulations (-1) and (-1, -1) are equivalent. For n > 1, however, the computations in (-1) 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 (-1) 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 (-1) as:


Reformulation in terms of increments as departures from  can be accomplished as follows:

The result for  using (-1) is obtained from (-1).

Notice that (-1) and (-1) 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 (-1) is computed with  midpoint steps, then the method has a symmetric error expansion ([G65], [S70]).


Reformulation of (-1) 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 (-1), the formulation (-1) is computed with  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 (-1) has an advantage over (-1) in finite-precision arithmetic because the values  , which typically have a larger magnitude than the increments  , 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  steps and ExplicitModifiedMidpoint uses  steps followed by the smoothing step (-1).

Stability regions

This illustrates the effect of the smoothing step on the linear stability domain (carried out using the package NumericalMath`OrderStar`).

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 (-1) is computed for  linearly implicit midpoint steps, then the method has a symmetric error expansion [BD83].


Reformulation of (-1) 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 (-1) rewritten in terms of increments becomes:

The required quantities are obtained when (-1) is run with  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  steps and LinearlyImplicitModifiedMidpoint uses  steps followed by the smoothing step (-1).

Polynomial extrapolation in terms of increments

We have seen how to modify  , the entries in the first column of the extrapolation table, in terms of increments.

However, for certain base integration methods, each of the  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 (-1) 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  in (-1) can be computed iteratively, starting from the initial quantities  that are obtained from the modified base integration schemes without adding the contribution from  .

The final desired value  can be recovered as  .

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  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  .

In order to interrupt computations in the computation of  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 (-1) differs from (-1) only in the second equation. It requires finding the solution for a different right-hand side but no extra function evaluation.

For the implicit midpoint method:  and  , which simply requires a few basic arithmetic operations.

Increments are a more accurate formulation for the implementation of both forms of stability check.


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,....



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.


This plots the step sizes taken in computing the solution.


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].



This invokes the Adams method using a bigfloat version of LSODA. The maximum order of these methods is twelve.



This invokes the Extrapolation method with ExplicitModifiedMidpoint as the base integration scheme.



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.



Here are the solutions obtained by the various methods at the end of the integration.



For comparison, here is a higher-precision reference solution obtained using Extrapolation.



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 (-1) are depicted in green.

The results for the increment formulation (-1) followed by standard extrapolation (-1) is depicted in blue.

The results for the increment formulation (-1) with extrapolation carried out on the increments using (-1) 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.

The increment formulation used throughout the extrapolation process produces rounding errors in  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.



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  constrains the choice of step size as follows. The step size chosen by the method for order  satisfies:

This includes both an order-dependent factor and an order-independent factor.


The option StepSizeRatioBounds->{


} specifies bounds on the next step size to take such that:


An important aspect in Extrapolation is the choice of order.

Each extrapolation of order  has an associated work estimate  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  and the expected new step size to take for a method of order  (computed from (-1)):  .

Comparing consecutive estimates,  allows a decision about when a different order method will be more efficient.

The option  specifies safety factors to be included in the comparison of estimates  .

An order decrease is made when   .

An order increase is made when  .

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

Options of the method Extrapolation.
Options of the method LinearlyImplicitEuler, LinearlyImplicitMidpoint, and LinearlyImplicitModifiedMidpoint.