MATH661

Lecture 22: Ordinary Differential Equations - Multistep Methods

1.Adams-Bashforth and Adams-Moulton schemes

Consider now the approximation of $ℱ = f$ in the first-order differential equation

y^{'} = f (t, y) .

(1)

Integration over a time step $[t_{i}, t_{i + 1}]$ gives

y (t_{i + 1}) - y (t_{i}) = \int_{t_{i}}^{t_{i + 1}} f (t, y (t)) d t,

and use of quadrature formulas leads to numerical solutions for solving (1). Consider for instance data $𝒟 = {(t_{i + 1 - k}, f_{i + 1 - k}), k = 1, \dots, s}$ going back $s$ intervals of size $h$ , $t_{i + 1 - k} = t_{i + 1} - k h$ . Any quadrature formula based on this data could be used, but the most often encountered approach is to use a polynomial approximant. This can be stated in Lagrange form as

f (t, y (t)) ≅ \sum_{k = 1}^{s} ℓ_{k} (t) f_{k}, f_{k} = f (t_{i + 1 - k}, y (t_{i + 1 - k})) ≅ f (t_{i + 1 - k}, y_{i + 1 - k}) .

The last approximate equality arises from replacing the exact value $y (t_{i + 1 - k})$ by its approximation $y_{k} ≅ y (t_{i + 1 - k})$ . The result is known as an Adams-Bashforth scheme

y_{i + 1} = y_{i} + \int_{t_{i}}^{t_{i + 1}} \sum_{k = 1}^{s} ℓ_{k} (t) f_{k} d t = y_{i} + \sum_{k = 1}^{s} (\int_{t_{i}}^{t_{i + 1}} ℓ_{k} (t) d t) f_{k} = y_{i} + h \sum_{k = 1}^{s} b_{k} f_{i + 1 - k},

with coefficients that are readily computed (cf. Table 1).

b_{k} = \frac{1}{h} (\int_{t_{i}}^{t_{i + 1}} ℓ_{k} (t) d t) .

\begin{array}{lllll} s & b_{1} & b_{2} & b_{3} & b_{4} \\ 1 & 1 \\ 2 & - \\ 3 & - \\ 4 & - & - \end{array}

Table 1. Adams-Bashforth scheme coefficients.

The $s = 1$ Adams-Bashforth scheme is identical to forward Euler and the above approach yields schemes that are explicit, i.e., the new value is directly obtained from knowledge of previous values.

Choosing data $𝒟 = {(t_{i + 1 - k}, f_{i + 1 - k}), k = 0, \dots, s - 1}$ that contains the point yet to be computed $(t_{i + 1}, y_{i + 1})$ gives rise to a class of implicit schemes known as the Adams-Moulton schemes (Table 2)

y_{i + 1} = y_{i} + h \sum_{k = 0}^{s - 1} b_{k} f_{i + 1 - k},

\begin{array}{lllll} s & b_{0} & b_{1} & b_{2} & b_{3} \\ 1 & 1 \\ 2 \\ 3 & - \\ 4 & - \end{array}

Table 2. Adams-Moulton scheme coefficients.

2.Simultaneous operator approximation - linear multistep methods

Approximation of both operators $ℒ = d / d t$ and $ℱ = f$ arising in $ℒ y = ℱ y$ , or $y^{'} = f (t, y)$ is possible. Combining previous computations, the resulting schemes can be stated as

\sum_{k = 0}^{s} a_{k} y_{i + k} = h \sum_{k = 0}^{s} b_{k} f_{i + k}, f_{i + k} = f (t_{i + k}, y_{i + k}) .

(2)

Both sides arise from linear approximants: of the derivative on the left, and of $f$ on the right.

3.Consistency, convergence, stability

Any of the above schemes defines a sequence ${y_{n}}_{n \in ℕ}$ that approximates the solution $y (t_{n})$ of the initial value problem

y^{'} = f (t, y), y (0) = y_{0},

over a time interval $[0, T]$ , $t_{n} = n h$ , $h = T / N$ . A scheme is said to be convergent if

{lim}_{\begin{array}{l} h \to 0 \\ N h = T \end{array}} y_{N} = y (T) .

The above states that in the limit of taking small step sizes while maintaining $N h = T$ for some finite time $T$ , the estimate at the endpoint converges to the exact value. Such a definition is rather difficult to apply directly, and an alternative characterization of convergence is desirable.

3.1.Model problem

To motivate the overall approach, consider first the following model problem

y^{'} = λ y, y (0) = y_{0}, λ ⩽ 0

(3)

The model problem arises from truncation of the general non-linear function $f$ to first order

y^{'} = f (y) = f (0) + f^{'} (0) y + \dots .

Since $f (0)$ is a constant that simply leads to linear growth, and the model problem captures the lowest-order non-trivial behavior. The exact solution is

y (t) = e^{λ t} y_{0} \Rightarrow y (t_{n}) = e^{n λ h} y_{0}^{},

giving $y (T) = e^{λ T} y_{0}$ . The restriction of $λ ⩽ 0$ in the model problem arises from consideration of the effect of a small perturbation in the initial condition representative of floating point representation errors. This leads to $\tilde{y} (T) = e^{λ T} (y_{0} + δ)$ , and the error $ε = \tilde{y} (T) - y (T) = e^{λ T} δ$ can only be maintained small if $λ ⩽ 0$ .

Applying the forward Euler scheme to the model problem (3) gives

y_{n + 1} = y_{n} + λ h y_{n} = (1 + z) y_{n},

with $z = λ h$ . After $N$ steps the numerical approximation is

y_{N} = {(1 + z)}^{N} y_{0} .

The exponential decay of the exact solution can only be recovered if which leads to a restriction on the allowable step size

- \frac{2}{λ} > h > 0 .

If the step size is too large, $h > - 2 / λ$ , inherent floating point errors are amplified by the forward Euler method, and the scheme is said to be unstable. This is avoided by choosing a subunitary parameter $z$ , $| z | = | λ h | ⩽ 1$ , which leads to a step size restriction $h < 1 / | λ |$ .

These observations on the simple case of the Euler forward method generalize to linear multistep methods. Applying (2) to the model problem (3) leads to the following linear finite difference equation

\sum_{k = 0}^{s} a_{k} y_{i + k} = z \sum_{k = 0}^{s} b_{k} y_{i + k} .

(4)

The above is solved using a procedure analogous to that for differential equations by hypothesizing solutions of the form

y_{n} = r^{n},

to obtain a characteristic equation

π (r; z) = ρ (r) - z σ (r) = 0,

where $ρ (r), σ (r)$ are polynomials

ρ (r) = \sum_{k = 0}^{s} a_{k} r^{k}, σ (r) = \sum_{k = 0}^{s} b_{k} r^{k} .

The above polynomials allow an operational assessment of algorithms of form (2). An algorithm (2) that recovers the ordinary differential equation (1) in the limit of $h \to 0$ is said to be consistent, which occurs if and only if

ρ (1) = 0, ρ^{'} (1) - σ (1) = 0 .

Furthermore an algorithm of form (2) that does not amplify inherent floating point errors is said to be stable, which occurs if the roots of $π (r; z)$ are subunitary in absolute value

| r_{j} | < 1, π (r_{j}; z) = 0 .

Theorem. An algorithm to solve (1) that is consistent and stable is convergent.

3.2.Boundary locus method

A convenient procedure to determine the stable range of step sizes is to consider $r$ of unit absolute value

r = e^{i θ},

and evaluate the characteristic equation

π (e^{i θ}; z) = ρ (e^{i θ}) - z σ (e^{i θ}) = 0 \Rightarrow z (θ) = \frac{ρ (e^{i θ})}{σ (e^{i θ})},

where $z (θ)$ is the boundary locus delimiting zones of stability in the complex plane (Fig 1).

A method is said to be $A$ -stable if its region of stability contains the entire left half-plane in $ℂ$ , and is said to be $L$ -stable if ${lim}_{ω \to \infty} ρ (ω e^{i θ}) / σ (ω e^{i θ}) = 0 .$