1 Predictor-corrector methods

We have seen that when using an implicit linear multistep method there is an additional difficulty because we cannot, in general, solve simply for the newest approximate $y$ -value $y_{n+k}$ . A general $k$ -step implicit method involves, at the $k^{\text{th}}$ time step,

$\begin{array}{cccc}\hfill {\alpha }_k{y}_{n+k}+\hfill & \hfill \cdots +{\alpha }_1{y}_1+{\alpha }_0{y}_0=\hfill & \hfill h\left({\beta }_k{f}_{n+k}+\right.\hfill & \hfill \left.\cdots +{\beta }_1{f}_1+{\beta }_0{f}_0\right)\hfill \\ \hfill \uparrow \hfill & \hfill \hfill & \hfill \uparrow \hfill \\ \hfill \text{the}\hfill & \hfill \hfill & \hfill y_{n+k}\text{ occurs}\hfill \\ \hfill \text{unknown}\hfill & \hfill \hfill & \hfill \text{here too}\hfill \end{array}$

and if $f$ depends on $y$ in a complicated way then it is not obvious how to dig $y_{n+k}$ out of $f_{n+k}=f\left(\left(n+k\right)h,y_{n+k}\right)$ .

One solution to this problem would be to only ever use explicit methods in which ${\beta }_k=0$ . But this is not a good solution, for implicit methods generally have better properties than the explicit ones (for example, the implicit trapezium is second order while the explicit Euler is only first order).

Another solution involves a so-called predictor-corrector method. This involves

1. The predictor step. We use an explicit method to obtain an approximation $y_{n+k}^P$ to $y_{n+k}$ .
2. The corrector step. We use an implicit method, but with the predicted value $y_{n+k}^P$ on the right-hand side in the evaluation of $f_{n+k}$ . We use $f_{n+k}^P$ to denote this approximate (predicted) value of $f_{n+k}$ .
3. We can then go on to correct again and again. At each step we put the latest approximation to $y_{n+k}$ in the right-hand side of the scheme ( via $f$ ) to generate a new approximation from the left-hand side.

(This is not unlike an implementation of Newton-Raphson. In that method we require an initial guess (we “predict") and then the Newton-Raphson approach tells us how to iterate (or “correct") our latest approximation. The main difference here is that we have a systematic way of obtaining the initial prediction.)

It is sufficient for our purposes to illustrate the idea of a predictor-corrector method using the simplest possible pair of methods. We use Euler’s method to predict and the trapezium method to correct.

Example 12

Suppose that $y=y\left(t\right)$ is the solution to the initial value problem

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

Use Euler’s method and the trapezium method as a predictor-corrector pair (with one correction at each time step). Take the time step to be $h=0.05$ so as to obtain approximations to $y\left(0.05\right)$ and $y\left(0.1\right)$ .

Solution

Euler’s method, $y_{n+1}=y_n+h{f}_n$ , is the explicit method so we use that to predict . For the first time step we require $f_0=f\left(0,y_0\right)=f\left(0,3\right)=3$ and therefore

$y_1^P=y_0+h{f}_0=3+0.05\times 3=3.15$

We now use this predicted value of $y_1$ to obtain a “predicted" value for $f_1$ which we can use in the implicit trapezium method. We find $f_1^P=f\left(h,y_1^P\right)=f\left(0.05,3.15\right)=3.2$ . We now correct using the trapezium method in the form

$y_1=y_0+\frac{h}{2}\left(f_0+f_1^P\right)=3+\frac{1}{2}\left(0.05\right)\left(3+3.2\right)=3.155$

This completes prediction and one correction for the first time step.

For the second time step we require $f_1=f\left(h,y_1\right)=f\left(0.05,3.155\right)=3.205$ and therefore

$y_2^P=y_1+h{f}_1=3.155+0.05\times 3.205=3.31525$

which is the predicted value for $y_2$ . We now correct it with

$y_2=y_1+\frac{h}{2}\left(f_1+f_2^P\right)=3.155+\frac{1}{2}\left(0.05\right)\left(3.205+3.41525\right)=3.320506$

We conclude that

If correction is repeated until the corrected values settle down to a converged number then the approximation inherits all the (nice) properties of the implicit scheme. So, in the example above we would have second order accurate results obtained by a procedure which gets around the implicit nature of the trapezium method. Of course in the hand-calculations done above we only corrected once, rather than repeatedly to convergence.

