Mathematical foundations of Euler & Heun Methods.
Leonhard Euler formulated the first numerical algorithm for differential equations in his 1768 treatise Institutiones Calculi Integralis. By approximating the instantaneous derivative dy/dt by a finite difference quotient Δy/h, Euler's method projects the solution along the tangent direction. While conceptually elementary, Euler's method introduces fundamental principles that govern all computational dynamics: the critical distinction between local truncation error (per step) and global error accumulation, the Dahlquist stability barrier, and the trade-offs between explicit and implicit solvers. In this laboratory, we implement explicit Forward Euler, implicit Backward Euler, and the second-order Heun predictor-corrector, map their stability regions, and verify the first-order O(h) global error convergence rate.
Viva voce questions & analytical theory
What is the difference between local truncation error and global truncation error in ODE solvers?
Local truncation error (LTE) is the error committed in a single step assuming the previous solution point was exact; for Forward Euler, LTE is O(h²). Global truncation error (GTE) is the cumulative accumulated error across N = (T - t_0)/h steps. Because errors accumulate over O(1/h) steps, GTE is one order lower: O(h²) · (1/h) = O(h).
How is the stability domain of Forward Euler derived?
Applying Forward Euler to the Dahlquist test equation y' = λy yields y_{n+1} = (1 + λh) y_n. For solutions to remain bounded as n → ∞, the amplification factor must satisfy |R(λh)| = |1 + λh| ≤ 1. In the complex z = λh plane, this corresponds to a disk of radius 1 centered at z = -1.
Why is Backward Euler unconditionally stable (A-stable) for stiff differential equations?
For Backward Euler, y_{n+1} = y_n + λh y_{n+1}, which rearranges to y_{n+1} = y_n / (1 - λh). The amplification factor is R(z) = 1 / (1 - z). Whenever Re(λ) ≤ 0, |1 - z| ≥ 1, so |R(z)| ≤ 1 for the entire left half-plane. This unconditional stability is known as A-stability.
What makes a differential equation 'stiff'?
A system of ODEs y' = f(t, y) is stiff when its Jacobian eigenvalues have negative real parts with vastly different timescales (|Re(λ_max)| ≫ |Re(λ_min)| ≫ 0). Explicit methods are forced to take impractically tiny time steps dictated by stability rather than accuracy, whereas implicit methods remain stable.
How does Heun's method achieve second-order O(h²) global accuracy?
Heun's method uses a predictor step (Forward Euler) to estimate y~_{n+1}, then evaluates an explicit trapezoidal corrector step using the average of slopes at t_n and t_{n+1}. This cancels the O(h) Taylor error term, yielding O(h³) LTE and O(h²) GTE.
What is the Lax Equivalence Theorem in the context of ODE solvers?
Lax's Equivalence Theorem (and Dahlquist's Equivalence Theorem for multistep methods) states that for a consistent numerical method (LTE → 0 as h → 0), zero-stability is both necessary and sufficient for global convergence to the true solution.
Why is Forward Euler rarely used in production scientific computing?
Forward Euler has very low first-order accuracy O(h) and a tiny conditional stability region. Achieving high precision requires millions of tiny time steps, causing excessive CPU time and floating-point roundoff accumulation compared to higher-order Runge-Kutta or multistep methods.
Algorithm precautions and stability rules
- Always check the stability restriction h < 2 / |λ_max| before running explicit Forward Euler on dissipative systems to prevent catastrophic divergence.
- For stiff differential equations with large negative eigenvalues, use implicit Backward Euler to avoid step size collapse.
- Ensure floating point precision remains 64-bit IEEE 754 during iterative time marching to minimize accumulated roundoff error.
- When comparing errors across step sizes, ensure the final evaluation time T is matched exactly (t = N · h).
- In Backward Euler, solve the implicit nonlinear equation using Newton-Raphson iteration with a convergence tolerance of 10⁻¹².
- Do not extrapolate ODE solutions beyond the domain where the right-hand vector field f(t, y) satisfies the Lipschitz continuity condition.
Theoretical limitations
- Forward Euler requires very small step sizes to maintain accuracy and stability, making it computationally inefficient for precision applications.
- Heun's method is only conditionally stable, still requiring small step sizes on stiff problems.
- Euler methods are not symplectic: when applied to Hamiltonian orbital mechanics, they cause artificial energy drift over long integration times.