28 Numerical Methods for ODEs
Numerical integration of the initial value problem (IVP): \[ \begin{align} \frac{dy}{dt} = f(t, y), \qquad y(t_0) = y_0. \end{align} \]
28.1 Basic Methods
Advance via forward difference: \[ \begin{align} y_{n+1} = y_n + h f(t_n, y_n). \end{align} \] Error: Local \(O(h^2)\), Global \(O(h)\). Stability: Unstable for stiff problems. Absolute stability requires \(|1 + h\lambda| \leq 1\).
Apply Forward Euler to the test equation \(y' = \lambda y\): \[ \begin{align} y_{n+1} &= y_n + h (\lambda y_n) \\ &= (1 + h\lambda) y_n. \end{align} \] Recursive application gives \(y_n = (1 + h\lambda)^n y_0\). For the solution to remain bounded, we require the amplification factor to satisfy \(|1 + h\lambda| \leq 1\). For real \(\lambda < 0\): \[ \begin{align} -1 &\leq 1 + h\lambda \leq 1 \\ -2 &\leq h\lambda \leq 0 \\ h &\leq \frac{2}{|\lambda|}. \end{align} \] If the step size \(h\) exceeds this limit, the numerical solution oscillates and grows exponentially.
Advance via backward difference: \[ \begin{align} y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}). \end{align} \] Rule: Requires solving a nonlinear system at each step (via Newton). Stability: A-Stable (unconditionally stable for all \(h\) if \(\text{Re}(\lambda) < 0\)). The choice for stiff systems.
28.2 Runge-Kutta Methods
Multistage methods that match higher-order Taylor expansions without requiring analytic derivatives of \(f\).
Evaluates \(f\) at four stages per step: \[ \begin{align} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + h/2, y_n + h k_1 / 2) \\ k_3 &= f(t_n + h/2, y_n + h k_2 / 2) \\ k_4 &= f(t_n + h, y_n + h k_3) \\ y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{align} \] Error: Local \(O(h^5)\), Global \(O(h^4)\). Standard default for non-stiff systems.
Uses an embedded pair to estimate local error \(e\) without extra stages (Dormand-Prince). Mechanism: \(h_{\text{new}} = h_{\text{old}} (\varepsilon / \|e\|)^{1/5}\). The algorithm adjusts the step size \(h\) dynamically to maintain the local error within a specified tolerance.
28.3 Symplectic Integrators
Preserve the geometric structure of Hamiltonian systems (Energy conservation).
For \(H = \frac{1}{2}v^2 + V(q)\), ensures no long-term energy drift:
\(v_{n+1/2} = v_n + \frac{h}{2} F_n\)
\(q_{n+1} = q_n + h v_{n+1/2}\)
\(v_{n+1} = v_{n+1/2} + \frac{h}{2} F_{n+1}\)
Property: Symplectic and time-reversible. Numerical energy oscillates but stays bounded.
Numerical Energy Drift. Standard explicit methods like Forward Euler or RK4 exhibit monotonic energy growth or decay in Hamiltonian systems. Symplectic methods are required for long-term physical simulations to ensure energy stability.
28.4 Exercises
Solve \(y' = -y, y(0)=1\) with Forward Euler for 4 steps. Compare error to \(e^{-1}\).
Verify the stability limit of Forward Euler for \(y' = \lambda y\) is \(h < 2/|\lambda|\).
Show that Forward Euler on a harmonic oscillator leads to unbounded energy growth.
Implement RK45 for \(y' = \cos(t^2)\) and plot the step-size \(h\) over time.