The example above is such that the dependence of $f\left(t,y\right)$ on $y$ is very simple and we could use the approach seen in Section 32.1 to implement the trapezium method. It turns out that the true trapezium method approximations to $y\left(0.05\right)$ and $y\left(0.1\right)$ are $y_1=3.155128$ and $y_2=3.320776$ respectively, to 6 decimal places. The predictor-corrector method will produce these values if enough corrections are taken.

As noted in the last paragraph, the example above was one in which it is possible to get around the implicit nature of the trapezium method easily because of the simple way in which the right-hand side of the differential equation depends on $y$ . This is not true of the next example.

Example 13

Suppose that $y=y\left(t\right)$ is the solution to the initial value problem

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

Use Euler’s method and the trapezium method as a predictor-corrector pair (with one correction at each time step). Take the time step to be $h=0.2$ so as to obtain approximations to $y\left(0.2\right)$ and $y\left(0.4\right)$ .

Solution

Euler’s method, $y_{n+1}=y_n+h{f}_n$ , is the explicit method so we use that to predict. For the first time step we require $f_0=f\left(0,y_0\right)=f\left(0,1\right)=-1.55741$ and therefore

$y_1^P=y_0+h{f}_0=1+0.2\times -1.55741=0.688518$

We now use this predicted value to obtain a “predicted" value for $f_1$ which we can use in the implicit trapezium method. We find $f_1^P=f\left(h,y_1^P\right)=f\left(0.2,0.688518\right)=-0.82285$ . We now correct using the trapezium method in the form

$y_1=y_0+\frac{h}{2}\left(f_0+f_1^P\right)=1+\frac{1}{2}\left(0.2\right)\left(-1.55741-0.822848\right)=0.761974$

This completes prediction and one correction for the first time step.

For the second time step we require $f_1=f\left(h,y_1\right)=f\left(0.2,0.761974\right)=-0.95422$ and therefore

$y_2^P=y_1+h{f}_1=0.76194+0.2\times -0.95422=0.571131$

which is the predicted value for $y_2$ . We now correct it with

$y_2=y_1+\frac{h}{2}\left(f_1+f_2^P\right)=0.761974+\frac{1}{2}\left(0.2\right)\left(-0.95422--0.64257\right)=0.602296$

We conclude that

Task!

Suppose that $y=y\left(t\right)$ is the solution to the initial value problem

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

Use Euler’s method and the trapezium method as a predictor-corrector pair (with one correction at each time step). Take the time step to be $h=0.1$ so as to obtain approximations to $y\left(0.1\right)$ and $y\left(0.2\right)$ .

Answer

Exercise

Suppose that $y=y\left(t\right)$ is the solution to the initial value problem

$\frac{dy}{dt}=1∕\left(1+y^2\right)y\left(0\right)=1$

Use Euler’s method and the trapezium method as a predictor-corrector pair (with one correction at each time step). Take the time step to be $h=0.25$ so as to obtain approximations to $y\left(0.25\right)$ and $y\left(0.5\right)$ .

Answer

$y_1^P=y_0+h{f}_0=1+0.25\times 0.5=1.125$

We now use this predicted value to obtain a “predicted" value for $f_1$ which we can use in the implicit trapezium method. We find $f_1^P=f\left(h,y_1^P\right)=f\left(0.25,1.125\right)=0.441379$ . We now correct using the trapezium method in the form

$y_1=y_0+\frac{h}{2}\left(f_0+f_1^P\right)=1+\frac{1}{2}\left(0.25\right)\left(0.5+0.441379\right)=1.117672$

This completes prediction and one correction for the first time step.

For the second time step we require $f_1=f\left(h,y_1\right)=f\left(0.25,1.117672\right)=0.444604$ and therefore

$y_2^P=y_1+h{f}_1=1.125+0.25\times 0.444604=1.228823$

which is the predicted value for $y_2$ . We now correct it with

$y_2=y_1+\frac{h}{2}\left(f_1+f_2^P\right)=1.117672+\frac{1}{2}\left(0.25\right)\left(0.444604+0.398405\right)=1.223049$

We conclude that  $y\left(0.25\right)\approx 1.117672,y\left(0.5\right)\approx 1.223049$  to six decimal places.