Euler Methods and Basic ODE Solvers

Numerical ODE solvers approximate the solution of the initial value problem $y' = f(t, y)$ , $y(t_0) = y_0$ . The simplest methods -- forward and backward Euler -- illustrate the fundamental trade-offs between accuracy, stability, and computational cost.

One-Step Methods

Definition7.1Forward Euler Method

The forward (explicit) Euler method advances from $y_n \approx y(t_n)$ to $y_{n+1} \approx y(t_{n+1})$ via $y_{n+1} = y_n + h f(t_n, y_n)$ where $h = t_{n+1} - t_n$ is the step size. The method is first-order: the local truncation error (LTE) is $\tau_{n+1} = \frac{h}{2}y''(\xi_n) = O(h)$ , and the global error is $|y(t_n) - y_n| = O(h)$ .

Definition7.2Backward Euler Method

The backward (implicit) Euler method is $y_{n+1} = y_n + h f(t_{n+1}, y_{n+1})$ . This is an implicit equation for $y_{n+1}$ that typically requires Newton's method or fixed-point iteration to solve. Despite being only first-order accurate, backward Euler is A-stable: it remains stable for all step sizes when applied to $y' = \lambda y$ with $\mathrm{Re}(\lambda) < 0$ .

Runge-Kutta Methods

Definition7.3Classical Fourth-Order Runge-Kutta

The RK4 method computes: $k_1 = f(t_n, y_n), \quad k_2 = f(t_n + h/2, y_n + hk_1/2)$ $k_3 = f(t_n + h/2, y_n + hk_2/2), \quad k_4 = f(t_n + h, y_n + hk_3)$ $y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$ RK4 has local truncation error $O(h^4)$ and global error $O(h^4)$ . It requires 4 function evaluations per step, matching the theoretical minimum for an explicit method of order 4.

ExampleRK4 Applied to the Pendulum Equation

For the nonlinear pendulum $\theta'' + \sin\theta = 0$ , written as $y_1' = y_2$ , $y_2' = -\sin(y_1)$ with $\theta(0) = \pi/4$ , $\theta'(0) = 0$ : RK4 with $h = 0.01$ gives 14-digit accuracy over 100 periods, while forward Euler with the same step size diverges after a few periods due to energy drift. With $h = 0.1$ , RK4 still achieves 6-digit accuracy.

Step Size Control

RemarkAdaptive Step Size Selection

Embedded Runge-Kutta methods (e.g., Dormand-Prince, RK45) compute two approximations of different order from the same function evaluations: $y_{n+1}$ (order $p$ ) and $\hat{y}_{n+1}$ (order $p+1$ ). The difference $\|y_{n+1} - \hat{y}_{n+1}\|$ estimates the local error. If this exceeds the tolerance, the step is rejected and $h$ is reduced. The optimal new step size is $h_{\text{new}} = h \cdot \left(\frac{\text{tol}}{\|y_{n+1} - \hat{y}_{n+1}\|}\right)^{1/(p+1)}$ , with safety factors to avoid oscillation.

ExampleStiffness and Implicit Methods

The system $y_1' = -1000 y_1 + y_2$ , $y_2' = y_1 - y_2$ has eigenvalues $\approx -1000$ and $-1$ . The fast mode $e^{-1000t}$ decays rapidly but forces explicit methods to use $h < 2/1000 = 0.002$ for stability, even when only the slow dynamics are of interest. Implicit methods (backward Euler, BDF) allow $h$ determined by accuracy rather than stability, offering speedups of $100\times$ or more for such stiff problems.