ODE Solvers

ODE solvers approximate ordinary differential equations by advancing data through small steps. In quantum mechanics they appear in radial equations, separated one-dimensional boundary-value problems, semiclassical trajectories, classical equations used in approximations, and finite-dimensional time-dependent systems.

This page owns the general numerical ODE machinery: Runge–Kutta methods, adaptive step-size control, stiffness warnings, and shooting methods. The special structure of real-time Schrödinger propagation is treated separately in Time-Stepping Methods and Matrix Exponentials Numerically.

Initial-Value Problems

The standard initial-value problem is

\frac{dy}{dt} = f(t,y), \qquad y(t_0)=y_0.

Here $y$ may be a scalar, a real vector, or a complex vector. A numerical solver produces approximations

y_n \approx y(t_n), \qquad t_n=t_0+nh.

A one-step method has the form

y_{n+1} = \Phi_h(t_n,y_n),

where $h$ is the step size and $\Phi_h$ is the numerical update rule.

Higher-order ODEs are usually converted to first-order systems before using a general solver. For example, the one-dimensional stationary Schrödinger equation

- \frac{\hbar^2}{2m} \frac{d^2\psi}{dx^2} + V(x)\psi = E\psi

can be written with

y_1=\psi, \qquad y_2=\psi',

so that

\frac{d}{dx} \begin{pmatrix} y_1\\ y_2 \end{pmatrix} = \begin{pmatrix} y_2\\ \dfrac{2m}{\hbar^2}\bigl(V(x)-E\bigr)y_1 \end{pmatrix}.

This conversion is useful for shooting methods, but it is not the only numerical approach. Finite-difference and spectral discretizations solve the same eigenvalue problem by constructing a matrix; see Finite Difference Methods and Spectral Methods.

Local and Global Error

An ODE method has local truncation error of order $p+1$ if one exact input step satisfies

y(t_{n+1}) - \Phi_h(t_n,y(t_n)) = O(h^{p+1}).

Under suitable smoothness and stability assumptions, the accumulated global error over a fixed interval is usually order $p$ :

y_n-y(t_n) = O(h^p).

This distinction matters in convergence tests. A fourth-order method generally has local error $O(h^5)$ and global error $O(h^4)$ , not $O(h^5)$ globally.

Error constants can still be large. A high-order method is not automatically accurate if the solution has sharp features, singular points, rapid phase oscillations, or unstable directions.

Explicit Euler as a Warning Example

The explicit Euler method is

y_{n+1} = y_n+h f(t_n,y_n).

It is first order and cheap, but it is rarely a good production solver. For the test equation

y'=\lambda y,

Euler gives

y_{n+1} = (1+h\lambda)y_n.

The amplification factor $1+h\lambda$ can grow even when the exact solution is bounded or decaying. This is why stability must be analyzed separately from formal order.

In real-time Schrödinger dynamics, explicit Euler is especially poor because it creates artificial norm growth for nonzero energy modes; see Time-Stepping Methods.

Runge–Kutta Methods

Runge–Kutta methods evaluate the right-hand side several times inside a step. The classical fourth-order method is

\begin{aligned} k_1&=f(t_n,y_n),\\ k_2&=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1\right),\\ k_3&=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}k_2\right),\\ k_4&=f(t_n+h,y_n+hk_3), \end{aligned}

followed by

y_{n+1} = y_n + \frac{h}{6} \left( k_1+2k_2+2k_3+k_4 \right).

This method is often a good first test for smooth, nonstiff initial-value problems. It is easy to implement and has global fourth-order accuracy when its assumptions are met.

For a linear autonomous problem $y'=Ay$ , a Runge–Kutta method can be described by a stability function $R$ :

y_{n+1} = R(hA)y_n.

For the scalar equation $y'=\lambda y$ , this becomes

y_{n+1} = R(h\lambda)y_n.

The method is stable on a given mode only if $R(h\lambda)$ has acceptable size and phase behavior for that mode. In quantum work this matters both for dissipative auxiliary equations and for oscillatory equations with nearly imaginary eigenvalues.

Adaptive Step-Size Control

Adaptive solvers estimate their own local error and change the step size. Many practical solvers use an embedded Runge–Kutta pair: two formulas share the same intermediate evaluations but have different orders.

Let $y_{n+1}^{(p)}$ and $y_{n+1}^{(p+1)}$ be two estimates. A basic error estimate is

e_n = y_{n+1}^{(p+1)} - y_{n+1}^{(p)}.

For a vector problem, compare components using absolute and relative tolerances:

s_i = \mathrm{atol}_i + \mathrm{rtol}_i \max\left( \lvert y_{n,i}\rvert, \lvert y_{n+1,i}\rvert \right).

A typical normalized error measure is

\epsilon_n = \left( \frac{1}{d} \sum_{i=1}^{d} \left\lvert \frac{e_{n,i}}{s_i} \right\rvert^2 \right)^{1/2}.

