1 General linear multistep methods

Euler’s method and the trapezium method are both special cases of a wider class of so-called linear multistep methods. The following Key Point gives the most general situation that we will look at.

The next Example shows one such choice.

Example 5

Write down the linear multistep scheme defined by the choices $k=1$ , ${\alpha }_0=-1$ , ${\alpha }_1=1$ , ${\beta }_0={\beta }_1=\frac{1}{2}$ .

Solution

Here $k=1$ so that

${\alpha }_1{y}_{n+1}+{\alpha }_0{y}_n=h\left({\beta }_1{f}_{n+1}+{\beta }_0{f}_n\right)$

and substituting the values in for the four coefficients gives

$y_{n+1}-y_n=h\left(\frac{1}{2}{f}_{n+1}+\frac{1}{2}{f}_n\right)$

which, as we know, is the trapezium method.

Task!

Write down the linear multistep scheme defined by the choices $k=1$ , ${\alpha }_0=-1$ , ${\alpha }_1=1$ , ${\beta }_0=1$ and ${\beta }_1=0$ .

Answer

Write down the linear multistep scheme defined by the choices $k=2$ , ${\alpha }_0=0$ , ${\alpha }_1=-1$ , ${\alpha }_2=1$ , ${\beta }_2=0$ , ${\beta }_1=\frac{3}{2}$ and ${\beta }_0=-\frac{1}{2}$ .

Answer

In the preceding Section we saw several examples implementing the Euler and trapezium methods.

The next Example deals with the explicit 2-step that was the subject of the Task above.

Example 6

A numerical scheme has been used to approximate the solution of

$\frac{dy}{dt}=t+y,y\left(0\right)=3$

and has produced the following estimates, to 6 decimal places,

$y\left(0.4\right)\approx 4.509822,y\left(0.45\right)\approx 4.755313$

Now use the 2-step, explicit linear multistep scheme

$y_{n+2}-y_{n+1}=h\left(1.5f_{n+1}-0.5f_n\right)$

to approximate $y\left(0.5\right)$ .

Solution

Evidently the value $h=0.05$ will serve our purposes and we seek $y_{10}\approx y\left(0.5\right)$ . The values we will need to use in our implementation of the 2-step scheme are $y_9=4.755313$ and

$f_9=f\left(0.45,y_9\right)=5.205313f_8=f\left(0.4,y_8\right)=4.909822$

to 6 decimal places since $f\left(t,y\right)=t+y$ . It follows that

And we conclude that $y\left(0.5\right)\approx 5.022966$ , where this approximation has been given to 6 decimal places.

Notice that in this implementation of a 2-step method we needed to use the values of the two $y$ values preceding the one currently being sought. Both $y_8$ and $y_9$ were used in finding $y_{10}$ .

Similarly, a $k$ -step method will use, in general, $k$ previous $y$ values at each time step.

This means that there is an issue to be resolved in implementing methods that are 2- or higher-step, because when we start we are only given one starting value $y_0$ . This issue will be dealt with towards the end of this Section. The following exercise involves a 2-step method, but (like the example above) it does not encounter the difficulty relating to starting values as it assumes that the numerical procedure is already underway.

Task!

A numerical scheme has been used to approximate the solution of

$\frac{dy}{dt}=t∕yy\left(0\right)=-2$

and has produced the following estimates, to 6 decimal places,

$y\left(0.24\right)\approx -2.013162,y\left(0.26\right)\approx -2.015546$

Now use the 2-step, explicit linear multistep scheme

$y_{n+2}-\frac{1}{2}{y}_{n+1}-\frac{1}{2}{y}_n=\frac{3}{2}h{f}_{n+1}$

to approximate $y\left(0.28\right)$ .

Answer

1.1 Zero stability

We now begin to classify linear multistep methods. Some choices of the coefficients give rise to schemes that work well, and some do not. One property that is required if we are to obtain reliable approximations is that the scheme be zero stable . A scheme that is zero stable will not produce approximations which grow unrealistically with $t$ .

We define the first characteristic polynomial

$\rho \left(z\right)={\alpha }_0+{\alpha }_{1}z+{\alpha }_2{z}^2+\dots {\alpha }_k{z}^k$

