1.4 Predictor-corrector methods

Predictor-corrector methods combine explicit and implicit multistep formulas to solve ordinary differential equations. The explicit formula generates a rough estimate (the prediction), and the implicit formula refines it (the correction). This two-phase approach gives you better accuracy than a purely explicit method while avoiding the full cost of solving the nonlinear system that a purely implicit method would require.

This section covers the main predictor-corrector families, how the predict-correct cycle works, and the stability and convergence theory behind them.

Predictor-corrector method overview

The core idea is straightforward: use an explicit multistep method to predict $y_{n+1}$, then feed that prediction into an implicit multistep method to correct it. Because the implicit formula already has a good initial guess from the predictor, you typically need only one or two iterations of the corrector rather than solving a full nonlinear system.

This makes predictor-corrector methods a natural bridge between explicit methods (cheap but less stable) and implicit methods (stable but expensive). They're a fundamental tool in Numerical Analysis II because they illustrate how you can combine methods strategically to get the best of both worlds.

Types of predictor-corrector methods

Adams-Bashforth-Moulton methods

The most widely used predictor-corrector family pairs an explicit Adams-Bashforth formula (predictor) with an implicit Adams-Moulton formula (corrector). Both formulas use previously computed solution values and derivative evaluations to estimate the next point.

• The predictor and corrector are typically chosen to be the same order, or the corrector is one order higher
• Only one new function evaluation is needed per step (at the predicted value), making these efficient
• Variable-order implementations exist that adapt the polynomial degree during integration for automatic error control
• Well-suited for initial value problems where the solution is reasonably smooth

Milne-Simpson methods

This pairing uses Milne's method as the predictor and Simpson's rule as the corrector. Both are based on integrating an interpolating polynomial through known function values.

• Simpson's corrector achieves fourth-order accuracy, which is attractive for its relatively low cost
• The downside: Milne-Simpson is weakly unstable for certain problems. Parasitic solutions can grow over long integrations, even when the step size is within the absolute stability region. This is a classic example of why stability analysis matters beyond just looking at order of accuracy.
• Requires more starting values than Adams methods of comparable order

Hamming's method

Hamming modified the Milne-Simpson approach specifically to fix its stability weakness.

• Introduces a weighted combination of predicted and corrected values that damps the parasitic solution Milne-Simpson suffers from
• Sacrifices a small amount of accuracy compared to Milne-Simpson in exchange for significantly better stability
• Provides a built-in error estimate from the difference between predicted and corrected values

Predictor step

Explicit methods

The predictor uses an explicit formula, meaning $y_{n+1}^p$ depends only on already-known quantities. The general Adams-Bashforth predictor takes the form:

$y_{n+1}^p = y_n + h\sum_{i=0}^{k-1} b_i f_{n-i}$

Here $f_{n-i} = f(t_{n-i}, y_{n-i})$ are derivative evaluations at previous steps, and the coefficients $b_i$ are determined by requiring the formula to be exact for polynomials up to a certain degree. For example, the fourth-order Adams-Bashforth predictor ($k=4$) has coefficients $b_0 = 55/24$, $b_1 = -59/24$, $b_2 = 37/24$, $b_3 = -9/24$.

The prediction is cheap to compute since no iteration is needed, but it's less accurate than the corrected value that follows.

Extrapolation techniques

Beyond standard Adams-Bashforth formulas, you can predict $y_{n+1}$ by extrapolating from known solution values using polynomial interpolation (Lagrange or Newton forms). Richardson extrapolation can further improve accuracy by combining predictions at different step sizes to cancel leading error terms. These techniques are more specialized and computationally expensive, but useful when higher accuracy is needed from the predictor.

Corrector step

Implicit methods

The corrector uses an implicit formula where the unknown $y_{n+1}$ appears on both sides. The general Adams-Moulton corrector is:

$y_{n+1}^c = y_n + h\sum_{i=-1}^{k-1} a_i f_{n-i}$

The $i = -1$ term means the sum includes $f_{n+1} = f(t_{n+1}, y_{n+1}^c)$, which is what makes this implicit. The key insight: instead of solving this nonlinear equation from scratch, you substitute the predicted value $y_{n+1}^p$ to evaluate $f_{n+1}$, giving you a corrected value directly.

Iteration process

The standard approach is called PECE (Predict-Evaluate-Correct-Evaluate):

1. P: Compute $y_{n+1}^p$ using the predictor formula
2. E: Evaluate $f_{n+1}^p = f(t_{n+1}, y_{n+1}^p)$
3. C: Compute $y_{n+1}^c$ using the corrector formula with $f_{n+1}^p$
4. E: Evaluate $f_{n+1}^c = f(t_{n+1}, y_{n+1}^c)$ for use in the next step

