14.5 Numerical methods
A one-step method advances \(y_{k+1} = y_k + h\Phi (t_k, y_k, h)\) on a grid \(t_k = t_0 + kh\). Euler’s method is \(\Phi = f(t, y)\); the classical fourth-order Runge–Kutta method combines four evaluations of \(f\). The local error is the error made in one step from the exact solution, and the method has order \(p\) if the local error is \(O(h^{p+1})\).
Theorem 14.5.1. Let \(f\) be Lipschitz in \(y\) with constant \(L\) and the exact solution \(y\) be twice continuously differentiable on \([t_0, T]\) with \(|y''| \leq M\). Euler’s method with step \(h\) has global error \[ |y_k - y(t_k)| \leq \frac {Mh}{2L}\bigl (e^{L(t_k - t_0)} - 1\bigr ) , \] so it converges with order \(1\). More generally, a one-step method of order \(p\) whose increment \(\Phi \) is Lipschitz in \(y\) has global error \(O(h^p)\).
Proof. By Taylor’s theorem \(y(t_{k+1}) = y(t_k) + hf(t_k, y(t_k)) + \tau _k\) with \(|\tau _k| \leq Mh^2/2\). Subtracting from \(y_{k+1} = y_k + hf(t_k, y_k)\), the error \(e_k = y_k - y(t_k)\) satisfies \(|e_{k+1}| \leq (1 + hL)|e_k| + Mh^2/2\). By induction, with \(e_0 = 0\), \(|e_k| \leq \frac {Mh^2}{2}\sum _{j<k}(1 + hL)^j = \frac {Mh}{2L}\bigl ((1 + hL)^k - 1\bigr ) \leq \frac {Mh}{2L}(e^{kLh} - 1)\), a discrete Gronwall inequality. The general case is the same with \(\tau _k = O(h^{p+1})\). □
The local errors, each of order \(h^{p+1}\), accumulate over \((T - t_0)/h\) steps into a global error of order \(h^p\), amplified by the factor \(e^{L(T - t_0)}\) of continuous dependence. Practical solvers such as solve_ivp’s default, the Dormand–Prince pair RK45, compute two solutions of orders \(4\) and \(5\) at each step from the same evaluations, use their difference to estimate the local error, and adjust \(h\) to keep it below the tolerance: adaptive step-size control.
Stiffness
Applied to the test equation \(y' = \lambda y\) with \(\Real \lambda < 0\), Euler’s method gives \(y_{k+1} = (1 + h\lambda )y_k\), which decays only if \(|1 + h\lambda | < 1\). For \(\lambda = -1000\) this requires \(h < 0.002\), even when the solution itself varies slowly. A problem is stiff when it contains such fast-decaying components alongside the slow ones of interest: every explicit method has a bounded stability region and must take tiny steps, not for accuracy but for stability. The implicit Euler method \(y_{k+1} = y_k + hf(t_{k+1}, y_{k+1})\) gives \(y_{k+1} = y_k/(1 - h\lambda )\), which decays for every \(h > 0\) when \(\Real \lambda < 0\): it is A-stable. Implicit methods (backward differentiation formulas, implicit Runge–Kutta methods such as Radau) cost more per step, since they solve an algebraic system with Newton’s method, but take far fewer steps on stiff problems (Software 14.7.3). Chemical kinetics, circuit simulation and the method of lines for parabolic partial differential equations all produce stiff systems.
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.