where the ${\alpha }_i$ are the coefficients of the linear multistep method as defined in Key Point 8 (page 21). This polynomial appears in the definition of zero stability given in the following Key Point.

The second characteristic polynomial is defined in terms of the coefficients on the right-hand side (the ${\beta }_j$ ), but its use is beyond the scope of this Workbook.

Example 7

Find the roots of the first characteristic polynomial for each of the examples below and determine whether or not the method is zero stable.

1. $y_{n+1}-y_n=h{f}_n$
2. $y_{n+1}-2y_n=h{f}_n$
3. $y_{n+2}+3y_{n+1}-4y_n=h\left(2f_{n+2}+f_{n+1}+2f_n\right)$
4. $y_{n+2}-y_{n+1}=\frac{3}{2}h{f}_{n+1}$
5. $y_{n+2}-2y_{n+1}+y_n=h\left(f_{n+2}-f_n\right)$
6. $y_{n+2}+2y_{n+1}+5y_n=h\left(f_{n+2}-f_{n+1}+2f_n\right)$

Solution

1. In this case $\rho \left(z\right)=z-1$ and the single zero of $\rho$ is $z=1$ . This is a simple (that is, not repeated) root with magnitude equal to 1, so the method is zero stable.
2. $\rho \left(z\right)=z-2$ which has one zero, $z=2$ . This has magnitude $2>1$ and therefore the method is not zero stable.
3. $\rho \left(z\right)=z^2+3z-4=\left(z-1\right)\left(z+4\right)$ . One root is $z=-4$ which has magnitude greater than 1 and the method is therefore not zero stable.
4. Here ${\alpha }_2=1$ , ${\alpha }_1=-1$ and ${\alpha }_0=0$ , therefore

   $\rho \left(z\right)=z^2-z=z\left(z-1\right)$

   which has two zeros, $z=0$ and $z=1$ . These both have magnitude less than or equal to 1 and there is no repeated zero with magnitude equal to 1, so the method is zero stable.

5. $\rho \left(z\right)=z^2-2z+1={\left(z-1\right)}^2$ . Here $z=1$ is not a simple root, it is repeated and, since it has magnitude equal to 1, the method is not zero stable.
6. $\rho \left(z\right)=z^2+2z+5$ and the roots of $\rho \left(z\right)=0$ can be found from the quadratic formula. In this case the roots are complex and are equal to $$ Zero-stability requires that the absolute values have magnitude less than or equal to 1. Consequently we conclude that the method is not zero stable.

Task!

Find the roots of the first characteristic polynomial for the linear multistep scheme

$y_{n+2}-2y_{n+1}+y_n=h\left(f_{n+2}+2f_{n+1}+f_n\right)$

and hence determine whether or not the scheme is zero stable.

Answer

At this stage, the notion of zero stability is rather abstract, so let us try using a zero unstable method and see what happens. We consider the simple test problem

$\frac{dy}{dt}=-y,y\left(0\right)=1$

which we know to have analytic solution $y\left(t\right)=e^{-t}$ , a quantity which decays with increasing $t$ . Implementing the zero unstable scheme

$y_{n+1}-2y_n=h{f}_n$

on a spreadsheet package with $h=0.05$ gives the following results

$\begin{array}{ccc}\hfill n& \hfill t=nh\hfill & \hfill y_n\approx y\left(nh\right)\\ ̲& ̲& ̲\\ \hfill 0& \hfill 0.00\hfill & \hfill 1.00000\\ \hfill 1& \hfill 0.05\hfill & \hfill 1.95000\\ \hfill 2& \hfill 0.10\hfill & \hfill 3.80250\\ \hfill 3& \hfill 0.15\hfill & \hfill 7.41488\\ \hfill 4& \hfill 0.20\hfill & \hfill 14.45901\\ \hfill 5& \hfill 0.25\hfill & \hfill 28.19506\\ \hfill 6& \hfill 0.30\hfill & \hfill 54.98037\\ \hfill 7& \hfill 0.35\hfill & \hfill 107.21172\\ \hfill 8& \hfill 0.40\hfill & \hfill 209.06286\\ \hfill 9& \hfill 0.45\hfill & \hfill 407.67258\\ \hfill 10& \hfill 0.50\hfill & \hfill 794.96153\\ \hfill 11& \hfill 0.55\hfill & \hfill 1550.17499\\ \hfill 12& \hfill 0.60\hfill & \hfill 3022.84122\\ \hfill 13& \hfill 0.65\hfill & \hfill 5894.54039\\ \hfill 14& \hfill 0.70\hfill & \hfill 11494.35376\\ \hfill 15& \hfill 0.75\hfill & \hfill 22413.98982\end{array}$

