Skip to article frontmatterSkip to article content

Implementation of multistep methods

We now consider some of the practical issues that arise when multistep formulas are used to solve IVPs. In this section we emphasize the vector IVP, u=f(t,u)\mathbf{u}'=\mathbf{f}(t,\mathbf{u}), and use boldface in the difference formula (6.6.2) as well.

6.7.1Explicit methods

As a concrete example, the AB4 method is defined by the formula

ui+1=ui+h(5524fi5924fi1+3724fi2924fi3),i=3,,n1.\mathbf{u}_{i+1} = \mathbf{u}_i + h\, ( \tfrac{55}{24}\mathbf{f}_i - \tfrac{59}{24} \mathbf{f}_{i-1} + \tfrac{37}{24}\mathbf{f}_{i-2} - \tfrac{9}{24}\mathbf{f}_{i-3}), \quad i=3,\ldots,n-1.

Function 6.7.1 shows a basic implementation of AB4.

Observe that Function 6.4.2 is used to find the starting values u1,u2,u3\mathbf{u}_1,\mathbf{u}_2,\mathbf{u}_3 that are needed before the iteration formula takes over. As far as RK4 is concerned, it needs to solve (the same step size as in the AB4 iteration). These results are then used to find f0,,f3\mathbf{f}_0,\ldots,\mathbf{f}_3 and get the main iteration started.

6.7.2Implicit methods

The implementation of an implicit multistep method is more difficult. Consider the second-order implicit formula AM2, also known as the trapezoid method. To advance from step ii to i+1i+1, we need to solve

z12hf(ti+1,z)=ui+12hf(ti,ui) \mathbf{z} - \tfrac{1}{2} h f(t_{i+1},\mathbf{z}) = \mathbf{u}_i + \tfrac{1}{2} h \mathbf{f}(t_i,\mathbf{u}_i)

for z\mathbf{z}. This equation can be written as g(z)=0\mathbf{g}(\mathbf{z})=\boldsymbol{0}, so the rootfinding methods of Chapter 4 can be used. The new value ui+1\mathbf{u}_{i+1} is equal to the root of this equation.

An implementation of AM2 using Function 4.6.3 from Quasi-Newton methods is shown in Function 6.7.2.

6.7.3Stiff problems

At each time step in Function 6.7.2, or any implicit IVP solver, a rootfinding iteration of uncertain expense is needed, requiring multiple calls to evaluate the function f\mathbf{f}. This fact makes the cost of an implicit method much greater on a per-step basis than for an explicit one. Given this drawback, you are justified to wonder whether implicit methods are ever competitive! The answer is emphatically yes, as Example 6.7.2 demonstrates.

Although the result of Example 6.7.2 may seem counter-intuitive, there is no contradiction. A fourth-order explicit formula is more accurate than a second-order implicit one, in the limit h0h\to 0. But there is another limit to consider, tt\to \infty with hh fixed, and in this one the implicit method wins.

Problems for which implicit methods are much more efficient than explicit counterparts are called stiff. A complete mathematical description will wait for Chapter 11, but a sure sign of stiffness is the presence of phenomena on widely different time scales. In Example 6.7.2, for instance, there are two slow periods during which the solution changes very little, interrupted by a very fast transition in the state. An explicit method “thinks” that the step size must always be dictated by the timescale of the fast transition, whereas an implicit method can take large steps during the slow periods.

6.7.4Adaptivity

As with RK methods, we can run two time stepping methods simultaneously in order to estimate the error and adjust the step size. For example, we could pair AB3 with AB4 as practically no cost because the methods differ only in how they include known information from the recent past. The more accurate AB4 value should allow an accurate estimate of the local error in the AB3 value, and so on.

Because multistep methods rely on the solution history, though, changing the step size is more algebraically complicated than for RK methods. If hh is changed, then the historical values ui1,ui2\mathbf{u}_{i-1},\mathbf{u}_{i-2}\ldots and fi1,fi2\mathbf{f}_{i-1},\mathbf{f}_{i-2}\ldots are no longer given at the right moments in time to apply the iteration formula. A typical remedy is to use interpolation to re-evaluate the historical values at the appropriate times. The details are important but not especially illuminating, and we do not give them here.

6.7.5Exercises