MATH 337, by T.

Lakoba, University of Vermont

21

Multistep, Predictor-Corrector, and Implicit methods

In this section, we will introduce methods that may be as accurate as high-order Runge-Kutta
methods but will require fewer function evaluations.
We will also introduce implicit methods, whose signicance will become clearer in a later
section.

3.1

Idea behind multistep methods

The gure on the right illustrates the (familiar) fact that if you know y (xi ), i.e. the slope
of y(x), then you can compute a rst-order
1st order
accurate approximation Yi+1
to the solution yi+1 .
Likewise, if you know the slope and the curvature of your solution at a given point, you
can compute a second-order accurate approx2nd order
imation, Yi+1
, to the solution at the next
step.

matches slope

matches slope
and curvature

Yi

ord
Y1st
i+1

ord
Y2nd
i+1

yi+1

Now, recall that curvature is proportional to y . This motivates the following.

Question: How can we nd approximation to yi using already computed values Yik ,
k = 0, 1, 2, . . . ?
Answer: Note that

y yi1
fi fi1
=
.
(3.1)
yi i
h
h
Here and below we will use the notation fi in two slightly dierent ways:
fi f (xi , yi )

or

fi f (xi , Yi )

(3.2)

whenever this does not cause any confusion.

Continuing with Eq. (3.1), we can state it more specically by writing
yi =

yi yi1
fi fi1
+ O(h) =
+ O(h) ,
h
h

(3.3)