where 5 decimal places have been given for $y_n$ . The dramatic growth in the values of $y_n$ is due to the zero instability of the method. (There are in fact other things than zero instability wrong with the scheme $y_{n+1}-2y_n=h{f}_n$ , but it is the zero instability that is causing the large numbers.)

1.2 Consistency and order

A scheme that is zero stable will produce approximations that do not grow in size in a way that is not present in the exact, analytic solution. Zero stability is a required property, but it is not enough on its own. There remains the issue of whether the approximations are close to the exact values.

The truncation error of the general linear multistep method is a measure of how well the differential equation and the numerical method agree with each other. It is defined by

${\tau }_j=\frac{1}{\beta }\left(\frac{c_0}{h}y\left(jh\right)+c_1{y}^{\prime }\left(jh\right)+c_{2}h{y}^{″}\left(jh\right)+c_3{h}^2{y}^{‴}\left(jh\right)+\dots \right)=\frac{1}{\beta h}{\sum }_{p=0}^{\infty }{c}_p{h}^p{y}^{\left(p\right)}\left(jh\right)$

where $\beta =\sum {\beta }_j$ is a normalising factor.

It is the first few terms in this expression that will matter most in what follows, and it helps us that there are formulae for the coefficients which appear

$c_0=\sum {\alpha }_j,c_1=\sum \left(j{\alpha }_j-{\beta }_j\right),c_2=\sum \left(\frac{j^2}{2}{\alpha }_j-j{\beta }_j\right),c_3=\sum \left(\frac{j^3}{3!}{\alpha }_j-\frac{j^2}{2}{\beta }_j\right)$

and so on, the general formula for $p\ge 2$ is $c_p=\sum \left(\frac{j^p}{\left(p\right)!}{\alpha }_j-\frac{j^{p-1}}{\left(p-1\right)!}{\beta }_j\right)$ .

Recall that the truncation error is intended to be a measure of how well the differential equation and its approximation agree with each other. We say that the numerical method is consistent with the differential equation if ${\tau }_j$ tends to zero as $h\to 0$ . The following Key Point says this in other words.

Example 8

Show that Euler’s method ( $y_{n+1}=y_n+h{f}_n$ ) is consistent.

Solution

In this case ${\alpha }_1=1$ , ${\alpha }_0=-1$ , ${\beta }_1=0$ and ${\beta }_0=1$ . It follows that

$c_0=\sum {\alpha }_j=1-1=0\text{and}{c}_1=\sum j{\alpha }_j-{\beta }_j=1{\alpha }_1-\left({\beta }_0+{\beta }_1\right)=1-\left(1+0\right)=0$

and therefore Euler’s method is consistent.

Task!

Show that the trapezium method ( $y_{n+1}=y_n+\frac{h}{2}\left(f_{n+1}+f_n\right)$ ) is consistent.

Answer

Task!

Determine the consistency (or otherwise) of the following 2-step linear multistep schemes

1. $y_{n+2}-2y_{n+1}+y_n=h\left(f_{n+2}-f_n\right)$
2. $y_{n+2}-y_{n+1}=h\left(f_{n+1}-2f_n\right)$
3. $y_{n+2}-y_{n+1}=h\left(2f_{n+2}-f_{n+1}\right)$

Answer

(Notice also that the first characteristic polynomial $\rho \left(z\right)$ , defined on page 6 of this Section, evaluated at $z=1$ is equal to ${\alpha }_0+{\alpha }_1+\cdots +{\alpha }_k=c_0$ . It follows that a consistent scheme must always have $z=1$ as one of the roots of its $\rho \left(z\right)$ .)

