5: Multistep methods.

The methods that we have considered so far are called single-step methods because the point \((x_{i+1}, y_{i+1})\) is computed using only information from the previous step, namely \((x_i, y_i)\). In particular, the methods have all depended on a formula of the form

\[\begin{equation*} y_{i+1} = y_i + h\phi(x_i, y_i, h). \end{equation*}\]

These methods are also called starting methods because the initial-value condition \(y(x_0)=y_0\) gives enough information to use these methods from the start. On the other hand, multistep methods use more than one of the previously computed \(y_j\). An \(m\)-step method would use the previous \(m\) values: \((x_i,y_i), (x_{i-1},y_{i-1}),\dots, (x_{i-m+1},y_{i-m+1})\). A linear \(m\)-step method takes the form

(14)\[y_{i+1} = \sum_{j=1}^m\alpha_jy_{i+1-j} + h\sum_{j=0}^m\beta_jf(x_{i+1-j}, y_{i+1-j}).\]

Clearly such a method cannot be used until \(m\) previous values are known. For this reason, these methods must be paired with a starting method and are themselves called continuing methods. However, there is an additional wrinkle if \(\beta_0\ne 0\). In particular, if \(\beta_0\ne 0\), then \(y_{i+1}\) appears on both sides of the defining equation. Therefore, we say that the method is explicit if \(\beta_0=0\) and implicit otherwise. If the method is implicit and the function \(f(x, y)\) is linear in \(y\), then it is easy to solve for \(y_{i+1}\) once the specific IVP is given, but in practice this is often not the case. We will consider implicit methods in the next lecture. For now, we concentrate on the case of explicit methods. Some of the most popular explicit linear multistep methods are the Adams–Bashforth methods, which set \(\beta_0=0\) and \(\alpha_j=0\) for \(j\ge 2\).

Adams–Bashforth 2-step method.

For a \(2\)-step Adams—Bashforth method, equation (14) simplifies to

(15)\[y_{i+1} = \alpha_1y_i + h\big(\beta_1f(x_i,y_i) + \beta_2f(x_{i-1},y_{i-1})\big).\]