where we will compute the O(h) term later. For now, we use (3.3) to approximate yi+1 as
follows:
h2
yi + hyi + yi + O(h3 )
2
)
(
h2 fi fi1
+ O(h) + O(h3 )
= yi + hfi +
2
h
(
)
3
1
=
yi + h
fi fi1 + O(h3 ) .
2
2

yi+1 = y(xi + h) =

Remark 1: To start the corresponding nite-dierence method, i.e.

)
(
1
3
fi fi1
Yi+1 = Yi + h
2
2

(3.4)

(3.5)

MATH 337, by T. Lakoba, University of Vermont

22

(now we use fi as f (xi , Yi )), one needs two initial points of the solution, Y0 and Y1 . These can
be computed, e.g., by the simple Euler method; this is discussed in more detail in Section 3.4.
Remark 2: Equation (3.4) becomes exact rather than approximate if y(x) = p2 (x) ax2 +bx+c
is a second-degree polynomial in x. Indeed, in such a case,
yi = 2axi + b,

and

yi = 2a =

yi yi1
;
h

(3.6)

(note the exact equality in the last formula). We will use this remark later on.
Method (3.5) is of the second order. If we want to obtain a third-order method along the
same lines, we need to use the third derivative of the solution:
yi =

yi 2yi1
+ yi2
+ O(h)
h2

(3.7)

(you will be asked to verify this equation in one of the homework problems). Then we proceed
as in Eq. (3.4), namely:
yi+1 = yi + hyi +

h2 h3
y + yi + O(h4 ) .
2 i
6

(3.8)

If you now try to substitute the expression on the r.h.s. of (3.3) for yi , you notice that you
actually need an expression for the O(h)-term there that would have accuracy of O(h2 ). Here
is the corresponding calculation:

yi yi1
y (xi ) y (xi1 )
=
h
h
[
]
2
yi yi hyi + h2 yi + O(h3 )
=
h
h
= yi yi + O(h2 ),
2

(3.9)

whence

yi yi1
h
+ yi + O(h2 ) .
(3.10)
h
2
To complete the derivation of the third-order nite-dierence method, we substitute Eqs.
(3.10), (3.7), and yi = fi etc. into Eq. (3.8). The result is:

yi =

Yi+1 = Yi +

h
[23fi 16fi1 + 5fi2 ] ;
12

(3.11)

the local truncation error of this method is O(h4 ). Method (3.11) is called the 3rd-order
AdamsBashforth method.
Similarly, one can derive higher-order AdamsBashforth methods. For example, the 4thorder AdamsBashforth method is
Yi+1 = Yi +

h
[55fi 59fi1 + 37fi2 9fi3 ] .
24

(3.12)

Methods like (3.5), (3.11), and (3.12) are called multistep methods. To start a multistep
method, one requires more than one initial point of the solution (in the examples considered
above, the number of required initial points equals the order of the method).

MATH 337, by T. Lakoba, University of Vermont

23

Comparison of multistep and RungeKutta methods

The advantage of multistep over single-step RK methods of the same accuracy is that the
multistep methods require only one function evaluation per step, while, e.g., the cRK method
requires 4, and the RKFehlberg method, 6, function evaluations.
The disadvantages of multistep methods are primarily related to the diculties that one
would face by making them use a varying step size. In contrast, for single-step RK methods,
this is a straightforward procedure, as we learned in Lecture 2. Below I list some of those
diculties for multistep methods.
1. Using a variable step size requires either an interpolation of the numerical solution from a
variable- onto a constant-step size grid, or using rather complicated coecients involving ratios
of step sizes, say hi /hi1 , instead of the constants like 3/2, 1/2 in formulae (3.5) etc.
2. If function f (x, y) changes abruptly at x = xchange , as in the problems in Homework # 2
involving an opening parachute, then a RK method would require only one step to adjust to the
new values of f beyond xchange . In contrast, a k-step method would take k steps to recover.
That is, over those k steps taken past the point xchange it would be using the information about
the pre-xchange values of f , which are irrelevant for what is happening beyond xchange .4
3. In Sec. 3.6 we will learn about error control in so-called predictorcorrector methods
(see Sec. 3.5), which is based on a key fact examplied by Eqs. (3.35), (3.36), and the pair
of equations below them. Namely, the local truncation errors of the predictor and corrector
equations must be proportional to each other. However, the derivation of this fact, part of
which is given in Sec. 3.8, hinges on the assumption that all steps have equal size.
Despite the aforementioned diculties, multistep methods with a variable step size have
been constructed and studied, although their theory goes far beyond the material covered in
this course.

3.2

An alternative way to derive formulae for multistep methods

Recall that the 2nd-order AdamsBashforth method (3.5) was exact on solutions y(x) that are
2nd-degree polynomials: y(x) = p2 (x) (see Remark 2 after Eq. (3.4)). Similarly, one expects
that the 3rd-order AdamsBashforth method should be exact for y(x) = p3 (x). We will now
use this observation to derive the formula for this method, Eq. (3.11), in a dierent manner
than in Sec. 3.1.
To begin, we take, according to the above observation, f (x, y) = y (x) = (p3 (x)) = p2 (x),
i.e. a 2nd-degree polynomial in x. We now integrate the dierential equation y = f (x, y) from
xi to xi+1 and obtain:

xi+1

yi+1 = yi +

f (x, y(x))dx .

(3.13)

xi

Let us approximate the integral by a quadrature formula, as follows:

xi+1
f (x, y(x))dx h(b0 fi + b1 fi1 + b2 fi2 )

(3.14)

xi

and require that the above equation hold exactly, rather than approximately, for any f (x, y(x)) =
p2 (x). This is equivalent to requiring that (3.14) hold exactly for f = 1, f = x, and f = x2 .
Without loss of generality5 , one can set xi = 0 and then rewrite Eq. (3.14) for the above three
4
5

This reason was suggested by Mr. Jacob WahlenStrothman, who took this course in 2012.
In a homework problem, you will be asked to show this.

MATH 337, by T. Lakoba, University of Vermont

24

forms of f :

for f = 1:

= h ( b0 1 + b1 1 + b2 1 )

1 dx

x dx

1 2
h
2

= h ( b0 0 + b1 (h) + b2 (2h) )

x2 dx =

1 3
h
3

= h ( b0 0 + b1 (h)2 + b2 (2h)2 ) .

for f = x:
0

for f = x :
0

(3.15)

Equations (3.15) constitute a linear system of 3 equations for 3 unknowns b0 , b1 , and b2 . Solving
it, we obtain
23
16
5
b0 = ,
b1 = ,
b2 = ,
12
12
12
which in combination with Eq. (3.14) yields the same method as (3.11). Methods of higher
order can be obtained similarly.

3.3

A more general form of multistep methods, with examples

The AdamsBashforth methods above have the following common form:

Yi+1 Yi = h

bk fik .

(3.16)

k=0

As has been shown in Sec. 3.2, the sum on the r.h.s. approximates
xi+1
f (x, y(x))dx .
xi

Let us now consider multistep methods of a more general form:

M

Yi+1

ak Yik = h

bk fik ,

(3.17a)

k=0

k=0

where

ak = 1.

(3.17b)

k=0

The last condition intuitively makes sense since the sum on the l.h.s. of (3.17a) replaces the
term 1 Yi in (3.16). Note that if we rewrite the l.h.s. of (3.17a) as
M

ak (Yi+1 Yik ),

k=0

where we have used (3.17b), then the r.h.s. of (3.17a) acquires an interpretation that is similar
to that of the r.h.s. of (3.16). Namely, it approximates
M

k=0

xi+1

f (x, y(x))dx .

ak
xik

MATH 337, by T. Lakoba, University of Vermont

25

In the next Lecture, we will discover that many methods of the form (3.17) have a serious
aw in them, but for now let us consider two particular examples, focusing only on the accuracy
of the following methods.
Simple center-dierence (Leap-frog) method
Recall that

yi yi1
= yi + O(h) .
h

(3.18)

However6 ,

yi+1 yi1
= yi + O(h2 ) .
(3.19)
2h
Thus, the l.h.s. of (3.19) provides a more accurate approximation to yi than does the l.h.s. of
(3.18). So we use Eq. (3.19) to produce a 2nd-order method:
Yi+1 = Yi1 + 2hfi ,

(3.20)

which is of the form (3.17). We need both Y0 and Y1 to start this method.
A divergent third-order method
(The term divergent will be explained in the next Lecture.)
Let us try to increase the order of method (3.20) from 2nd to 3rd by including extra terms
into the scheme:
Yi+1 (a0 Yi + a1 Yi1 + a2 Yi2 ) = b0 h fi ,
(3.21)
where we now require that the local truncation error of (3.21) be O(h4 ). We can follow the
derivation found either in Sec. 3.1 (Taylor-series expansion) or Sec. 3.2 (requiring that (3.21)
hold true for y = p3 (x)) to obtain the values of the coecients a0 through a2 , and b0 . The
result is:
3
1
Yi+1 + Yi 3Yi1 + Yi2 = 3h fi .
(3.22)
2
2
Supposedly, method (3.22) is more accurate than the Leap-frog method (3.20). However, we will
show in the next Lecture that method (3.22) is completely useless for numerical computations.

3.4

Starting a multistep method

To start any of the single-step methods, considered in Lectures 1 and 2, one only needs to know
the initial condition, Y0 = y0 , at x = x0 . To start any multistep method, one needs to know
the numerical solution at several points. For example, to start an AdamsBashforth method of
order m, one would need the values Y0 , . . . , Ym1 (see Eqs. (3.5), (3.11), and (3.12)). That is,
to start an mth-order method, one needs to know the solution at the rst m points. We will
now address the following question:
Suppose that we want to start a multistep method of order m using the values Y1 , . . . , Ym1
that have been computed by a starting (single-step) method of order n. What should the
order n of the starting method be so as not to compromise the order m of the multistep
method?
First, it is clear that if n m, then the local truncation error (LTE) made in the
computation of Y1 , . . . , Ym1 and of the terms on the r.h.s. of (3.16) and (3.17) will be at least
6

You will be asked to verify this.

MATH 337, by T. Lakoba, University of Vermont

26

as small (in the order of magnitude sense) as the LTE of the multistep method. So, using a
starting method whose order is no less than the order of the multistep method will not degrade
the accuracy of the latter method. But is it possible to use a starting method with n < m for
the same end result?
We will now show, using method (3.16) as an example, that it is possible to take n = m 1
(i.e., the starting methods order may be one less than the multistep methods order).7
The LTEs of Y1 through Ym1 are O(hn+1 ). Then the error contributed to Ym from the
second term (i.e., from Yi with i = m 1) on the l.h.s. of (3.16), is O(hn+1 ):
error of l.h.s. of (3.16) = O(hn+1 ).

(3.23)

Next, if fi through fiN on the r.h.s. were calculated using the exact solution y(x), then the
error of the r.h.s. would have been O(hm+1 ). Indeed, this error is just the LTE of method
xi+1
N

(3.16) that arises due to the approximation of

f (x, y(x))dx by h
bk fik . However,
xi

k=0

the fik s are calculated using values Y1 through Ym1 which themselves have been obtained
with the error O(hn+1 ) of the starting method. Then the error of each fik is also O(hn+1 ). 8
Therefore, combining the two underlined errors in the text above, one has:
error of r.h.s. of (3.16) = O(hm+1 ) + h O(hn+1 ) = max{O(hn+2 ), O(hm+1 )} .

(3.24)

Let us summarize. The error on the l.h.s. of (3.16), which occurs due to Yi (with i = m 1)
being computed by an nth-order starting method, is O(hn+1 ). The error on the r.h.s of (3.16)
is given by (3.24). Therefore:
combined error in (3.16) = O(hn+1 ) from l.h.s. + max{O(hn+2 ), O(hm+1 )} from r.h.s. = O(hn+1 ) .
(3.25)
(recall that we are only interested in the situation where n < m).
Now, in order not to decrease the accuracy of the multistep method, this error must satisfy
two criteria:
(i) It must have the same order of magnitude as the global error at the end of the computation, i.e., O(hm ); and in addition,
(ii) It may propagate to the next computed solution, i.e., to Yi+2 , but it must not accumulate
at each step with other errors of the same magnitude. (See an explanation below.)
One can easily see that criterion (i) is indeed satised for n+1 = m, i.e., when n = m1. As for
criterion (ii), it is also satised, although at rst sight that might seem strange. Indeed, we have
seen in Lecture 1 that if a LTE is O(hn+1 ), then the global error should be O(1/h) O(hn+1 ) =
O(hn ). However, this is not so for the situation considered here. While it is possible to show
this by a painstaking analysis of the orders of terms in Eqs. (3.16) and (3.23)(3.25), it is much
easier (although perhaps not as rigorous) to illustrate the situation with a gure, as shown
below.
7

Unfortunately, I was unable to find any detailed published proof of this result, and so the derivation found
below is my own. As such, it is subject to mistakes . However, a set of Matlab codes accompanying this
Lecture where the 3rd-order AdamsBashforth method (3.11) can be started using the modified Euler, Midpoint,
or simple Euler method, shows that if not this derivation itself, than at least its result is probably correct.
8
This can be seen by the following calculation. Let j i k and Yj y(xj ) = O(hn+1 ). Then |f (xj , Yj )
f (xj , y(xj ))| L|Yj y(xj )| = O(hn+1 ), where the inequality holds in view of the Lipshitz condition.

MATH 337, by T. Lakoba, University of Vermont

In this gure, y(x) is the exact solution, Y is the numerical solution obtained with a multistep method
of order m (where in the gure, m = 2 + 1 = 3), and
Ystart is the solution obtained by the starting method
for the rst m steps (recall that Y0 is the initial condition and thus does not need to be computed by the
starting method). The dashed line shows the ctitious solution, denoted as Y , which starts at (Ystart )m
and
is computed under the hypothesis that the term

27

y(x)

O(h m)
O(h n+1)

xi+1

Ystart

f (x, y(x))dx is computed exactly rather than

xi

being approximated by h

bk fik . The magnitude

O(h n+1)

i= 0 1 2

k=0

order of the error generated by each of the methods

is shown next to the double arrows.
One can see that Y is just the exact solution of the ODE except that its initial condition
is shifted by Ystart (xm ) y(xm ) = O(hn+1 ). Consequently, it is intuitive to expect that the
corresponding nal solutions will dier by a similar amount. Thus, y(xfinal ) Y (xfinal ) =
O(hn+1 ), and hence at x = xfinal ,
y(x) Y = (y(x) Y ) + (Y Y ) = O(hn+1 ) + O(hm ).
Thus, for n = m 1, the global error will be O(hm ), as required.
Let us summarize. We have shown that in order to have a multistep method (3.16) of order
m, the order n of the starting single-step method should be no lower than (m 1). However,
in Sec. 3.6 we will see that there is a reason why one may want to use a single-step method of
order m (as opposed to (m 1)) to start a multistep method of order m.

3.5

Predictorcorrector methods: General form

Let us recall the Modied Euler method introduced in Lecture 1 and write it here using slightly
dierent notations:
p
Yi+1
= Yi + hfi
(
)
p
c
)
Yi+1
= Yi + 12 h fi + f (xi+1 , Yi+1
(3.26)
c
.
Yi+1 = Yi+1
We can interpret the above as follows: We rst predict the new value of the solution Yi+1 by
the rst equation, and then correct it by the second equation. Methods of this kind are called
predictorcorrector (PC) methods.
Question: What is the optimal relation between the orders of the predictor and corrector
equations?
Answer: The example of the Modied Euler method suggests that the order of the corrector should be one higher than that of the predictor. More precisely, the following theorem
holds:
Theorem If the order of the corrector equation is n, then the order of the corresponding
PC method is also n, provided that the order of the predictor equation is no less than n 1.
Proof We will assume that the global error of the corrector equation by itself is O(hn )
and the global error of the predictor equation by itself is O(hn1 ). Then we will prove that the
global error of the combined PC method is O(hn ).

MATH 337, by T. Lakoba, University of Vermont

28

The general forms of the predictor and corector equations are, respectively:
p
Yi+1

Predictor:

= YiQ + h

pk fik ,

(3.27)

p
ck fik + hc1 f (xi+1 , Yi+1
).

(3.28)

k=0

Corrector:

c
Yi+1

= YiD + h

k=0

In the above two equations, Q, D, N, M are some integer nonnegative numbers. (One of the
questions at the end of this Lecture asks you to represent Eq. (3.26) in the form (3.27), (3.28),
i.e. to give values for Q, D, N, M and the coecients pk s and ck s.)
As we have done in previous derivations, let us assume that all computed values Yik ,
k = 0, 1, 2, . . . coincide with the exact solution at the coresponding points: Yik = yik . Then
we can use the identity
xi+1
xi+1
see (3.13)

yi+1 = yiQ + (yi+1 yiQ ) = yiQ +

y (x)dx = YiQ +
f (x, y(x))dx
xiQ

xiQ

and rewrite Eq. (3.27) as:

p
Yi+1
= YiQ +

xi+1

f (x, y(x))dx+ h
xiQ

pk fik

xi+1

f (x, y(x))dx

p
Yi+1
= yi+1 +EP .

xiQ

k=0

(3.29)
Here EP is the error made by replacing the exact integral
xi+1
f (x, y(x))dx
xiQ

by the linear combination of fik s, found on the r.h.s. of (3.27). Since, by the condition of the
Theorem, the global error of the predictor equation is O(hn1 ), then the LTE EP has the order
of O(h(n1)+1 ) = O(hn ).
Similarly, Eq. (3.28) is rewritten as
(
)
p
c
Yi+1
= yi+1 + EC + hc1 f (xi+1 , Yi+1
) f (xi+1 , yi+1 ) .
(3.30)
Here EC is the error obtained by replacing the exact integral
xi+1
f (x, y(x))dx
xiD

by the quadrature formula

h

ck fik

k=1

(note that the lower limit of the summation is dierent from that in (3.28)!). The last term on
p
the r.h.s. of (3.30) occurs because, unlike all previously computed Yik s, the Yi+1
= yi+1 .
c
yi+1 = O(hn+1 ) in (3.30). By the
To complete the proof,9 we need to show that Yi+1
condition of the Theorem, the corrector equation has order n, and hence the LTE EC = O(hn+1 ).
9

At this point, you have probably forgotten what we are proving. Pause, re-read the Theorems statement,
and then come back to finish the reading.

MATH 337, by T. Lakoba, University of Vermont

29

Then all that remains to be estimated is the last term on the r.h.s. of (3.30). To that end, we
recall that f satises the Lipschitz condition with respect to y, whence
p
p
|f (xi+1 , Yi+1
) f (xi+1 , yi+1 )| L|Yi+1
yi+1 | = L |EP |,

(3.31)

where L is the Lipschitz constant. Combining Eqs. (3.30) and (3.31) and using the triangle
inequality (1.5), we nally obtain
c
|Yi+1
yi+1 | |EC | + hL|EP | = O(hn+1 ) + h O(hn ) = O(hn+1 ) ,

(3.32)

which proves that the PC method has the LTE of order n + 1, and hence is the nth-order
method.
q.e.d.
We now present two PC pairs that in applications are sometimes preferred10 over the
Modied Euler method. The rst pair is:
p
Predictor : Yi+1
= Yi + 12 h(3fi fi1 )
(
)
p
c
Corrector : Yi+1
,
= Yi + 12 h fi + fi+1

(3.33)

p
p
where fi+1
= f (xi+1 , Yi+1
). The order of the PC method (3.33) is two.
The other pair is:

Predictor: 4th-order AdamsBashforth

p
Yi+1
= Yi +

1
h(55fi
24

59fi1 + 37fi2 9fi3 ) .

Corrector: 4th-order AdamsMoulton

c
Yi+1
= Yi +

p
1
h(9fi+1
24

(3.34)

+ 19fi 5fi1 + fi2 ) ,

This PC method as a whole has the same name as its corrector equation: the 4th-order
AdamsMoulton.

3.6

Predictorcorrector methods: Error monitoring

An observation one can make from Eqs. (3.33) is that both the predictor and corrector equations have the order two (i.e. the LTEs of O(h3 )). In view of the Theorem of the previous
subsection, this may seem to be unnecessary. Indeed, the contribution of the predictors error
to the LTE is h O(h3 ) = O(h4 ) (see Eq. (3.32)), while the LTE of the corrector equation
itself (which determines that of the entire PC method) is only O(h3 ). There is, however, an
important consideration because of which method (3.33) may be preferred over the Modied
Euler. Namely, one can monitor the error size in (3.33), whereas the Modied Euler does not
give its user such a capability. Below we explain this statement in detail. A similar treatment
can be applied to the AdamsMoulton method (3.34).
The key fact is that the LTEs of the predictor and correction equations (3.33) are proportional to each other in the leading order:
p
yi+1 Yi+1
=

yi+1
10

c
Yi+1

5 3
h yi + O(h4 ),
12
1 3
12
h yi + O(h4 ) .

In the next Section we will explain why this is so.

(3.35)
(3.36)

MATH 337, by T. Lakoba, University of Vermont

30

For the readers information, the analogues of the above estimates for the AdamsMoulton
method (3.34) are:
4th-order AdamsMoulton method:
251 5 (5)
h yi + O(h6 ),
720
19 5 (5)
h yi + O(h6 ) .

720

p
yi+1 Yi+1
=
c
yi+1 Yi+1
=

We derive (3.36) in Appendix 1 to this Lecture, while the derivation of (3.35) is left as an
exercise. Here we only note that the derivation of (3.36) hinges upon the fact that
yi Yip = O(h3 ), which is guaranteed by (3.35). Otherwise, i.e. if yi Yip = O(h2 ), as
in the predictor for the Modified Euler method, the term on the r.h.s. of (3.36)
would not have had such a simple form. See Remarks 2 and 3 after Eq. (3.40) for details.
We will now explain how (3.35) and (3.36) can be used together to monitor the error of the
PC method (3.33). From (3.36) we obtain the error of the corrector equation:
|ci+1 |

1 3
h |yi | .
12

On the other hand, from Eqs. (3.35) and (3.36) together, we have
(
)
5
1
p
c
|Yi+1 Yi+1 |
+
h3 |yi |.
12 12

(3.37)

(3.38)

Thus, from (3.37) and (3.38) one can estimate the error via the dierence of the predicted and
corrected values of the solution:
1 p
c
|ci+1 | |Yi+1
Yi+1
|.
6

(3.39)

Moreover, Eqs. (3.35) and (3.36) can also be used to obtain a higher-order method than
(3.33), because they imply that
yi+1 =

)
1( p
c
Yi+1 + 5Yi+1
+ O(h4 ) .
6

Hence

)
1( p
c
Yi+1 + 5Yi+1
(3.40)
6
p
c
alone. (Note
produces a more accurate approximation to the solution than either Yi+1
or Yi+1
a similarity with the Romberg extrapolation described in Lecture 1.)
Thus, Eqs. (3.33), (3.39), and (3.40) can be used to program a PC method that allows the
user to monitor the error. Namely:
1. Compute Yip and Yic from (3.33). Compute the improved solution from (3.40).
c
using (3.39).
2. Estimate the error of Yi+1
The above procedure produces a 3rd-order-accurate solution (3.40) while monitoring the
error size of the associated 2nd-order method (3.33). This is analogous to the situation for the
adaptive RK methods described in Lecture 2, which computed an nth-order accurate solution
while controlling the error of a related method of the lower order, (n 1).
The following Remarks address various issues that arise when at least one of the methods
in the PC pair is multistep.
Yi+1 =

MATH 337, by T. Lakoba, University of Vermont

31

Remark 1: As we have just discussed, using the predictor and corrector equations of the
same order has the advantage of allowing one to monitor the error. However, it may have a
disadvantage of making such schemes less stable compared to schemes where the predictors
order is one less than that of the corrector. We will study the concept of stability in the next
Lecture.
Remark 2: Suppose we plan to use a PC method with the predictor and corrector equations having the same order, say m, so as to monitor the error, as described above. Furthermore,
suppose that the predictor and/or corrector equations are based on an AdamsBashforth multistep method, i.e. that of the form (3.16). (For example, the predictor equations in (3.33) and
(3.34) are of this form.) Let us now re-examine the question addressed in Sec. 3.4. Namely,
what order starting method should we use in the predictor equation so that the PC method
would be able to monitor the error?
In Sec. 3.4 we showed that a starting method of order (m 1) would suce to make the
predictor method have order m. More precisely, the global error of the predictor method is
O(hm ) in this case. However, the local truncation error of such a predictor method is also
O(hm ) and not O(hm+1 )! 11 For example, if we start the 2nd-order predictor method in (3.33)
by the simple Euler method, then the LTE of O(h2 ) made by the simple Euler in computing
Y1p , will propagate up to Yi for any i. Such an error will invalidate both Eqs. (3.35) and
(3.36); this can be seen from the derivation of (3.36) found in Appendix 1. This, in turn, will
make impossible the error monitoring similar to (3.39), as well as the implementation of the
(m + 1)th-order method (see (3.40)) based on the mth-order PC pair. On the other hand,
if the starting single-step methods LTE is O(hm+1 ), then its contribution to the LTE of the
multistep AdamsBashforth-type method is O(h hm+1 ) = O(hm+2 ) (see (3.32)). This would
not aect the (m + 1)th order of the multistep methods LTE.
Thus, we conclude: If you want to be able to monitor the error in a PC method
where the predictor and corrector equations are AdamsBashforth-type multistep
methods of order m, you need to start the predictor equation by an mth-order
starting method.
Remark 3:12 Let us now suppose that at least one of the multistep methods in the PC
pair is not of the AdamsBashforth type but instead has the more general form (3.17). We
will now argue that such a method, if it is of order m, is to be started with a single-step
method of order (m + 1) in order for the PC pair to be able to monitor the error (as well
as to be able to produce a solution of order (m + 1), as in (3.40)). The key is still to have
the LTEs of the predictor and corrector equations being of order O(hm+1 ) and also being
proportional to each other, as in (3.35) and (3.36). The second of these conditions implies that
the contribution of the LTE of the starting method must be of order O(hm+2 )! Indeed, if it
were of order O(hm+1 ), it would not translate into an O(h hm+1 ) = O(hm+2 ) quantity, as it
did for AdamsBashforth-type methods, but would stay as O(hm+1 ) due to the terms ak Yik
in (3.17a).
Thus, we conclude: If you want to be able to monitor the error in a PC method
where the predictor and corrector equations are multistep methods of the form
(3.17) and of order m, you need to start the predictor equation by an (m + 1)th-order
starting method.
Remark 4: At the end of Sec. 3.1 we listed disadvantages of using multistep methods with
11

This error made by the starting method would propagate, but not accumulate, to the last computed value
of the numerical solution, as we showed in Sec. 3.4.
12
This fact was pointed out by Mr. Tyler Gray, who took this course in Spring 2016.

MATH 337, by T. Lakoba, University of Vermont

32

a variable step size. Let us reiterate here one of them. Namely, the r.h.s.s of (3.35) and (3.36)
are not simply proportional to y but contain other terms if the step size is varied. Therefore,
the error monitoring for PC methods with a variable step size is a complicated (if feasible at
all) task.
Appendix 2 (Sec. 3.9) discusses a repeated use of the corrector equation within a given PC
method for the goal of increasing the overall accuracy of the PC method. This Appendix can
be skipped on the rst reading.
To summarize on the PC methods:
1. The PC methods may provide both high accuracy and the capability of error monitoring,
all at a potentially lower computational cost than RK-Fehlberg or RK-Merson methods.
For example, the AdamsMoulton method (3.34) has the error of the same (fourth) order
as the aforementioned RK methods, while requiring k + 1 function evaluations, where k
is the number of times one has to iterate the corrector equation. If k < 4, then AdamMoulton requires fewer function evaluations than either RK-Merson or RK-Fehlberg.
2. The adjustment of the step size in PC methods is awkward (as it is in all multistep
methods); see the end of Sec. 3.1 and Remark 3 above.
3. One may ask, why not just halve the step size of the AdamsBashforth method (which
would reduce the global error by a factor of 24 = 16, i.e. a lot) and then use it alone
without the AdamsMoulton corrector formula? The answer is this. First, one will
then lose the ability to monitor the error. Second, the AdamsBashforth may sometimes
produce a numerical solution which has nothing to do with the exact solution, while
the PC AdamsMoultons solution will stay close to the exact one. This issue will be
discussed in detail in the next Lecture.

3.7

Implicit methods

We noted in Lecture 1 that the simple Euler method is analogous to the left Riemann sums
when integrating the dierential equation y = f (x).
right Riemann sums

The method analogous to the right Riemann

sums is:
Yi+1 = Yi + hf (xi+1 , Yi+1 ) .

(3.41)

It is called the implicit Euler, or backward

Euler, method. This is a rst-order method:
Its global error is O(h) and the LTE is O(h2 ).
x

x1

x2

x3

We note that if f (x, y) = a(x)y +b(x), then the implicit equation (3.41) can be easily solved:
Yi+1 =

Yi + hbi+1
.
1 hai+1

(3.42)

MATH 337, by T. Lakoba, University of Vermont

33

However, for a general nonlinear f (x, y), equation (3.41) cannot be solved exactly, and its
solution then has to be found numerically, say, by the Newton-Raphson method.
Question: Why does one want to use the implicit Euler, which is so much harder to solve
than the simple Euler method?
Answer: Implicit methods have stability properties that are much better than those of
explicit methods (like the simple Euler). We will discuss this in the next Lecture.
Note that the last remark about AdamsBashforth vs. AdamsMoulton, found at the end
of the previous Section, is also related to the stability issue. Indeed, the AdamsBashforth
method (the rst equation in (3.34)) is explicit, and thus according to the above, it should be
not as stable as the AdamsMoulton method (the second equation in (3.34)), which is implicit
p
c
if one treats Yi+1
in it as being approximately equal to Yi+1
.
Finally, we present equations for the Modied implicit Euler method:
h
Yi+1 = Yi + (f (xi , Yi ) + f (xi+1 , Yi+1 )) .
(3.43)
2
This is a second-order method.

3.8

Appendix 1: Derivation of (3.36)

Here we derive the LTE of the corrector equation in the method (3.33). Assuming, as usual,
p
that Yi = yi , and using Yi+1
= yi+1 + O(h3 ) (since the order of the predictor equation is two),
one obtains from the corrector equation of (3.33):
yi + 12 h (yi + f (xi+1 , yi+1 + O(h3 )) )
yi + 12 h (yi + f (xi+1 , yi+1 ) + O(h3 ) )
(
)

=
yi + 12 h yi + yi+1
+ O(h3 )
(
[
]
)
= yi + 12 h yi + yi + hyi + 12 h2 yi + O(h3 ) + O(h3 )
=
yi + hyi + 21 h2 yi + 14 h3 yi + O(h4 ) .

c
Yi+1
=
=

On the other hand, for the exact solution we have the usual Taylor series expansion:
1
1
yi+1 = yi + hyi + h2 yi + h3 yi + O(h4 ) .
2
6
Subtracting (3.44) from (3.45), we obtain
1
c
yi+1 Yi+1
= h3 yi + O(h4 ) ,
12
which is (3.36).

3.9

(3.44)

(3.45)

Appendix 2: Repeated application of the corrector equation in

a PC method

Note that in a PC method, one can apply the corrector formula more than once. For example,
for the method (3.33), we will then have:
p
Yi+1

c,1
Yi+1

c,2
Yi+1

Yi+1

Yi + 21 h(3fi fi1 )
(
)
p
Yi + 12 h fi + f (xi+1 , Yi+1
) ,
(
)
c,1
Yi + 12 h fi + f (xi+1 , Yi+1
) ,

etc.
c,k
= Yi+1

(3.46)

MATH 337, by T. Lakoba, University of Vermont

34

Question: How many times should we apply the corrector equation?

We need to strike a compromise here. If we apply the corrector too many times, then we
will waste computer time if each iteration of the corrector changes the solution by less than the
LTE of the method. On the other hand, we may have to apply it more than once in order to
c,k
c,k1
make the dierence |Yi+1
Yi+1
| between the last two iterations much smaller than the LTE
of the corrector equation (since the latter error is basically the error of the method; see Eq.
(3.36)).
Ideally, one would like to know the conditions under which it is sucient to apply the
corrector equation only once, so that no benets would be gained by its successive applications.
Below we derive such a sucient condition for the method (3.33)/(3.46). For another PC
method, e.g., AdamsMoulton, an analogous condition can be derived along the same lines.
Suppose the maximum allowed global error of the solution is glob . The allowed LTE is then
about loc = hglob (see Sec. 2.2). We impose two requirements:
(i) The LTE of our solution should not exceed loc ;
c,2
c,1
(ii) The dierence |Yi+1
Yi+1
| should be much smaller than loc .
Requirement (i) is necessary to satisfy in order to obtain the required accuracy of the numerical
solution. Requirement (ii) is necessary to satisfy in order to use the corrector equation only
once.
Requirement (i) along with Eq. (3.39) yields
1 p
c,1
|Y Yi+1
| < loc ,
6 i+1

p
c,1
|Yi+1
Yi+1
| < 6loc .

(3.47)

If at some xi condition (3.47) does not hold, the step size needs to be reduced in accordance
p
c,1
with the errors order |Yi+1
Yi+1
| = O(h3 ).
Requirement (ii) implies:
[
(
)]
c,1
c,2
c,1
)
|Yi+1
Yi+1
| =  Yi + 12 h fi + f (xi+1 , Yi+1
[
(
)]
p
Yi + 12 h fi + f (xi+1 , Yi+1
) 


= 1 h f (xi+1 , Y c,1 ) f (xi+1 , Y p )
2

Lipschitz

i+1

 c,1

p 
1
hL Yi+1
Yi+1
2

(3.47)

1
hL
2

i+1

6loc .

(3.48)

c,2
c,1
Thus, a sucient condition for |Yi+1
Yi+1
| loc to hold is

1
hL 6loc loc ,
2

hL

or

1
.
3

(3.49)

If condition (3.49) is satised, then a single application of the corrector equation is adequate.
If, however, the step size is not small enough, we may require two iterations of the corrector.
Then a second application of (3.48) would produce the condition:
c,3
c,2
|Yi+1
Yi+1
| loc

)2
1
hL 6loc loc ,
2

which is less restrictive than (3.49).

or

hL

2
0.82 ,
3

(3.50)

MATH 337, by T. Lakoba, University of Vermont

3.10

35

Questions for self-assessment

1. Make sure you can reproduce the derivation of Eq. (3.4).

2. What is the idea behind the derivation of Eq. (3.5)?
3. Derive Eqs. (3.9) and (3.10).
4. Derive Eq. (3.11) as indicated in the text.
5. Describe two alternative ways to derive formulae for multistep methods.
6. Verify Eq. (3.19).
7. For a multistep method of order m, what should the order of the starting method be?
8. Convince yourself that method (3.26) is of the form (3.27) and (3.28).
9. What is the origin of the error EP in Eq. (3.29)?
10. What is the origin of the error EC in Eq. (3.30)?
11. How should the orders of the predictor and corrector equations be related? Why?
12. Is there a reason to use a predictor as accurate as the corrector?
13. What are the advantages and disadvantages of the PC methods compared to the RK
methods?
14. What is the reason one may want to use an implicit method?
15. What is the signicance of Requirements (i) and (ii) found before Eq. (3.47)?
16. Make sure you can explain the derivations of (3.48) and (3.49).