If $\epsilon_n\le 1$ , the step is accepted. If $\epsilon_n\gt1$ , the step is rejected and retried with a smaller $h$ . A common step-size update has the form

h_{\mathrm{new}} = \eta h \epsilon_n^{-1/(p+1)},

with a safety factor $\eta\lt1$ .

Adaptive methods are convenient, but the tolerance is not a physical accuracy guarantee. If the wrong variables are scaled, if a boundary condition is not checked, or if a conserved quantity drifts, a small local error estimate can still produce an untrustworthy quantum result.

Stability and Stiffness

Stability asks whether numerical errors are amplified by the algorithm. The scalar test equation

y'=\lambda y

is a useful first diagnostic. If the exact solution decays because $\operatorname{Re}\lambda\lt0$ , a method should not turn that decay into numerical growth.

An ODE is stiff when stable and unstable time scales are widely separated, forcing an explicit method to take very small steps for stability rather than accuracy. Stiffness can appear in:

imaginary-time propagation, where high-energy modes decay rapidly;
master-equation or dissipative models with fast relaxation channels;
radial equations with singular coefficients;
multiscale semiclassical equations;
boundary-layer problems.

Implicit methods, backward differentiation formulas, and implicit Runge–Kutta methods are often used for stiff systems. They require nonlinear or linear solves inside each step, so their accuracy depends on solver tolerances as well as the nominal ODE step size.

For closed real-time Schrödinger evolution, stiffness is not the only concern. A method can be stable yet nonunitary, or norm-preserving yet phase-inaccurate. That structure is handled in Time-Stepping Methods.

Shooting Methods

A shooting method turns a boundary-value problem into a root-finding problem. Suppose an ODE depends on an unknown parameter $E$ and must satisfy boundary conditions at $a$ and $b$ . Choose initial data at $a$ , integrate to $b$ , and measure the boundary mismatch.

For a one-dimensional bound-state problem, one may define a residual such as

F(E) = \psi_E(b),

when the desired right boundary condition is $\psi(b)=0$ . Eigenvalues are then roots:

F(E)=0.

In practice, shooting is more delicate than this formula suggests:

exponentially growing forbidden-region components can dominate the numerical solution;
normalization is usually imposed after the eigenvalue is found;
even and odd parity states require different initial data;
node counts help identify which eigenvalue branch is being followed;
singular endpoints need asymptotic initial conditions rather than arbitrary finite values.

For example, an even bound state in a symmetric potential can be shot from the origin with

\psi(0)=1, \qquad \psi'(0)=0,

while an odd state uses

\psi(0)=0, \qquad \psi'(0)=1.

The arbitrary scale is harmless because the Schrödinger equation is linear. After the correct energy has been found, normalize the wavefunction using Numerical Quadrature.

Numerov for Second-Order Linear Equations

For equations of the form

y''(x)=q(x)y(x),

the Numerov method gives a high-order three-point recurrence:

\left( 1-\frac{h^2}{12}q_{n+1} \right)y_{n+1} = 2\left( 1+\frac{5h^2}{12}q_n \right)y_n - \left( 1-\frac{h^2}{12}q_{n-1} \right)y_{n-1}.

For the time-independent Schrödinger equation,

q(x) = \frac{2m}{\hbar^2} \bigl(V(x)-E\bigr).

Numerov is useful for smooth one-dimensional or radial bound-state equations, especially inside shooting workflows. It is not a universal ODE solver: it assumes a special second-order form, fixed grid spacing in its basic version, and enough smoothness for the high-order cancellation to work.

Boundary-Value Problems Are Not Just Initial-Value Problems

Quantum eigenvalue problems usually impose boundary conditions at more than one point. An initial-value solver by itself cannot decide the allowed energy; it only propagates data for a chosen energy.

The mathematical problem is closer to

\mathcal L\psi = Ew\psi, \qquad B_a\psi=0, \qquad B_b\psi=0,

where $B_a$ and $B_b$ encode boundary data. Shooting attacks this by repeated IVP solves. Finite-difference, spectral, finite-element, and Galerkin approaches attack it by forming a finite eigenvalue problem directly. The operator viewpoint is developed in Eigenvalue Problems and Sturm–Liouville Theory.

Conservation and Structure Checks

ODE solvers should be tested against structures the differential equation is supposed to preserve.

For a closed quantum state, check norm conservation:

\langle\psi_n,\psi_n\rangle \approx \langle\psi_0,\psi_0\rangle.

For a time-independent closed system, check energy drift:

\langle\psi_n,H\psi_n\rangle \approx \langle\psi_0,H\psi_0\rangle.

For classical Hamiltonian trajectories, check conserved energy and, when relevant, symplectic structure. Generic high-order ODE solvers can be accurate over short times while accumulating long-time qualitative drift. Structure-preserving methods are often worth the extra care when the qualitative invariant is the point of the calculation.

Validation Workflow

A reliable ODE calculation should include several independent checks:

