6: Predictor–corrector methods.#
The trouble with implicit methods.#
The numerical methods for solving ODEs that we have considered so far have all been explicit, which is to say that the formula defining the \((i+1)\)st approximation can be explicitly solved for \(y_{i+1}\). Implicit methods are characterized by formula which are not able to solved in the general case. Suppose as usual that we are trying to numerically solve the first order IVP
Linearization at the \(i\)th mesh point \(x_i\) and projecting forward along the tangent line yields the (forward) Euler method
which we have already discussed in detail. On the other hand, linearization at the \((i+1)\)st mesh point \(x_{i+1}\) and projecting backward along the tangent line yields the backward Euler method
Since the precise form of the “right-hand side” \(f(x,y)\) is not known in advance, we cannot explicitly solve this equation for \(y_{i+1}\) in the general case – that is without knowing \(f(x,y)\). If the function \(f(x,y)\) turns out to be linear in \(y\), then equation (19) can easily be solved for \(y\). However, there are many important applications where \(f(x,y)\) is nonlinear in \(y\). Thus, in order to use an implicit formula such as (19), it is necessary to incorporate some nonlinear equation solver such as Newton’s method or fixed-point iteration to solve for \(y_{i+1}\).
Predictor–corrector methods.#
Since fixed-point iteration techniques require an initial guess, it is common to use an explicit method to “predict” the value of \(y_{i+1}\) before employing the implicit method to “correct” the guess. Thus, the pairing of an explicit predictor method with an implicit corrector method is often referred to as a predictor-corrector method. The question of how many times to apply the fixed-point iteration of the corrector method can be a delicate matter. Each iteration brings us closer to convergence but at the cost of more evaluations of the function \(f(x,y)\). In theory, one might choose to iterate the corrector until the value \(y_{i+1}\) converges to within some tolerance. Since iteration of corrector method such as (19) only brings us closer to the solution of (19) and not necessarily the true solution \(y(x_{i+1})\), it is often not worth the additional cost to apply more than 1 iteration. If higher accuracy is needed, it is probably more efficient to reduce the step size \(h\).
Advantages and disadvantages.#
The main advantage of a predictor-corrector method is that the implicit (corrector) method tends to give the method better stability properties. A numerical method is said to be stable if small perturbations in the initial data do not cause the resulting numerical solution to diverge away from the original as \(x\to\infty\). The obvious disadvantage of a predictor-corrector method is that each iteration of the implicit method costs additional functional evaluations which translates to more work. To mitigate this effect, the number of corrector iterations per step is generally kept low.
Adams–Moulton 1-step (implicit) method.#
As usual, consider the first order IVP
Choosing a model of the form
where \(y_i'=f(x_i,y_i)\) for each \(i\ge 0\), we proceed by the method of undetermined coefficients. In particular, we force the model to be exact for the first three monomials \(y(x)=1\), \(y(x)=x\), and \(y(x)=x^2\). Thus, by a calculation that is precisely similar to the derivation that we gave for the Adams–Bashforth two-step (explicit) method, we ultimately arrive at the formula for the Adams–Moulton one-step (implicit) method, namely,
This method is often paired with the order 2 Adams–Bashforth two-step (explicit) method to create the order 2 Adams–Bashforth–Moulton (predictor-corrector) method (ABM2)
where \(y_i'=f(x_i,y_i)\) and \(\hat{y_{i+1}}' = f(x_{i+1}, \hat{y_{i+1}})\) for each \(i\ge 0\). Here \(\hat{y_{i+1}}\) is referred to as the predicted value and \(y_{i+1}\) the corrected value.
Adams–Moulton 3-step (implicit) method.#
In a similar manner, one may derive the Adams–Moulton three-step (implicit) method, which is defined by the formula
and has local truncation error that is \(O(h^5)\). It is, therefore, considered an order 4 method. This method is commonly paired with the order 4 Adams–Bashforth four-step (explicit) method to create the order 4 Adams–Bashforth–Moulton (predictor-corrector) method (ABM4)
where \(y_i'=f(x_i,y_i)\) and \(\hat{y_{i+1}}' =f(x_{i+1},\hat{y_{i+1}})\) for each \(i\ge 0\).
A Python implementation of ABM2.#
The following Python implementation of ABM2 is included in the math263 module.
import numpy as np
def abm2(f, a, b, y0, n):
"""
numerically solves the IVP
y' = f(x,y), y(a)=y0
over the interval [a, b] via n steps of second order Adams–Bashforth–Moulton
predictor-corrector method
"""
h = (b - a) / n
x = np.linspace(a, b, num=n + 1)
y = np.empty((x.size, np.size(y0)))
y[0] = y0
# starter method: Heun's modified Euler method
k1 = f(x[0], y[0])
k2 = f(x[1], y[0] + h * k1)
y[1] = y[0] + h * (k1 + k2) / 2
# continuing method: ABM2 predictor-corrector
f2 = f(x[0], y[0])
for i in range(1, n):
# predict with AB2
f1 = f(x[i], y[i])
yhat = y[i] + h * (3 * f1 - f2) / 2
# correct with AM1
f0 = f(x[i + 1], yhat)
y[i + 1] = y[i] + h * (f0 + f1) / 2
f2 = f1
return x, y
Example.#
We now compare the Adams–Bashforth 2-step method (AB2) and the Adams–Bashforth–Moulton 2-step method (ABM2) for our test IVP
over the interval \([a,b]=[0,2]\) with \(n=10\) steps.
import numpy as np
import sympy
from tabulate import tabulate
import math263
# define IVP parameters
f = lambda x, y: y - x**2 + 1
a, b = 0, 2
y0 = 1 / 2
# solve the IVP symbolically with the sympy library
x = sympy.Symbol("x")
y = sympy.Function("y")
ode = sympy.Eq(y(x).diff(x), f(x, y(x)))
soln = sympy.dsolve(ode, y(x), ics={y(a): y0})
sym_y = sympy.lambdify(x, soln.rhs, modules=["numpy"])
# numerically solve the IVP with n=10 steps of AB2 and n=10 steps of ABM2
n = 10
x, y_ab2 = math263.ab2(f, a, b, y0, n)
x, y_abm2 = math263.abm2(f, a, b, y0, n)
# tabulate the results
print(f"Comparison of global errors for AB2 and ABM2.")
table = np.c_[x, abs(sym_y(x) - y_ab2[:, 0]), abs(sym_y(x) - y_abm2[:, 0])]
hdrs = ["i", "x_i", "AB2 global error", "ABM2 global error"]
print(
tabulate(
table,
hdrs,
tablefmt="mixed_grid",
floatfmt=["0.0f", "0.2f", "0.5e", "0.5e"],
showindex=True,
)
)
Comparison of global errors for AB2 and ABM2.
┍━━━━━┯━━━━━━━┯━━━━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━━━┑
│ i │ x_i │ AB2 global error │ ABM2 global error │
┝━━━━━┿━━━━━━━┿━━━━━━━━━━━━━━━━━━━━┿━━━━━━━━━━━━━━━━━━━━━┥
│ 0 │ 0.00 │ 0.00000e+00 │ 0.00000e+00 │
├─────┼───────┼────────────────────┼─────────────────────┤
│ 1 │ 0.20 │ 3.29862e-03 │ 3.29862e-03 │
├─────┼───────┼────────────────────┼─────────────────────┤
│ 2 │ 0.40 │ 2.28765e-03 │ 4.30765e-03 │
├─────┼───────┼────────────────────┼─────────────────────┤
│ 3 │ 0.60 │ 2.00600e-04 │ 5.57120e-03 │
├─────┼───────┼────────────────────┼─────────────────────┤
│ 4 │ 0.80 │ 2.95246e-03 │ 7.18297e-03 │
├─────┼───────┼────────────────────┼─────────────────────┤
│ 5 │ 1.00 │ 7.50351e-03 │ 9.23551e-03 │
├─────┼───────┼────────────────────┼─────────────────────┤
│ 6 │ 1.20 │ 1.39116e-02 │ 1.18450e-02 │
├─────┼───────┼────────────────────┼─────────────────────┤
│ 7 │ 1.40 │ 2.27729e-02 │ 1.51575e-02 │
├─────┼───────┼────────────────────┼─────────────────────┤
│ 8 │ 1.60 │ 3.48556e-02 │ 1.93565e-02 │
├─────┼───────┼────────────────────┼─────────────────────┤
│ 9 │ 1.80 │ 5.11477e-02 │ 2.46721e-02 │
├─────┼───────┼────────────────────┼─────────────────────┤
│ 10 │ 2.00 │ 7.29153e-02 │ 3.13931e-02 │
┕━━━━━┷━━━━━━━┷━━━━━━━━━━━━━━━━━━━━┷━━━━━━━━━━━━━━━━━━━━━┙
Below we include a comparison of the global errors at the right-most endpoint of the interval. It is evident that each is an order 2 method. However, the ABM2 predictor-corrector method is more accurate for this example.
# compute abs errors at right endpoint for various step-sizes
base = 10
max_exp = 7
num_steps = [base**j for j in range(1, max_exp)]
h = [(b - a) / n for n in num_steps]
ab2_errors = [
abs(math263.ab2(f, a, b, y0, n)[1][:, 0][-1] - sym_y(b)) for n in num_steps
]
abm2_errors = [
abs(math263.abm2(f, a, b, y0, n)[1][:, 0][-1] - sym_y(b)) for n in num_steps
]
# tabulate the results
print(f"Comparison of global errors |y_n - y({b})| for various step-sizes.")
table = np.c_[h, ab2_errors, abm2_errors]
hdrs = ["step-size", "AB2 global error", "ABM2 global error"]
print(tabulate(table, hdrs, tablefmt="mixed_grid", floatfmt=["0.6f", "0.6e", "0.6e"]))
Comparison of global errors |y_n - y(2)| for various step-sizes.
┍━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━━━┑
│ step-size │ AB2 global error │ ABM2 global error │
┝━━━━━━━━━━━━━┿━━━━━━━━━━━━━━━━━━━━┿━━━━━━━━━━━━━━━━━━━━━┥
│ 0.200000 │ 7.291525e-02 │ 3.139310e-02 │
├─────────────┼────────────────────┼─────────────────────┤
│ 0.020000 │ 1.179610e-03 │ 2.536390e-04 │
├─────────────┼────────────────────┼─────────────────────┤
│ 0.002000 │ 1.226335e-05 │ 2.470403e-06 │
├─────────────┼────────────────────┼─────────────────────┤
│ 0.000200 │ 1.230993e-07 │ 2.463756e-08 │
├─────────────┼────────────────────┼─────────────────────┤
│ 0.000020 │ 1.231354e-09 │ 2.463727e-10 │
├─────────────┼────────────────────┼─────────────────────┤
│ 0.000002 │ 1.189537e-11 │ 2.408740e-12 │
┕━━━━━━━━━━━━━┷━━━━━━━━━━━━━━━━━━━━┷━━━━━━━━━━━━━━━━━━━━━┙