The idea then is to force this formula to be exact for the first 3 terms of the Taylor series of the true solution \(y\). This will ensure that the local truncation error for the method is \(O(h^3)\) as \(h\to 0\). Note that with only \(3\) parameters to determine (viz, \(\alpha_1, \beta_1, \beta_2)\), this is the best that can be done. By linearity, it is enough to make it exact for the cases \(y(x)=1\), \(y(x)=x\), and \(y(x)=x^2\). If \(y(x)=1\), then \(f(x,y(x)) = y'(x) = 0\), and so the equation (15) becomes

(16)\[1=\alpha_1.\]

Similarly, if \(y(x)=x\), then \(y'(x)=1\), and (15) becomes

(17)\[x_{i+1}=\alpha_1 x_i + h(\beta_1+\beta_2).\]

Finally, if \(y(x)=x^2\), then \(y'(x)=2x\), and (15) becomes

(18)\[x_{i+1}^2 = \alpha_1x_i^2+ h(2\beta_1x_i + 2\beta_2x_{i-1}).\]

Now, equations (16), (17), and (18) must hold for all values of \(x_i\) and \(h\). Choosing \(x_{i}=1\) and \(h=1\) so that \(x_{i-1}=0\) and \(x_{i+1}=2\), we arrive at the system of equations

\[\begin{align*} 1&=\alpha_1,\\ 2&=\alpha_1+\beta_1+\beta_2,\\ 4&=\alpha_1+2\beta_1. \end{align*}\]

Solving the system yields \(\alpha_1=1\), \(\beta_1=3/2\), and \(\beta_2=-1/2\). The Adams–Bashforth two-step method (AB2) is therefore defined by the recurrence

\[\begin{equation*} y_{i+1} = y_i + \frac{h}{2}\big(3y_i' - y_{i-1}'\big), \end{equation*}\]

where \(y_i' = f(x_i, y_i)\) for each \(i\ge 1\).

As we noted above, the local truncation error for the method is \(O(h^3)\) as \(h\to 0\). This makes AB2 an order \(2\) method. Thus, it should be paired with a starter method whose order of accuracy is at least \(2\) such as Heun’s modified Euler method.

Adams–Bashforth 4-step method.

In a similar manner, one may derive the Adams–Bashforth four-step method (AB4), which is defined by the recurrence

\[\begin{equation*} y_{i+1}=y_i+\frac{h}{24}\big(55y_i'-59y_{i-1}'+37y_{i-2}'-9y_{i-3}'\big) \end{equation*}\]

and has local truncation error that is \(O(h^5)\) as \(h\to 0\). Thus, it should be paired with an order \(4\) starting method such as the classical RK4 method.

Advantages and disadvantages.

The disadvantage of a multistep method is that it is typically slightly more difficult to implement than a single-step method. It even requires a single-step method to get started. However, there tends to be a major gain in order of accuracy per work. We usually measure the work per step of a numerical method based on the number of times that we need to evaluate the “right-hand side” function \(f(x,y)\). The typical single-step Runge-Kutta method of order 2 (e.g., Heun’s modified Euler method) requires 2 functional evaluations per step. We compare this to the Adams–Bashforth 2-step method which has the same order of accuracy and only requires one additional functional evaluation per step since it saves the evaluations from previous steps until it is finished with them. Furthermore, since the error \(|y(x_i)-y_i|\) tends to increase with \(i\), it is reasonable to hope that a multistep may have better accuracy (not just order of accuracy) by mere fact that it explicitly incorporates more accurate previous data when approximating \(y_{i+1}\).

Python implementation.

Our Python implementation of the Adams–Bashforth 2-step method uses Heun’s modified Euler method as starter. The code is included in the math263 module.

import numpy as np


def ab2(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 method
    """
    h = (b - a) / n
    x = np.linspace(a, b, num=n + 1)
    y = np.empty((x.size, np.size(y0)))
    y[0] = y0
    # take first step with Heun's MEM
    k1 = f(x[0], y[0])
    k2 = f(x[1], y[0] + h * k1)
    y[1] = y[0] + h * (k1 + k2) / 2
    # begin multistepping
    f2 = f(x[0], y[0])
    for i in range(1, n):
        f1 = f(x[i], y[i])
        y[i + 1] = y[i] + h * (3 * f1 - f2) / 2
        f2 = f1
        # step f-vals down to get ready for next step
    return x, y

Example.

We now conduct an empirical comparison of AB2 with Heun’s modified Euler method (MEM) using the IVP

\[\begin{align*} y'&=y-x^2+1,\\ y(0)&=1/2 \end{align*}\]

as a test case. First we use both methods to solve the problem over the interval \([0,2]\) in \(n=10\) steps. We then present the global truncation errors \(e_i=|y(x_i) - y_i|\) for each method at each step in a table.

from timeit import timeit

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 forward Euler and n=10 steps of RK4
n = 10
x, y_mem = math263.mem(f, a, b, y0, n)
x, y_ab2 = math263.ab2(f, a, b, y0, n)

# tabulate the results
print(
    f"Comparison of global errors for MEM and AB2\
 across interval for step-size h = {(b - a)/n}."
)
table = np.c_[x, abs(sym_y(x) - y_mem[:, 0]), abs(sym_y(x) - y_ab2[:, 0])]
hdrs = ["i", "x_i", "MEM global error", "AB2 global error"]
print(tabulate(table, hdrs, tablefmt="mixed_grid", floatfmt="0.5f", showindex=True))

Comparison of global errors for MEM and AB2 across interval for step-size h = 0.2.
┍━━━━━┯━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━━┑
│   i │     x_i │   MEM global error │   AB2 global error │
┝━━━━━┿━━━━━━━━━┿━━━━━━━━━━━━━━━━━━━━┿━━━━━━━━━━━━━━━━━━━━┥
│   0 │ 0.00000 │            0.00000 │            0.00000 │
├─────┼─────────┼────────────────────┼────────────────────┤
│   1 │ 0.20000 │            0.00330 │            0.00330 │
├─────┼─────────┼────────────────────┼────────────────────┤
│   2 │ 0.40000 │            0.00717 │            0.00229 │
├─────┼─────────┼────────────────────┼────────────────────┤
│   3 │ 0.60000 │            0.01170 │            0.00020 │
├─────┼─────────┼────────────────────┼────────────────────┤
│   4 │ 0.80000 │            0.01699 │            0.00295 │
├─────┼─────────┼────────────────────┼────────────────────┤
│   5 │ 1.00000 │            0.02317 │            0.00750 │
├─────┼─────────┼────────────────────┼────────────────────┤
│   6 │ 1.20000 │            0.03036 │            0.01391 │
├─────┼─────────┼────────────────────┼────────────────────┤
│   7 │ 1.40000 │            0.03871 │            0.02277 │
├─────┼─────────┼────────────────────┼────────────────────┤
│   8 │ 1.60000 │            0.04839 │            0.03486 │
├─────┼─────────┼────────────────────┼────────────────────┤
│   9 │ 1.80000 │            0.05956 │            0.05115 │
├─────┼─────────┼────────────────────┼────────────────────┤
│  10 │ 2.00000 │            0.07242 │            0.07292 │
┕━━━━━┷━━━━━━━━━┷━━━━━━━━━━━━━━━━━━━━┷━━━━━━━━━━━━━━━━━━━━┙

We now compare the global error for AB2 at the right endpoint of the interval with that of MEM as we shrink the step-size.

# 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]
mem_errors = [
    abs(math263.mem(f, a, b, y0, n)[1][:, 0][-1] - sym_y(b)) 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
]
# compare size of error to size at previous step-size
mem_cutdown = [mem_errors[i + 1] / mem_errors[i] for i in range(len(num_steps) - 1)]
mem_cutdown.insert(0, None)
ab2_cutdown = [ab2_errors[i + 1] / ab2_errors[i] for i in range(len(num_steps) - 1)]
ab2_cutdown.insert(0, None)

# tabulate the results
print(f"Comparison of global errors |y_n - y({b})| for various step-sizes.")
table = np.c_[h, mem_errors, mem_cutdown, ab2_errors, ab2_cutdown]
hdrs = ["step-size", "MEM global error", "cutdown", "AB2 global error", "cutdown"]
print(
    tabulate(table, hdrs, tablefmt="mixed_grid", floatfmt=["0.7f", "e", "f", "e", "f"])
)

Comparison of global errors |y_n - y(2)| for various step-sizes.
┍━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━┑
│   step-size │   MEM global error │   cutdown │   AB2 global error │   cutdown │
┝━━━━━━━━━━━━━┿━━━━━━━━━━━━━━━━━━━━┿━━━━━━━━━━━┿━━━━━━━━━━━━━━━━━━━━┿━━━━━━━━━━━┥
│   0.2000000 │       7.241732e-02 │           │       7.291525e-02 │           │
├─────────────┼────────────────────┼───────────┼────────────────────┼───────────┤
│   0.0200000 │       7.797255e-04 │  0.010767 │       1.179610e-03 │  0.016178 │
├─────────────┼────────────────────┼───────────┼────────────────────┼───────────┤
│   0.0020000 │       7.846676e-06 │  0.010063 │       1.226335e-05 │  0.010396 │
├─────────────┼────────────────────┼───────────┼────────────────────┼───────────┤
│   0.0002000 │       7.851535e-08 │  0.010006 │       1.230993e-07 │  0.010038 │
├─────────────┼────────────────────┼───────────┼────────────────────┼───────────┤
│   0.0000200 │       7.851320e-10 │  0.010000 │       1.231354e-09 │  0.010003 │
├─────────────┼────────────────────┼───────────┼────────────────────┼───────────┤
│   0.0000020 │       7.706724e-12 │  0.009816 │       1.189537e-11 │  0.009660 │
┕━━━━━━━━━━━━━┷━━━━━━━━━━━━━━━━━━━━┷━━━━━━━━━━━┷━━━━━━━━━━━━━━━━━━━━┷━━━━━━━━━━━┙

We also compare the empirical average running time for our implementations as the number of steps increases.

num_trials = 10
mem_times = [
    timeit(lambda: math263.mem(f, a, b, y0, base**j), number=num_trials) / num_trials
    for j in range(1, max_exp)
]
ab2_times = [
    timeit(lambda: math263.ab2(f, a, b, y0, base**j), number=num_trials) / num_trials
    for j in range(1, max_exp)
]

# tabulate the results
print(f"Comparison of average compute time for various step-sizes.")
table = np.c_[num_steps, mem_times, ab2_times]
hdrs = ["num steps", "MEM time (secs)", "AB2 time (secs)"]
print(tabulate(table, hdrs, tablefmt="mixed_grid", floatfmt=["0.0f", "0.5f", "0.5f"]))

Comparison of average compute time for various step-sizes.
┍━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━━━┑
│   num steps │   MEM time (secs) │   AB2 time (secs) │
┝━━━━━━━━━━━━━┿━━━━━━━━━━━━━━━━━━━┿━━━━━━━━━━━━━━━━━━━┥
│          10 │           0.00005 │           0.00004 │
├─────────────┼───────────────────┼───────────────────┤
│         100 │           0.00044 │           0.00033 │
├─────────────┼───────────────────┼───────────────────┤
│        1000 │           0.00443 │           0.00319 │
├─────────────┼───────────────────┼───────────────────┤
│       10000 │           0.04361 │           0.03188 │
├─────────────┼───────────────────┼───────────────────┤
│      100000 │           0.43594 │           0.31782 │
├─────────────┼───────────────────┼───────────────────┤
│     1000000 │           4.34698 │           3.19615 │
┕━━━━━━━━━━━━━┷━━━━━━━━━━━━━━━━━━━┷━━━━━━━━━━━━━━━━━━━┙