Approaches to Numerical Methods for ODEs

This note introduces numerical methods for the initial value problem (IVP) of ordinary differential equations (ODEs), organized by the approach each one comes from. The list is not exhaustive, and different approaches sometimes lead to the same method, but grouping them this way makes the underlying ideas easier to follow.

Numerical Approaches to the Initial Value Problem

Let $y$ be a function of time with $y(t)\in\mathbb{R}^N$. Given the differential equation $\frac{dy}{dt} = f(t, y)$ together with an initial value $\quad y(t_0) = y_0$, the problem of finding $y(t)$ is called the initial value problem (IVP). The approaches to solving it numerically are listed below. Time is discretized as $t_n = t_0 + nh$, and $y_{n+1}$ is computed from $y_n(=y(t_n))$. Here $h$ is the step size, often taken to be a small value such as $10^{-3}$.
Throughout, $y$ denotes the exact solution and $u$ denotes the numerical solution.

1. Approximating by Taylor Expansion

This approach approximates $y(t)$ using the Taylor expansion

$$y(t_{n+1}) = y_n + f(t_n, y_n)\underbrace{(t_{n+1} - t_n)}_{=h} + \frac{1}{2}f'(t_n, y_n)\underbrace{(t_{n+1} - t_n)^2}_{=h^2} + \cdots.$$

Dropping the higher-order terms gives an approximation of $y(t)$ in terms of the right-hand side $f$ of the differential equation and its derivatives.

• Truncating at first order gives $u_{n+1} = u_n + h f(t_n, u_n),$ which is the method known as the Euler method. Since it approximates $u_{n+1}$ linearly using only $f(t_n, u_n)$, the algorithm and its implementation are simple, but both accuracy and stability are low, as discussed later.

• Truncating at second order gives $u_{n+1} = u_n + h f(t_n, u_n) + \frac{h^2}{2}f'(t_n, u_n),$ which requires computing

  $$f'(t_n, u_n) = \frac{\partial f}{\partial t}(t_n, u_n) + \frac{\partial f}{\partial u}(t_n, u_n)\underbrace{\frac{du}{dt}}_{=f(t_n, u_n)} \quad \text{(chain rule)}.$$

  Here $Df(t_0, u_0)$ is the Jacobian matrix of $f$ with respect to $u$. Putting this together,

  $$u_{n+1} \approx u_n + h f(t_n, u_n) + \frac{h^2}{2}\left(\frac{\partial f}{\partial t}(t_n, u_n) + Df(t_n, u_n)f(t_n, u_n)\right),$$

  but computing the Jacobian is often expensive, so this is not a popular method.

2. Approximating the Integral Equation by Quadrature

Rewriting the differential equation as an integral equation gives

$$y_{n+1} = y_n + \int_{t_n}^{t_{n+1}} f(t, y(t)) dt.$$

This approach approximates the integral by a quadrature rule.

• Forward Euler method Using $$\int^b_a g(t) dt \approx (b-a)g(a)$$ gives $$u_{n+1} = u_n + h f(t_n, u_n),$$ the same formula as the Euler method obtained above from the first-order Taylor expansion.
• Backward Euler method Using $$\int^b_a g(t) dt \approx (b-a)g(b)$$ gives $u_{n+1} = u_n + h f(t_{n+1}, u_{n+1}).$ This is also a first-order approximation, and since $u_{n+1}$ appears on both sides, a nonlinear equation has to be solved.
  A method in which the unknown $u_{n+1}$ appears on the right-hand side like this is called an implicit method. A method in which the unknown does not appear on the right-hand side, such as Forward Euler, is called an explicit method. Implicit methods have better stability, but solving the nonlinear equation often makes them computationally expensive. For example, one solves $$F(u_{n+1}) = u_{n+1} - u_n - h f(t_{n+1}, u_{n+1}) = 0$$ by Newton’s method. Stability is high, but so is the cost.
• Trapezoidal rule Using the trapezoidal rule $$\int^b_a g(t) dt \approx \frac{b-a}{2}(g(a) + g(b))$$ gives $u_{n+1} = u_n + \frac{h}{2}(f(t_n, u_n) + f(t_{n+1}, u_{n+1})).$ This is a second-order approximation. Since $u_{n+1}$ appears on both sides, a nonlinear equation must be solved with Newton’s method or similar, but stability is high (though lower than Backward Euler).
• Midpoint method Using $$\int^b_a g(t) dt \approx (b-a)g\left(\frac{a+b}{2}\right)$$ gives $u_{n+1} = u_n + h f\left(t_n + \frac{h}{2}, \underbrace{u(t_n + \frac{h}{2})}_{?}\right).$ Since $u(t_n + \frac{h}{2})$ is unknown, it is replaced by something like $\frac{u_n + u_{n+1}}{2}$. This yields $$u_{n+1} = u_n + h f\left(t_n + \frac{h}{2}, \frac{u_n + u_{n+1}}{2}\right),$$ which again requires solving a nonlinear equation, but has high stability (though lower than Backward Euler).