compare against an exactly solvable problem, such as the harmonic oscillator or a two-level system;
refine the step size and verify the expected error trend;
vary absolute and relative tolerances separately;
check boundary residuals for shooting problems;
verify normalization with the correct quadrature weights;
compare shooting eigenvalues against a matrix discretization when possible;
monitor conserved quantities or monotone quantities;
test sensitivity to the finite domain and boundary placement.

One clean benchmark is the harmonic oscillator. A shooting calculation should reproduce the low-lying energies

E_n = \hbar\omega \left( n+\frac{1}{2} \right),

and should produce wavefunctions with the correct parity and node count.

Choosing a Solver

Situation	Common first choice	Main check
smooth nonstiff IVP	adaptive Runge–Kutta	step refinement and conserved quantities
highly oscillatory phase	method tailored to oscillations or exponentials	phase error and invariant drift
stiff dissipative equation	implicit or stiff-aware method	solver tolerance and stability
one-dimensional bound state	shooting with bracketing or Numerov	boundary residual and node count
large discretized quantum dynamics	structure-aware time step or exponential action	unitarity, sparsity, and spectral range
singular radial equation	asymptotic start plus controlled integration	endpoint behavior and normalization

No method is best in isolation. The right solver is chosen for the equation, boundary data, stiffness, required observable, and validation strategy.

Common Mistakes

Treating an eigenvalue boundary problem as if one arbitrary initial-value solve determines the energy.
Judging a solver only by its formal order while ignoring stability.
Using adaptive tolerances without checking physical invariants.
Forgetting that a wavefunction’s arbitrary shooting scale must be normalized afterward.
Using explicit methods on stiff dissipative problems and interpreting tiny step sizes as physical necessity.
Shooting across a classically forbidden region without controlling exponential growth.
Ignoring singular endpoint behavior in radial equations.
Comparing wavefunctions before fixing phase, normalization, and grid conventions.
Assuming local error estimates measure the error in the final observable.

Cross-Links

References

E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer, 1993.
E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed., Springer, 1996.
U. M. Ascher and L. R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, SIAM, 1998.
J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, 3rd ed., Springer, 2002.
J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.
T. Pang, An Introduction to Computational Physics, 2nd ed., Cambridge University Press, 2006.

Exercises

Convert the stationary Schrödinger equation to a first-order system.

For

- \frac{\hbar^2}{2m}\psi''(x) + V(x)\psi(x) = E\psi(x),

define $y_1=\psi$ and $y_2=\psi'$ . Write the first-order system.

Solution

The first equation is

y_1'=y_2.

Solving the Schrödinger equation for $\psi''$ gives

\psi'' = \frac{2m}{\hbar^2} \bigl(V(x)-E\bigr)\psi.

Thus

y_2' = \frac{2m}{\hbar^2} \bigl(V(x)-E\bigr)y_1.

The system is

\frac{d}{dx} \begin{pmatrix} y_1\\ y_2 \end{pmatrix} = \begin{pmatrix} y_2\\ \dfrac{2m}{\hbar^2}\bigl(V(x)-E\bigr)y_1 \end{pmatrix}.

Show that explicit Euler has first-order global accuracy for $y'=\lambda y$ in the small-step expansion.

Solution

One Euler step gives

y_{n+1} = (1+h\lambda)y_n.

The exact one-step factor is

e^{h\lambda} = 1+h\lambda+\frac{h^2\lambda^2}{2}+O(h^3).

The one-step difference is $O(h^2)$ , so the local truncation error is second order. Over $O(1/h)$ steps on a fixed interval, this usually accumulates to first-order global error, assuming the problem is stable.

Why do adaptive solvers combine absolute and relative tolerances?

Solution

A purely relative tolerance becomes meaningless when a component passes near zero, because the allowed error would also approach zero. A purely absolute tolerance can be too loose for large components. The scale

s_i = \mathrm{atol}_i + \mathrm{rtol}_i \max\left( \lvert y_{n,i}\rvert, \lvert y_{n+1,i}\rvert \right)

allows small components to be controlled by an absolute floor and large components to be controlled proportionally.

In a symmetric potential, what shooting initial data would you use for even and odd states at the origin?

Solution

For an even state, the derivative vanishes at the origin, so a convenient arbitrary scale is

\psi(0)=1, \qquad \psi'(0)=0.

For an odd state, the value vanishes at the origin, so one can use

\psi(0)=0, \qquad \psi'(0)=1.

The constants $1$ are not physical normalizations. They only set a nonzero scale for the shooting integration.

Why can a locally accurate adaptive ODE solution still give a wrong bound-state energy?

Solution

The adaptive solver controls the local error for a chosen initial value and chosen energy. A bound-state energy is determined by boundary conditions at both ends. If the shooting residual, domain size, singular endpoint behavior, node count, or normalization is wrong, a locally accurate IVP integration can still correspond to the wrong boundary-value problem.