You can iterate the corrector (PECECE...) for additional refinement, using fixed-point iteration. Convergence is guaranteed when $h|{\partial f}/{\partial y}|$ is sufficiently small. In practice, one correction is often enough because the predictor already provides a good starting guess, and iterating further doesn't improve the method's order of accuracy.

Order of accuracy

Local truncation error

The local truncation error (LTE) is the error introduced in a single step, assuming all previous values are exact. For a $p$-th order method, the LTE is:

$\text{LTE} = O(h^{p+1})$

You can derive this by expanding the true solution in a Taylor series and comparing it to the numerical formula. The LTE of the predictor and corrector are typically different. Their difference provides a convenient, cheap error estimate:

$\text{Error estimate} \approx C(y_{n+1}^c - y_{n+1}^p)$

where $C$ is a known constant depending on the specific predictor-corrector pair.

Global truncation error

The global truncation error accumulates local errors over the entire integration interval $[t_0, t_f]$. For a $p$-th order method:

$\text{Global error} = O(h^p)$

This is one order lower than the local error because errors accumulate over roughly $(t_f - t_0)/h$ steps. The global error is what ultimately determines how accurate your final answer is.

Stability analysis

Absolute stability

Stability analysis determines whether small perturbations grow or decay as integration proceeds. You test this using the model equation $y' = \lambda y$, where $\lambda$ is a complex constant.

Substituting into the predictor-corrector pair gives a recurrence relation, and the method is absolutely stable for a given $h\lambda$ if all solutions of that recurrence remain bounded. The set of all stable $h\lambda$ values forms the stability region in the complex plane.

• Adams-Bashforth-Moulton methods have finite stability regions (they are not A-stable)
• Higher-order Adams methods have smaller stability regions, which is a practical limitation
• Milne-Simpson's stability region has a gap that allows parasitic growth, which is why Hamming's modification exists

Relative stability

Relative stability goes beyond asking "do errors stay bounded?" to asking "does the numerical solution preserve the qualitative behavior of the true solution?" For instance, if the true solution oscillates, does the numerical solution oscillate at roughly the right frequency? If the true solution decays, does the numerical solution decay at a comparable rate?

This is analyzed using the logarithmic norm of the Jacobian and is particularly important for oscillatory or multi-scale problems where you need the numerical method to track the solution's character, not just its magnitude.

Implementation considerations

Starting values

Multistep methods need $k$ previous values to begin, but an initial value problem only provides $y_0$. You need to generate $y_1, y_2, \ldots, y_{k-1}$ using another method. The standard approach:

1. Use a single-step method (typically a Runge-Kutta method of matching order) to compute the first $k-1$ values
2. Ensure the starting method has at least the same order of accuracy as the predictor-corrector pair
3. Once enough values are available, switch to the multistep predictor-corrector

If your starting values are less accurate than the multistep method, that lower accuracy can contaminate the entire solution.

Step size selection

Choosing the step size $h$ involves balancing accuracy against computational cost.

• Initial step size: Often estimated from the problem's timescale or by taking a trial step and checking the error
• Adaptive control: The error estimate from $y_{n+1}^c - y_{n+1}^p$ drives automatic step size adjustment. If the estimated error exceeds a tolerance, the step is rejected and retried with a smaller $h$. If the error is well below tolerance, $h$ can be increased.
• Changing step size in multistep methods is nontrivial because the formulas assume equally spaced points. You either need to re-derive coefficients for the new spacing or use the Nordsieck representation (discussed below).

Advantages vs disadvantages

Comparison with Runge-Kutta methods

Efficiency considerations

The main efficiency advantage of predictor-corrector methods is fewer function evaluations per step. Each PECE cycle requires only 2 evaluations of $f$, compared to 4 for a classical fourth-order Runge-Kutta step. For problems where evaluating $f$ is expensive (large systems, complex physics), this savings adds up.

The tradeoff is higher memory usage (you must store $k$ previous solution values and derivative evaluations) and the added complexity of startup and step size changes.

Applications in ODEs

Initial value problems

Predictor-corrector methods are most commonly applied to initial value problems of the form $y' = f(t, y)$, $y(t_0) = y_0$. They're particularly effective for problems with smooth, non-stiff solutions over long integration intervals, such as:

• Orbital mechanics (planetary motion, satellite trajectories)
• Chemical kinetics with moderate timescale separation
• Population dynamics and compartmental models in epidemiology

Boundary value problems

Predictor-corrector methods can be adapted for boundary value problems through shooting methods: you convert the BVP into a sequence of IVPs by guessing unknown initial conditions, solving forward with a predictor-corrector integrator, and iterating on the guess until the boundary conditions are satisfied. This approach is common in heat transfer, fluid dynamics, and structural analysis problems.