3. Predictor-Corrector Methods

The methods above require solving a nonlinear equation because $u_{n+1}$ also appears on the right-hand side. Predictor-corrector methods remove that need by replacing the $u_{n+1}$ on the right-hand side with an approximation, such as an Euler step, built from known values like $u_n$. Some examples follow.

Midpoint Method with an Euler Predictor

Approximating the midpoint with the Euler method,

$$u_{n+\frac{1}{2}} = u_n + \frac{h}{2}f(t_n, u_n)$$

makes the update formula

$$u_{n+1} = u_n + h f\left(t_n + \frac{h}{2}, u_{n+\frac{1}{2}}\right)$$

explicit.

Trapezoidal Rule with an Euler Predictor

Approximating the $u_{n+1}$ on the right-hand side with the Euler method,

$$u^*_{n+1} = u_n + h f(t_n, u_n)$$

makes the update formula

$$u_{n+1} = u_n + \frac{h}{2}(f(t_n, u_n) + f(t_{n+1}, u^*_{n+1}))$$

explicit.

Four-Stage Fourth-Order Runge-Kutta Method

Starting from Simpson’s rule

$$\int^b_a g(t) dt \approx \frac{b-a}{6}(g(a) + 4g\left(\frac{a+b}{2}\right) + g(b)),$$

both step $n+1/2$ and step $n+1$ are approximated by predictors to make the method explicit. The update formula is

$$\begin{aligned} k_1 &= f(t_n, u_n) \\ k_2 &= f\left(t_n + \frac{h}{2}, u_n + \frac{h}{2}k_1\right) \quad \left(\approx f\left(u\left(t_n + \frac{h}{2}\right)\right)\right)\\ k_3 &= f\left(t_n + \frac{h}{2}, u_n + \frac{h}{2}k_2\right) \quad \left(\approx f\left(u\left(t_n + \frac{h}{2}\right)\right)\right)\\ k_4 &= f(t_n + h, u_n + hk_3) \quad \left(\approx f\left(u\left(t_n + h\right)\right)\right)\\ u_{n+1} &= u_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4). \end{aligned}$$

This is a fourth-order approximation with good stability, but it costs more than the lower-order methods above.

4. Approximating by Polynomial Interpolation at Mesh Points

So far, every method computed $u_{n+1}$ from $u_n$ alone (raising the order by interpolating intermediate points). There are also multistep methods that use earlier points such as $u_{n-1}, u_{n-2},\cdots$. This section covers the Adams methods and BDF, the representative linear multistep methods derived from polynomial interpolation.

Adams Methods

Adams methods build a polynomial $p(t)$ interpolating the values $f_n, f_{n-1}, \cdots$ of the right-hand side $f$ at $t_n, t_{n-1}, \cdots$, and integrates it as

$$y_{n+1} = y_n + \int_{t_n}^{t_{n+1}} p(t) dt.$$

The explicit version is the Adams-Bashforth method, and the implicit version, which also uses $f_{n+1}$, is the Adams-Moulton method.

• The Adams-Bashforth method uses the Newton interpolation polynomial. First, the divided differences are defined as $$\begin{aligned} f[t_n, t_{n-1}] &= \frac{f_n - f_{n-1}}{t_n - t_{n-1}} \\ f[t_n, t_{n-1}, t_{n-2}] &= \frac{f[t_n, t_{n-1}] - f[t_{n-1}, t_{n-2}]}{t_n - t_{n-2}}. \end{aligned}$$ That is, the divided difference $f[t_n, t_{n-1}]$ is (change in $f$) / (change in $t$), and $f[t_n, t_{n-1}, t_{n-2}]$ is the divided difference of $f[t_n, t_{n-1}]$ and $f[t_{n-1}, t_{n-2}]$. (Geometrically, these correspond to slope and curvature.) Using these, the Newton interpolation polynomial of $f$ is¹ $$\begin{aligned} p_k(t) = f_n + f[t_n, t_{n-1}](t - t_n) + f[t_n, t_{n-1}, t_{n-2}](t - t_n)(t - t_{n-1}) + \cdots\\ + f[t_n, t_{n-1}, \cdots, t_{n-k+1}](t - t_n)(t - t_{n-1})\cdots(t - t_{n-k+1}). \end{aligned}$$ Integrating this, using the fact that the step size is a constant $h$, gives the formula for $u_{n+1}$. For example:
  • Two-step ($k=2$) $$u_{n+1} = u_n + \frac{h}{2}(3f_n - f_{n-1}).$$
  • Three-step ($k=3$) $$u_{n+1} = u_n + \frac{h}{12}(23f_n - 16f_{n-1} + 5f_{n-2}).$$
