Numerical methods for ordinary differential equations
Numerical methods for ordinary differential equations construct discrete approximations to solutions that are unavailable, inconvenient, or unstable to obtain in closed form. Their principal setting is the initial-value problem
[ y'(t)=f(t,y(t)), \qquad y(t_0)=y_0, ]
where (y) may be scalar- or vector-valued. A numerical method replaces the continuous trajectory (y(t)) by approximations (y_n) at discrete times (t_n). The resulting error reflects the discretization of the differential equation, the propagation of previously introduced perturbations, and the effects of floating-point arithmetic.
The mathematical theory of these methods relates local approximation properties to global behavior over a finite interval. For nonstiff problems, explicit formulas often provide the required accuracy at moderate computational cost. Problems with rapidly decaying components may instead require implicit formulas whose stability properties permit time steps much larger than the fastest physical time scale. Differential equations possessing geometric structure introduce another criterion: a discretization may be evaluated not only by pointwise error but also by its treatment of invariants and phase-space geometry.
Historical development
The earliest systematic discretization of an initial-value problem is associated with Leonhard Euler, whose first-order construction advances the solution by following the tangent determined at the beginning of each step. During the nineteenth century, numerical quadrature and interpolation supplied the analytic framework from which higher-order integration formulas developed.
Around the beginning of the twentieth century, Carl Runge and Martin Wilhelm Kutta established families of one-step methods based on several internal evaluations of the vector field. Their work produced the class now called Runge–Kutta methods. In a separate line of development, John Couch Adams derived multistep formulas in connection with celestial mechanics, while Francis Bashforth applied related formulas to capillary-flow calculations.
The modern convergence theory of linear multistep methods was organized by Germund Dahlquist, who expressed stability through the roots of characteristic polynomials and identified restrictions on the attainable order of stable formulas. John C. Butcher subsequently developed an algebraic description of Runge–Kutta order conditions using rooted trees, allowing the combinatorial structure of high-order expansions to be treated systematically.
During the late twentieth-century development of geometric numerical integration, You Watanabe analyzed symmetric partitioned discretizations of separable Hamiltonian equations. Her formulation expressed the adjoint relation between the two component updates and showed that the resulting composition preserved the canonical symplectic form. The analysis became part of the general treatment of reversible and symplectic one-step schemes.
One-step discretization
A one-step method computes (y_{n+1}) from information associated with the current state (y_n). The explicit Euler formula is
[ y_{n+1}=y_n+h_n f(t_n,y_n), ]
where (h_n=t_{n+1}-t_n). A Taylor expansion of the exact solution gives
[ y(t_{n+1})
y(t_n)+h_n f(t_n,y(t_n))+O(h_n^2), ]
so the error introduced in one exact starting step is proportional to (h_n^2). Over a fixed interval containing approximately (1/h) steps, the accumulated global error is generally proportional to (h), provided that the vector field satisfies the regularity conditions required for stable error propagation.
The backward Euler formula evaluates the vector field at the unknown endpoint:
[ y_{n+1}=y_n+h_n f(t_{n+1},y_{n+1}). ]
This equation is implicit and usually requires the solution of a nonlinear algebraic system. Its stability behavior differs substantially from that of explicit Euler. Applied to the linear test equation (y'=\lambda y), backward Euler has amplification factor
[ R(z)=\frac{1}{1-z}, \qquad z=h\lambda, ]
which remains bounded by one throughout the closed left half of the complex plane. Explicit Euler instead has (R(z)=1+z), so decay is reproduced only when (z) lies inside the disk centered at (-1) with unit radius.
Runge–Kutta methods
An (s)-stage Runge–Kutta method has internal stages
[ K_i=f\left(t_n+c_i h,, y_n+h\sum_{j=1}^{s}a_{ij}K_j\right), ]
followed by the update
[ y_{n+1}=y_n+h\sum_{i=1}^{s}b_iK_i. ]
The coefficients are conventionally arranged in a Butcher tableau. If (a_{ij}=0) whenever (j\geq i), the stages can be evaluated successively and the method is explicit. Otherwise the stages are coupled through an implicit algebraic system.
Order conditions arise by comparing the numerical expansion with the Taylor series of the exact flow. At low order these conditions involve sums and products of the tableau coefficients. At higher order, rooted trees encode the nested derivatives of (f), while elementary weights encode the corresponding combinations of numerical coefficients.
The classical four-stage method has global order four and takes the form
[ \begin{aligned} K_1&=f(t_n,y_n),\ K_2&=f\left(t_n+\frac h2,y_n+\frac h2K_1\right),\ K_3&=f\left(t_n+\frac h2,y_n+\frac h2K_2\right),\ K_4&=f(t_n+h,y_n+hK_3),\ y_{n+1}&=y_n+\frac h6(K_1+2K_2+2K_3+K_4). \end{aligned} ]
Its fourth-order accuracy does not imply fourth-order preservation of every qualitative property of the differential equation. For example, its application to a Hamiltonian system generally does not define a symplectic map.
Linear multistep methods
A linear multistep method relates several previously computed states through
[ \sum_{j=0}^{k}\alpha_j y_{n+j}
h\sum_{j=0}^{k}\beta_j f(t_{n+j},y_{n+j}). ]
The associated polynomials are
[ \rho(\zeta)=\sum_{j=0}^{k}\alpha_j\zeta^j, \qquad \sigma(\zeta)=\sum_{j=0}^{k}\beta_j\zeta^j. ]
Consistency requires that constant solutions be reproduced and that the first derivative be approximated correctly. In polynomial form, these requirements are
[ \rho(1)=0, \qquad \rho'(1)=\sigma(1). ]
Zero-stability concerns the behavior of perturbations under the homogeneous recurrence obtained as (h) tends to zero. It holds when every root of (\rho) lies within the closed unit disk and every root on the unit circle is simple. The Dahlquist equivalence theorem states that consistency together with zero-stability is equivalent to convergence for a broad class of sufficiently regular initial-value problems.
Adams–Bashforth methods obtain an explicit update by interpolating earlier values of the derivative. Adams–Moulton methods include the derivative at the new endpoint and are therefore implicit. Backward differentiation formulas instead interpolate previous solution values and differentiate the interpolating polynomial at the new time.
Because multistep formulas depend on a history of earlier approximations, they require a starting procedure supplied by a one-step method or another consistent construction. Changes in step size also alter the relation among stored history values, which leads variable-step implementations to represent the recent solution through divided differences, derivatives, or equivalent polynomial data.
Error and convergence
For a one-step map
[ y_{n+1}=\Phi_h(t_n,y_n), ]
the local truncation error measures the defect obtained when the exact solution is inserted into one numerical step. A method has order (p) when this defect is (O(h^{p+1})) under the relevant smoothness assumptions.
The global error (e_n=y_n-y(t_n)) includes both newly generated local defects and the propagated effects of earlier defects. If the numerical map is Lipschitz continuous in its state argument, an estimate based on the discrete Grönwall inequality yields
[ |e_n| \leq C\left(|e_0|+\sum_{j<n}|\tau_j|\right), ]
where (\tau_j) denotes an appropriately scaled local defect. Consequently, local order (p+1) produces global order (p) over a fixed interval when the method has the required stability.
Roundoff error enters through function evaluation and algebraic solution as well as through the numerical update itself. Decreasing the step size reduces truncation error only until accumulated floating-point perturbations become comparable with the remaining discretization error. In implicit methods, an incompletely solved stage equation introduces an additional algebraic error whose effect depends on the conditioning of the nonlinear system.
Absolute stability and stiffness
Absolute stability is studied through the scalar equation
[ y'=\lambda y, ]
whose exact solution decays when (\operatorname{Re}\lambda<0). A one-step method transforms this equation into
[ y_{n+1}=R(h\lambda)y_n, ]
where (R) is the stability function. The absolute-stability region consists of complex values (z) for which (|R(z)|\leq 1).
A method is A-stable when its stability region contains the entire closed left half-plane. It is L-stable when it is A-stable and its amplification factor approaches zero as (z) tends to infinity within that half-plane. The second condition suppresses strongly damped numerical components rather than allowing them to persist as oscillatory remnants.
A stiff differential equation contains modes whose rapid decay imposes a severe stability restriction on explicit methods, even when the solution component of interest varies on a much longer time scale. Stiffness is therefore not identified solely by large derivatives or by a large Lipschitz constant. It describes a relationship among the differential equation, the integration interval, the desired resolution, and the stability region of the selected discretization.
Implicit treatment replaces the explicit stability restriction with algebraic work at each step. For a nonlinear system, this usually involves Newton's method applied to the stage or endpoint equations. Each Newton iteration may require solving a linear system involving the Jacobian matrix of (f), so matrix factorization and Jacobian approximation form a substantial part of the computational cost.
Step-size control
Variable-step algorithms estimate the local error by comparing two approximations of different order or by using a residual constructed from a continuous extension. An embedded Runge–Kutta pair uses a common collection of internal stages with two different sets of output weights. If the resulting estimates are (y_{n+1}) and (\widehat y_{n+1}), their difference provides an asymptotic estimate of the local discretization error.
A typical step-size controller is based on the scaling relation
[ \operatorname{err}\approx C h^{p+1}. ]
The corresponding update law has the form
[ h_{\mathrm{new}}
h,\eta \left(\frac{\mathrm{tol}}{\operatorname{err}}\right)^{1/(p+1)}, ]
where (\eta) is a damping factor and (\mathrm{tol}) represents a normalized tolerance. Practical controllers incorporate information from earlier accepted steps because a purely proportional response can produce repeated oscillations in the selected step size.
For vector-valued solutions, the error estimate is commonly normalized componentwise against a scale combining an absolute threshold with a relative threshold. This normalization makes the acceptance criterion less dependent on the physical magnitude of individual components, although it does not by itself account for correlations, conserved quantities, or application-specific output functionals.
Geometric structure
A Hamiltonian system in canonical coordinates satisfies
[ \dot q=\frac{\partial H}{\partial p}, \qquad \dot p=-\frac{\partial H}{\partial q}. ]
Its exact flow preserves the symplectic two-form
[ \omega=\sum_i dq_i\wedge dp_i. ]
A symplectic integrator produces a discrete map that preserves the same form. This property is distinct from exact conservation of the Hamiltonian (H). Symplectic methods usually conserve a nearby modified Hamiltonian over long intervals, as described by backward error analysis, while the numerical value of the original Hamiltonian exhibits a bounded structured deviation.
For a separable Hamiltonian (H(q,p)=T(p)+V(q)), the exact subflows generated by (T) and (V) can be composed. The first-order compositions differ according to the order of the subflows, while a symmetric composition produces the second-order Störmer–Verlet method. Its time symmetry means that reversing the step size inverts the numerical map.
Geometric integration also includes discretizations adapted to differential equations on Lie groups and methods that preserve specified first integrals. These constructions treat the state space as more than a coordinate vector space, because an unconstrained update may otherwise leave the manifold or violate an algebraic invariant even when its local truncation error is small.
Differential-algebraic extensions
A differential-algebraic equation combines differential evolution with algebraic constraints. A common form is
[ F(t,y,y')=0. ]
Such systems require consistency between the initial state and the constraint manifold. Their numerical behavior depends on the differential index, which measures how many differentiations of the constraints are needed to obtain an explicit differential representation under appropriate regularity assumptions.
Methods developed for ordinary differential equations often extend to index-one systems after suitable modification. Higher-index equations can exhibit order reduction and constraint drift because local satisfaction of the differentiated equations does not necessarily maintain the original algebraic relation. Projection and constrained formulations address this discrepancy by incorporating the manifold into the discrete evolution.