Error estimation techniques

Richardson extrapolation

Richardson extrapolation estimates the error by solving the same problem with two different step sizes ($h$ and $h/2$) and combining the results. If the method has order $p$, the leading error term cancels:

$y_{\text{improved}} = \frac{2^p \, y_{h/2} - y_h}{2^p - 1}$

This gives both a more accurate solution and a reliable error estimate, at the cost of roughly doubling the computational work.

Embedded formulas

Predictor-corrector methods have a natural embedded error estimate: the difference $y_{n+1}^c - y_{n+1}^p$ between the corrected and predicted values. Since the corrector is typically one order higher than the predictor, this difference approximates the local truncation error of the predictor (up to a known constant). This estimate comes essentially for free and is what drives adaptive step size control.

Adaptive step size control

Error per step vs error per unit step

Two common strategies for controlling the local error:

• Error per step (EPS): require $|\text{LTE}| \leq \text{tol}$. This bounds the error at each step but can over-constrain the problem when step sizes vary widely.
• Error per unit step (EPUS): require $|\text{LTE}|/h \leq \text{tol}$. This normalizes by step size, giving more uniform accuracy per unit of the independent variable. EPUS is generally preferred for problems with varying timescales.

Step size adjustment algorithms

Once you have an error estimate $e_n$ at step $n$, the new step size is typically chosen as:

$h_{\text{new}} = h_{\text{old}} \cdot \left(\frac{\text{tol}}{|e_n|}\right)^{1/(p+1)}$

A safety factor (e.g., 0.8 or 0.9) is applied to avoid repeatedly overshooting the tolerance. More advanced controllers like PI (Proportional-Integral) control also use the error from the previous step to smooth out step size oscillations:

$h_{\text{new}} = h_{\text{old}} \cdot \left(\frac{\text{tol}}{|e_n|}\right)^{k_I} \cdot \left(\frac{|e_{n-1}|}{\text{tol}}\right)^{k_P}$

where $k_I$ and $k_P$ are tuning parameters. This prevents the "step size hunting" that simple controllers can exhibit.

Multistep methods connection

Relationship to linear multistep methods

Predictor-corrector methods are specific implementations of the general linear multistep method (LMM):

$\sum_{i=0}^k \alpha_i y_{n+i} = h \sum_{i=0}^k \beta_i f_{n+i}$

The predictor is an LMM with $\beta_k = 0$ (explicit), and the corrector is an LMM with $\beta_k \neq 0$ (implicit). Adams methods are the special case where $\alpha_0 = -1$, $\alpha_{k} = 1$, and all other $\alpha_i = 0$, meaning only $y_n$ and $y_{n+k}$ appear on the left side.

Nordsieck form

The Nordsieck vector stores the solution and its scaled derivatives at the current point:

$z_n = \left[y_n, \; h y_n', \; \frac{h^2}{2!} y_n'', \; \ldots, \; \frac{h^q}{q!} y_n^{(q)}\right]^T$

Advancing one step becomes a matrix-vector multiply, and changing the step size only requires rescaling the vector entries. This makes variable step size and variable order implementations much cleaner than reformulating the Adams coefficients directly. The Nordsieck form is the basis for many production ODE solvers.

Convergence analysis

Consistency conditions

A method is consistent if it reproduces the ODE exactly as $h \to 0$. For a linear multistep method, consistency requires two conditions on the coefficients:

• $\sum_{i=0}^k \alpha_i = 0$ (the method sums to zero when applied to a constant)
• $\sum_{i=0}^k i \, \alpha_i = \sum_{i=0}^k \beta_i$ (the method is exact for $y = t$)

Higher-order consistency (exactness for $y = t^2, t^3, \ldots$) gives additional equations that determine the coefficients for higher-order methods.

Zero-stability requirements

Zero-stability ensures the method doesn't amplify errors even for the trivial equation $y' = 0$. It's analyzed through the first characteristic polynomial $\rho(\zeta) = \sum_{i=0}^k \alpha_i \zeta^i$:

• All roots of $\rho(\zeta)$ must satisfy $|\zeta| \leq 1$
• Any root with $|\zeta| = 1$ must be simple (not repeated)

This is the root condition. The Dahlquist equivalence theorem ties everything together: a consistent, zero-stable linear multistep method is convergent. Consistency alone is not enough, and zero-stability alone is not enough. You need both.

Milne-Simpson satisfies the root condition (it is zero-stable), but its stability region has a problematic structure that allows parasitic solutions to grow for certain $h\lambda$ values. This is a different issue from zero-stability and illustrates why absolute stability analysis is also necessary in practice.