• The Adams-Moulton method is the implicit method that uses $f_{n+1}$. Its polynomial is $$\begin{aligned} p_k(t) = f_{n+1} + f[t_{n+1}, t_n](t - t_{n+1}) + f[t_{n+1}, t_n, t_{n-1}](t - t_{n+1})(t - t_n) + \cdots\\ + f[t_{n+1}, t_n, \cdots, t_{n-k+1}](t - t_{n+1})(t - t_n)\cdots(t - t_{n-k+2}). \end{aligned}$$ Integrating this with step size $h$ gives the formula for $u_{n+1}$. For example:
  • Two-step ($k=2$) $$u_{n+1} = u_n + \frac{h}{12}(5f_{n+1} + 8f_n - f_{n-1}).$$
  • Three-step ($k=3$) $$u_{n+1} = u_n + \frac{h}{24}(9f_{n+1} + 19f_n - 5f_{n-1} + f_{n-2}).$$

Backward Differentiation Formula (BDF)

This is also polynomial interpolation, but it interpolates the numerical solution $u(t)$ rather than the right-hand side $f$ of the differential equation (and includes $u_{n+1}$ among the nodes). The interpolating polynomial $p(t)$ is differentiated and evaluated at $t_{n+1}$, giving the equation

$$p'(t_{n+1}) = f(t_{n+1}, u_{n+1}),$$

which is solved for $u_{n+1}$. For the two-step method BDF2 ($k=2$), for instance, one builds the polynomial $p(t)$ through $u_{n+1}, u_n, u_{n-1}$ and differentiates it to get

$$p'(t_{n+1}) = \frac{3u_{n+1} - 4u_n + u_{n-1}}{2h}.$$

Therefore $u_{n+1}$ is obtained by solving

$$\frac{3u_{n+1} - 4u_n + u_{n-1}}{2h} = f(t_{n+1}, u_{n+1}).$$

BDF methods have high stability, but require solving a nonlinear equation.

5. Method of Undetermined Coefficients

Linear multistep methods (LMMs), including the Adams family and BDF, are written in the general form

$$u_{n+1} + \alpha_0 u_n + \cdots + \alpha_{k-1} u_{n-k+1} = h (\beta_{-1} f_{n+1} + \beta_0 f_n + \cdots + \beta_{k-1} f_{n-k+1}),$$

where $\alpha_i, \beta_i$ are coefficients specific to each method. Approach 4 obtained the coefficients by integrating or differentiating an interpolating polynomial. The method of undetermined coefficients instead assumes this form from the start and determines $\alpha_i, \beta_i$ by expanding $u(t)$ in a Taylor series about $t_n$.

Classifying the Methods

The approaches above can be grouped from several points of view.

1. Explicit and implicit:
   • Explicit: methods where $u_{n+1}$ does not appear on the right-hand side of $u_{n+1}= \ldots$ are cheap to compute but less stable.
   • Implicit: methods that require solving an equation with Newton’s method or similar at every step are expensive but more stable.
2. One-step and multistep:
   • One-step: $u_{n+1}$ is computed from $u_n$ alone. Stability is high.
   • Multistep: $u_{n+1}$ is computed using earlier points such as $u_{n-1}, u_{n-2}, \cdots$ as well. Efficient, but less stable. As an aside, multistep methods need a one-step method to produce the starting values for the first few steps.
3. Linear and nonlinear:
   • Linear: $u_{n+1}$ is expressed as a linear combination of past values of $u$. This includes linear multistep methods such as the Adams methods and BDF.
   • Nonlinear: methods that use predictors, such as Runge-Kutta, feed values of $f$ back into the arguments of $f$, and are therefore nonlinear in general.

References

• Maria Cameron, Numerical Methods for Solving Ordinary Differential Equations, AMSC 661 Lecture Notes, University of Maryland.

Footnotes

1. The Lagrange interpolation polynomial is equivalent, but as the number of points $k$ grows, the Lagrange form requires updating every coefficient, whereas the Newton form only needs a new higher-order term appended to the existing polynomial. ↩