Extrapolation IntroductionExtrapolation methods are a class of arbitraryorder 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 higherorder 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 RungeKuttatype methods. For high precision, however, the arbitrary order means that they can be arbitrarily faster than fixedorder 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. In[3]:= 
ExtrapolationThe 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 higherorder 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
by
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 AitkenNeville 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 sequencesAny extrapolation sequence can be specified in the implementation. Some common choices are as follows. This is the Romberg sequence. In[5]:= 
Out[5]=

This is the Bulirsch sequence. In[6]:= 
Out[6]=

This is the Harmonic sequence. In[7]:= 
Out[7]=

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. In[8]:= 
Out[8]=

Here is an example of adding a function to define the Harmonic sequence is where the method order is an optional pattern. In[9]:= 
Rounding error accumulationFor highorder extrapolation an important consideration is the accumulation of rounding errors in the AitkenNeville 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 AitkenNeville algorithm (1) on the initial data using the Harmonic sequence and a symmetric ordertwo base integration method.
Hence, for an ordersixteen 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 reductionIt seems worthwhile to look for approaches that can reduce the effect of rounding errors in highorder 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 baseintegration method in order to reduce the magnitude of the rounding errors in floatingpoint operations. This approach, based on ideas that dated back to [G51] (see also [K65], [M65a], [M65b], [V79]), are explained next. Base methodsThe following methods are the most common choices for base integrators in extrapolation.
ExplicitEuler
ExplicitMidpoint
ExplicitModifiedMidpoint (Gragg smoothing step)
LinearlyImplicitEuler
LinearlyImplicitMidpoint (BaderDeuflhard formulation without smoothing step) LinearlyImplicitModifiedMidpoint (BaderDeuflhard 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 stepsGiven , and H, consider a succession of integration steps with step size carried out using Euler's method: where . Correspondence with explicit RungeKutta methodsIt is well known that, for certain base integration schemes, the entries in the extrapolation table produced from (1) correspond to explicit RungeKutta methods (see Exercise 1, Section II.9 in [HNW93]). For example, (1) is equivalent to an stage explicit RungeKutta method: where the coefficients are represented by the Butcher table: ReformulationLet = . Then the integration (1) can be rewritten to reflect the correspondence with an explicit RungeKutta method (1, 1) as: where terms in the righthand 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 finiteprecision floatingpoint arithmetic. Multiple linearly implicit Euler stepsNext consider a succession of integration steps carried out using the linearly implicit Euler method. Increments arise naturally in the description of many semiimplicit 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 righthand side. The desired result is obtained from (1) as:
ReformulationReformulation 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]). ReformulationReformulation of (1) can be accomplished in terms of increments as: Gragg's smoothing stepThe 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 finiteprecision 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 regionsThis illustrates the effect of the smoothing step on the linear stability domain (carried out using the package NumericalMath`OrderStar`).
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: 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). 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 RungeKutta method. Therefore, it appears that the correspondence has not yet been fully exploited and further refinement is possible. Since the AitkenNeville algorithm (1) involves linear differences, the entire extrapolation process can be carried out using increments. This leads to the following modification of the AitkenNeville 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. There are a number of important implementation issues that should be considered, some of which are mentioned here. 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. 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]. 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. 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 righthand 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. One of the simplest nonlinear equations describing a circuit is van der Pol's equation. In[11]:= 
This solves the equations using Extrapolation with the LinearlyImplicitEuler base method with the default subHarmonic sequence 2,3 4,.... In[14]:= 
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. In[15]:= 
This plots the step sizes taken in computing the solution. In[16]:= 
Select the Lorenz equations. In[17]:= 
This invokes a bigfloat, or software floatingpoint number, embedded explicit RungeKutta method of order 9(8) [V78]. In[18]:= 
Out[18]=

This invokes the Adams method using a bigfloat version of LSODA. The maximum order of these methods is twelve. In[19]:= 
Out[19]=

This invokes the Extrapolation method with ExplicitModifiedMidpoint as the base integration scheme. In[20]:= 
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. In[21]:= 
Out[23]=

Here are the solutions obtained by the various methods at the end of the integration. In[24]:= 
Out[24]//DisplayForm=

For comparison, here is a higherprecision reference solution obtained using Extrapolation. In[25]:= 
Out[26]=

Workerror comparisonFor comparing different extrapolation schemes consider an example from [HW96]. In[27]:= 
The exact solution is given by y(t) = 1/cos(t). Increment formulationThis example involves an eighthorder extrapolation of ExplicitEuler with the Harmonic sequence. Approximately two digits of accuracy are gained by using the incrementbased 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 incrementbased 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 doubleprecision floatingpoint 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 comparisonThis 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 orderdependent factor and an orderindependent factor. The option StepSizeRatioBounds>{
,
} specifies bounds on the next step size to take such that:
OrderSafetyFactorsAn 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 summaryOptions of the method Extrapolation.Options of the method LinearlyImplicitEuler, LinearlyImplicitMidpoint, and LinearlyImplicitModifiedMidpoint.