Assuming that the method is consistent, the order of the scheme tells us how quickly the truncation error tends to zero as $h\to 0$ . For example, if $c_0=0$ , $c_1=0$ , $c_2=0$ and $c_3\ne 0$ then the first non-zero term in ${\tau }_j$ will be the one involving $h^2$ and the linear multistep method is called second-order . This means that if $h$ is small then ${\tau }_j$ is dominated by the $h^2$ term (because the $h^3$ and subsequent terms will be tiny in comparison) and halving $h$ will cause ${\tau }_j$ to decrease by a factor of approximately $\frac{1}{4}$ . The decrease is only approximately known because the $h^3$ and other terms will have a small effect.

We summarise the general situation in the following Key Point.

Combining the last two Key Points gives us another way of describing consistency: “A linear multistep method is consistent if it is at least first order".

Example 9

Find the order of

1. Euler’s method
2. The trapezium method.

Solution

1. We have already found that $c_0=c_1=0$ so the first quantity to calculate is

   $c_2=\sum \left(\frac{j^2}{2}{\alpha }_j-j{\beta }_j\right)=\frac{1}{2}{\alpha }_1-{\beta }_1=\frac{1}{2}$

   which is not zero and therefore Euler’s method is of order 1. (Or, in other words, Euler’s method is first order.)

2. We have already found that $c_0=c_1=0$ so the first quantity to calculate is

   $c_2=\sum \left(\frac{j^2}{2}{\alpha }_j-j{\beta }_j\right)=\frac{1}{2}{\alpha }_1-{\beta }_1=\frac{1}{2}-\frac{1}{2}=0$

   this is equal to zero, so we must calculate the next coefficient

   $c_3=\sum \left(\frac{j^3}{3!}{\alpha }_j-\frac{j^2}{2}{\beta }_j\right)=\frac{1}{6}{\alpha }_1-\frac{1}{2}{\beta }_1=\frac{1}{6}-\frac{1}{4}=-\frac{1}{12}$

   which is not zero. Hence the trapezium method is of order 2 (that is, it is second order).

This finally explains some of the results we saw in the first Section of this Workbook. We saw that the errors incurred by the Euler and trapezium methods, for a particular test problem, were roughly proportional to $h$ and $h^2$ respectively. This behaviour is dictated by the first non-zero term in the truncation error which is the one involving $c_{2}h$ for Euler and the one involving $c_3{h}^2$ for trapezium.

We now apply the method to another linear multistep scheme.

Example 10

Find the order of the 4-step, explicit linear multistep scheme

$y_{n+4}-y_{n+3}=\frac{h}{24}\left(55f_{n+3}-59f_{n+2}+37f_{n+1}-9f_n\right)$

Solution

In the established notation we have ${\alpha }_4=1$ , ${\alpha }_3=-1$ , ${\alpha }_2=0$ , ${\alpha }_1=0$ and ${\alpha }_0=0$ . The $\beta$ terms similarly come from the coefficients on the right hand side (remembering the denominator of 24).

Now

$c_0=\sum {\alpha }_j=0\text{and}{c}_1=\sum j{\alpha }_j-{\beta }_j=0$

from which we conclude that the method is consistent.

We also find that

(The exact value of $c_5$ is $\frac{251}{720}$ .)

Because $c_5$ is the first non-zero coefficient we conclude that the method is of order 4.

So the scheme in Example 10 has the property that the truncation error will tend to zero proportional to $h^4$ (approximately) as $h\to 0$ . This is a good thing, as it says that the error will decay to zero very quickly, when $h$ is decreased.

Task!

Find the order of the 2-step linear multistep scheme

$y_{n+2}-y_{n+1}=\frac{h}{12}\left(f_{n+2}+8f_{n+1}-f_n\right)$

Answer

1.3 Convergence

The key result concerning linear multistep methods is given in the following Key Point.

The proof of this result lies beyond the scope of this Workbook. It is worth pointing out that this is not the whole story. The convergence result is useful, but only deals with $h$ as it tends to zero. In practice we use a finite, non-zero value of $h$ and there are ways of determining how big an $h$ it is possible to “get away with" for a particular linear multistep scheme applied to a particular initial value problem.

If, when implementing the methods described above, it is found that the numerical approximations behave in an unexpected way (for example, if the numbers are very large when they should not be, or if decreasing $h$ does not seem to lead to results that converge) then one topic to look for in further reading is that of “absolute stability".