3 Numerical Solutions for Initial Value Problems
Learning objectives
After this chapter, you should be able to:
- distinguish local truncation error from global error;
- derive the explicit and implicit Euler updates for an initial value problem;
- explain why implicit methods require the solution of an algebraic or nonlinear system at every step; and
- relate convergence and absolute stability to the choice of a time step.
In many cases, finding an analytical solution to initial value problems (IVPs) is either impossible or impractical due to the complexity of the integrals involved. Therefore, numerical methods are commonly used to approximate solutions. These methods can be applied not only to single IVPs but also extend naturally to systems of differential equations. Additionally, since higher-order differential equations can be reformulated as systems of first-order equations, we will focus on numerical methods for first-order systems.
Definition 3.1 (Time Step) On a time interval I = [t_0, t_0 + T], we partition the interval into n subintervals, known as time steps. Each subinterval
I_k = [t_{k-1}, t_k]
has a step size \Delta t_k = t_k - t_{k-1}. The total partitioning into n steps ensures that t_n = t_0 + T. The term k-th time step may refer either to the subinterval I_k or the specific time point t_k, depending on the context. For simplicity, we often assume evenly spaced time steps, where \Delta t_k = h for all k.
Numerical methods for solving IVPs are divided into two primary categories: one-step methods and multi-step methods. In one-step methods, the next approximation z_{k+1} depends only on the current step z_k. In contrast, multi-step methods calculate z_{k+1} using information from several previous steps, including z_{k-1}, z_{k-2}, \dots, z_{k-n}.
Definition 3.2 (Explicit Methods) An explicit method (also called forward integration) calculates approximations sequentially. Starting with an initial condition z_0 at time t_0, it computes the sequence z_1, z_2, \dots, z_n using the following update formula:
\bm{z}_{k+1} = \bm{z}_{k} + \Delta t_k \bm{F}(t_k, \bm{z}_{k}, \dots, \bm{z}_{k-p}; \Delta t_k)
Here, \bm{F} is called the increment function, which determines the change from one time step to the next. The parameter p refers to the order of the numerical multi-step method. If p = 1, the method is a one-step method.
Definition 3.3 (Implicit Methods) An implicit method (also called backward integration) is similar in structure to explicit methods but requires solving for z_k implicitly. The update formula is given by:
\bm{z}_{k+1} = \bm{z}_k + \Delta t_k \bm{F}(t_{k+1}, \bm{z}_{k+1}, \bm{z}_{k},\dots, \bm{z}_{k-p};\Delta t_h)
In this case, the increment function \bm{F} depends on \bm{z}_{k+1}, meaning the equation must be solved iteratively or using numerical techniques to find \bm{z}_{k+1}. Implicit methods often provide better stability, especially for stiff differential equations, but are computationally more intensive because of the need to solve for \bm{z}_{k+1} at each step.
In summary, the distinction between explicit and implicit methods is primarily in how the next step is computed. Explicit methods are generally simpler but may require smaller time steps for stability, whereas implicit methods can handle larger time steps but at the cost of increased computational effort.
3.1 Euler Method
The Euler method is the simplest one-step numerical method for solving initial value problems (IVPs). One intuitive way to derive this method is through the truncated Taylor series expansion.
Consider the following IVP:
\left\{ \begin{aligned} \bm{z}'(t) &= \bm{f}(t,\bm{z}), \\[0.25em] \bm{z}(t_0) &= \bm{z}_0, \end{aligned} \right.
The Taylor series expansion of the solution \bm{z}(t) around t_0 gives us:
\bm{z}(t_0 + h) = \bm{z}(t_0) + \Delta t\bm{z}'(t_0)+ \mathcal{O}(\Delta t^2),
where \bm{z}'(t_0) = \bm{f}(t_0,\bm{z}_0). Neglecting higher-order terms, we approximate \bm{z}(t_0 + \Delta t) as:
\bm{z}(t_0 + \Delta t) \approx \bm{z}(t_0) + \Delta t \bm{z}'(t_0),
The difference between this approximation and the exact solution \bm{z}(t_0 + \Delta t) is proportional to \Delta t^2, meaning the error decreases as \Delta t becomes smaller.
General Form of the Euler Method
To generalize this for any time interval I = [t_0, t_0 + T], we divide the interval into subintervals and apply the approximation at each step. For the k-th subinterval
I_k = [t_{k-1}, t_k]
the explicit Euler method for the k-th step is given by:
\bm{z}_{k+1} = \bm{z}_k + \Delta t\bm{f}(t_k,\bm{z}_k), \quad k = 0, 1,2,\dots, n,
where \Delta t is the step size.
On the other hand, the implicit Euler method for the k-th step is defined as:
\bm{z}_{k+1} = \bm{z}_k + \Delta t\bm{f}(t_{k+1}, \bm{z}_{k+1}), \quad k = 0, 1,2,\dots, n-1.
In this case, the update depends on the unknown value \bm{z}_{k+1}, so an algebraic equation (nonlinear when \bm{f} is nonlinear) must be solved at every step.
Local and global error
Two different errors must not be confused. The local truncation error is the defect produced by one numerical step when its input is the exact solution. For a method of global order p, it is typically \mathcal{O}(\Delta t^{p+1}). The global error is the difference between the numerical and exact solutions after many steps; on a fixed time interval it is typically \mathcal{O}(\Delta t^p), provided the method propagates perturbations in a stable way. The loss of one power of \Delta t reflects the approximately T/\Delta t steps taken over an interval of length T.
Convergence and Stability
Both the explicit and implicit Euler methods are first-order accurate, meaning that their local truncation error is \mathcal{O}(\Delta t^2), leading to a global error that scales as \mathcal{O}(\Delta t). As the step size \Delta t\to 0, both methods converge to the exact solution, with the error at any time
t \in I = [t_0, t_0 + T],
bounded by C \Delta t, where C is a constant. However, in practice, this constant can grow exponentially with T, following
C \propto e^{L T},
where L is the Lipschitz constant of the function \bm{f}(t, \bm{z}). This growth in error amplification is particularly problematic for large-time integration.
Stability Considerations
One of the key differences between the explicit and implicit Euler methods lies in their stability properties. The explicit Euler method has a stability constraint dictated by the step size \Delta t, requiring \Delta t to be sufficiently small to ensure numerical stability. Specifically, for a linear test equation of the form
\frac{\mathrm{d}z}{\mathrm{d}t} = \lambda z, \quad \lambda \in \Bbb{C}
the explicit Euler method is stable only if h satisfies the stability condition:
|1 + \Delta t\lambda| \leq 1.
For stiff problems, where solutions evolve on vastly different time scales and exhibit rapid transient dynamics, this condition can impose severely restrictive time-step limitations, making explicit Euler impractical.
In contrast, the implicit Euler method is unconditionally stable for all values of h when applied to the same test equation. This property makes it particularly advantageous for stiff problems, where it allows for significantly larger time steps without numerical instability. However, this comes at the cost of solving a nonlinear equation at each step, increasing computational expense.
Thus, while both methods share the same convergence rate of \mathcal{O}(\Delta t), their stability properties determine their applicability:
- Explicit Euler: Suitable for non-stiff problems but requires small \Delta t for stability.
- Implicit Euler: More robust for stiff problems, allowing larger \Delta t, but requires solving implicit equations at each step.
3.2 Heun and Collatz Methods
Heun’s Method
Heun’s method, also known as the improved Euler method, is a second-order numerical method for solving initial value problems (IVPs). It refines the accuracy of the explicit Euler method by incorporating an additional correction step.
Consider the IVP:
\left\{ \begin{aligned} \bm{z}'(t) &= \bm{f}(t,\bm{z}), \\[0.25em] \bm{z}(t_0) &= \bm{z}_0, \end{aligned} \right.
Similar to the Euler method, Heun’s method begins by estimating the next step using the explicit Euler approximation:
\tilde{\bm{z}}_{k+1} = \bm{z}_k + \Delta t\bm{f}(t_k, \bm{z}_k).
This intermediate value \tilde{\bm{z}}_{k+1} is then used to refine the estimate by averaging the slopes at t_k and t_{k+1}:
\bm{z}_{k+1} = \bm{z}_k + \frac{\Delta t}{2} \left[ \bm{f}(t_k, \bm{z}_k) + \bm{f}(t_{k+1}, \tilde{\bm{z}}_{k+1}) \right].
This two-stage approach makes Heun’s method more accurate than explicit Euler, achieving a global error of \mathcal{O}(\Delta t^2) instead of \mathcal{O}(\Delta t).
Convergence and Stability
Heun’s method is second-order accurate, meaning that its local truncation error is \mathcal{O}(\Delta t^3), leading to a global error of \mathcal{O}(\Delta t^2). Compared to the explicit Euler method, it provides improved accuracy for the same step size.
In terms of stability, Heun’s method is conditionally stable and has a larger stability region than explicit Euler. However, it still suffers from stability constraints in stiff problems, requiring smaller time steps to maintain numerical stability.
Collatz Method
The Collatz method, is another improvement over the Euler method. It introduces multiple stages to achieve better accuracy by approximating the solution using additional intermediate steps.
Formulation of the Collatz Method
For a given step size \Delta t, the Collatz method approximates the solution as follows:
Compute an initial (half) Euler step: \tilde{\bm{z}} = \bm{z}_k + \frac{\Delta t}{2} \bm{f}\left(t_k, \bm{z}_k\right).
Compute a midpoint correction: \bm{z}_{k+1} = \bm{z}_k + \Delta t\bm{f}\left(t_k + \frac{h}{2}, \tilde{\bm{z}} \right).
This results in an overall second-order accuracy, improving both the accuracy and stability of the method over the simple Euler approach.
Convergence and Stability
Like Heun’s method, the Collatz method achieves a global error of \mathcal{O}(\Delta t^2), making it more accurate than the first-order Euler method. Its stability region is also significantly larger, making it more suitable for problems where Euler’s method struggles.
However, similar to Heun’s method, the Collatz method is still not unconditionally stable and may require small step sizes for stiff problems.
Comparison of Methods
| Method | Order of Accuracy | Stability | Advantages | Disadvantages |
|---|---|---|---|---|
| Explicit Euler | \mathcal{O}(\Delta t) | Conditionally stable | Simple and fast | Requires small \Delta t for stability |
| Heun’s Method | \mathcal{O}(\Delta t^2) | Improved but not unconditionally stable | More accurate than Euler | Still requires small h for stiff problems |
| Collatz Method | \mathcal{O}(\Delta t^2) | Improved but not unconditionally stable | Better than Heun for some problems | Still requires small h for stiff problems |
Both Heun’s and Collatz methods improve upon the Euler method by incorporating corrections that enhance accuracy and stability, making them more suitable for practical numerical integration.