ODE Integrators II: Implicit and Symplectic Methods

Convergence of Integrators¶

When solving ordinary differential equations (ODEs) numerically, errors inevitably enter through truncation (e.g., approximating derivatives with finite differences) and round-off (finite precision arithmetic). If these errors grow uncontrollably, the computed solution may diverge from the true behavior, even if the method is very accurate for short times.

It is useful to distinguish between several related ideas:

• Accuracy: How close the numerical solution is to the exact solution at a given time step.

• Consistency: The numerical method reproduces the original ODE as $\Delta t \to 0$. Formally, the local truncation error (the error made in a single step assuming exact input) should vanish as $\Delta t \to 0$.

• Stability: Errors introduced during computation do not grow uncontrollably as steps are repeated.

• Convergence: The global numerical solution approaches the exact solution as $\Delta t \to 0$.

These concepts are related by a central result of numerical analysis:

This is known as the Lax Equivalence Theorem. The theorem is really for linear finite difference methods for partial different equation, but it is still useful to discuss it in ODE integrator. It gives us a practical recipe:

• First check that the method is consistent (usually straightforward).

• Then analyze stability, we usually use the linear test equation:

  $$\begin{align} \frac{dx}{dt} = \lambda x, \quad \lambda \in \mathbb{C}. \end{align}$$

  (2)

  Its exact solution is $x(t) = x_0 e^{\lambda t}$. Applying a numerical method produces an update:

  $$\begin{align} x_{n+1} = R(z) x_n, \quad z = \lambda \Delta t, \end{align}$$

  (3)

  where $R(z)$ is the amplification factor. A method is stable if

  $$\begin{align} |R(z)| \leq 1, \end{align}$$

  (4)

  so that errors do not amplify step by step. The set of $z$ satisfying this condition defines the stability region.

• Convergence then follows automatically.

Explicit Methods¶

Forward Euler Stability and Consistency¶

The Forward Euler method for $dx/dt = f(x)$ is

Expanding the true solution with a Taylor series:

The difference between the true solution and the Forward Euler step is $\mathcal{O}(\Delta t^2)$. Thus, the local truncation error is $\mathcal{O}(\Delta t^2)$, meaning the method is first-order consistent.

Applying the test equation to forward Euler method gives:

so the amplification factor is

The stability condition requires

which is the interior of the circle centered at $(-1,0)$ with radius 1 in the complex plane.

Combined with stability, forward Euler is first-order convergent, i.e., the global error scales as $\mathcal{O}(\Delta t)$, within its stable region. Outside the stable region, the method diverges.

Exponential Decay¶

To illustrate stability and convergence, let’s recall the forward Euler scheme we implemented last time, except we allow explicit time dependency in the right hand side of the function.

We apply it to solve the equation $dx/dt = -x$, which has $\lambda = -1$ and solution $\exp(-t)$.

The error for $\Delta t = 4$ grows unbounded. Even at the critical stability step size $\Delta t = 2$, the oscillating behavior is fundamentally wrong. Nevertheless, as $\Delta t \rightarrow 0$, especially when $\Delta t < -2/\lambda$ the solution does convergence.

However, if $\lambda$ is positive, there is no positive $\Delta t$ that can make forward Euler stable. And we call forward Euler unconditionally unstable.

This suggests that,

• Consistency: Forward Euler has local error $\mathcal{O}(\Delta t^2)$.

• Stability: Stable only inside the unit disk centered at $(-1,0)$.

• Convergence: By Lax’s theorem, Forward Euler converges with global error $\mathcal{O}(\Delta t)$.

Implications:

• For $\lambda < 0$ (decaying solutions), stability requires $\Delta t \leq -2/\lambda$.

• For $\lambda > 0$ (growing solutions), Forward Euler is unconditionally unstable: no choice of $\Delta t$ can control error growth.

• For stiff systems, where $|\lambda|$ is very large, the step size restriction is too severe to be practical.

This motivates the need for implicit integrators (Backward Euler, trapezoidal rule) and later symplectic integrators for Hamiltonian systems.

Simple Harmonic Oscillator¶

The exponential test equation matches the linear test equation and is a good illustrative example, but many physical systems of interest are naturally written as systems of coupled ODEs.

Recalling again from last lecture, an important example is the simple harmonic oscillator:

Here $\theta(t)$ is the angular displacement, $\Omega(t)$ the angular velocity, and $g/l$ sets the natural frequency $\omega^2 = g/l$.

Applying the Forward Euler method gives:

This can be written compactly as:

Here $A$ is the amplification matrix, the multi-dimensional analogue of the scalar amplification factor $R(z)$. Stability requires that the eigenvalues of $A$ satisfy $|\lambda| \leq 1$.

Let’s compute them. The characteristic polynomial of $A$ is:

Setting this to zero, $(1 - \lambda)^2 = -(g/l)\Delta t^2$. Thus the eigenvalues are:

The magnitude of each eigenvalue is:

This is always larger than 1 for any nonzero time step. Therefore, Forward Euler applied to the harmonic oscillator is unconditionally unstable.

The true solution of the harmonic oscillator is bounded and oscillatory. However, the forward Euler discretization produces eigenvalues with modulus larger than 1. This means that the numerical solution grows exponentially in amplitude while also oscillating. This “blow up” happens even for very small $\Delta t$.

For the harmonic oscillator, Forward Euler is inconsistent with the bounded, energy-conserving nature of the physics. In fact,

• Explicit (forward) Euler introduces artificial growth of energy.

• Implicit (backward) Euler introduces artificial damping of energy.

• Symplectic integrators (e.g., Leapfrog, see below) preserve the energy structure and produce qualitatively correct long-term behavior.

This example connects stability analysis directly to physics: the choice of integrator is not just about error size, but about preserving the qualitative features of the dynamics.

Stability Regions for RK2 and RK4¶

So far we have analyzed the forward Euler method, which is simple but has a very restrictive stability region. In practice, physicists often use higher-order explicit Runge-Kutta (RK) methods, such as RK2 and RK4. These methods improve accuracy while also enlarging the stability region, allowing for larger time steps in many problems.

To analyze stability, we again apply the methods to the linear test equation:

Each Runge-Kutta scheme can be written in the form:

where $R(z)$ is the stability function of the method. Stability requires $|R(z)| \leq 1$.

• Forward Euler (for comparison):

  $$\begin{align} R_{\text{fE}}(z) = 1 + z \end{align}$$

  (18)

• RK2 (Heun’s method):

  $$\begin{align} R_{\text{RK2}}(z) = 1 + z + \tfrac{1}{2} z^2 \end{align}$$

  (19)

• RK4 (classical Runge–Kutta):

  $$\begin{align} R_{\text{RK4}}(z) = 1 + z + \tfrac{1}{2} z^2 + \tfrac{1}{6} z^3 + \tfrac{1}{24} z^4 \end{align}$$

  (20)

Notice: $R_{\text{RK4}}(z)$ is simply the first four terms of the Taylor expansion of $e^z$, which explains why RK4 is accurate to fourth order.

Some observations from the above plot:

• Forward Euler:

  • Stability region: disk of radius 1 centered at $(-1,0)$.

  • Very restrictive: only stable for very small $\Delta t$ when $\lambda$ has a negative real part.

• RK2

  • Stability region: significantly larger, extending farther left into the negative real axis.

  • Stable for $\text{Re}(z) \gtrsim -2$.

  • This allows larger time steps for decaying problems, while still being conditionally stable.

• RK4

  • Stability region: even larger, extending roughly to $\text{Re}(z) \approx -2.8$ on the real axis.

  • Includes a more extensive area in the complex plane, meaning oscillatory systems can tolerate larger $\Delta t$ before instability.

  • Most widely used in physics and engineering because of its balance between accuracy and stability.

In general,

• Order vs. Stability: Higher-order RK methods expand the stability region, allowing larger $\Delta t$ while maintaining stability.

• Still Conditionally Stable: Neither RK2 nor RK4 are A-stable. They cannot remain stable for arbitrarily stiff problems (where $\lambda$ is very negative).

• Physics Implication:

  • For non-stiff problems (e.g. orbital dynamics, wave propagation), RK4 is an excellent choice.

  • For stiff problems (e.g. radiative cooling, chemical networks), even RK4 will fail unless $\Delta t$ is extremely small. These problems require implicit methods.

Implicit Methods¶

So far we have focused on explicit methods such as forward Euler, RK2, and RK4. They are simple, intuitive, and accurate for many problems. However, we also discovered:

• Forward Euler fails on oscillatory systems like the harmonic oscillator (unconditionally unstable).

• Explicit RK methods extend the stability region, but they are still conditionally stable. They cannot handle arbitrarily stiff problems.

This raises a natural question:

Can we design methods that remain stable even when $\lambda \Delta t$ is very large and negative?

This leads us to implicit methods.

In explicit methods, the new state is computed directly from known information,

Example includes forward Euler $x_{n+1} = x_n + \Delta t f(x_n, t_n)$. This makes them cheap and easy per step, but with limited stability.

In implicit methods, the new state appears on both sides of the equation

Hence, root finders are required at each step. Example includes backward Euler $x_{n+1} = x_n + \Delta t f(x_{n+1}, t_{n+1})$. This makes them expensive per step but has much larger stability regions. They are ideal for stiff systems.

Backward Euler Method¶

The simplest implicit scheme is backward Euler, the implicit counterpart to Forward Euler.

This requires solving an equation for $x_{n+1}$ at every step.

For its stability analysis with the linear test equation, we apply backward Euler to $dx/dt = \lambda x$ and define $z = \lambda \Delta t$. Given the update rule,

Rearrange:

the amplification factor is:

The stability condition requires

• If $\text{Re}(z) \leq 0$ (stable in the continuous system), this inequality always holds.

• Therefore, backward Euler is A-stable: stable for the entire left half-plane.

This is in stark contrast with forward Euler and explicit RK methods.

Implementing the backward Euler method involves solving an implicit equation at each time step. The approach depends on whether the ODE is linear or nonlinear.

• Linear ODEs: If $f(x, t)$ is linear in $x$, the implicit equation can be solved directly without iterative methods. For example, for $dx/dt = \lambda x$:

  $$\begin{align} x_{n+1} = (1 - \lambda \Delta t)^{-1} x_n. \end{align}$$

  (28)

  Using matrix inverse, the above formulation even works for system of linear ODEs.

• Nonlinear ODEs: If $f(x, t)$ is nonlinear, solving for $x_{n+1}$ typically requires iterative methods such as the Newton-Raphson method. Python’s fsolve() from the scipy.optimize library can also be used for this purpose.

Stiff ODEs¶

Let’s implement the backward Euler method for a simple linear ODE and compare it with the forward Euler method.

Consider:

The exact solution is:

This ODE is stiff because it contains terms with vastly different decay rates ($e^{-1000t}$ vs. $e^{-t}$).

The comparison between forward Euler and backward Euler highlights their complementary strengths and weaknesses:

• Forward Euler is simple and inexpensive per step, but it fails badly on stiff problems. Even with small time steps, it can oscillate and diverge when fast-decaying modes are present.

• Backward Euler, in contrast, remains stable regardless of the step size for problems with $\text{Re}(\lambda) \leq 0$. This A-stability makes it ideally suited for stiff systems, where explicit schemes would require prohibitively small $\Delta t$. The price we pay is higher computational cost per step (solving implicit equations) as mentioned above, and potentially reduced accuracy.

In practice, implicit schemes are the method of choice whenever stability is the dominant concern. For example, in problems with a wide separation of timescales such as chemical kinetics, radiative cooling, or diffusion processes. They allow us to take much larger time steps while still capturing the slow dynamics faithfully.

However, backward Euler introduces artificial damping in oscillatory systems, suppressing dynamics that should remain bounded. This motivates the search for implicit methods that retain the stability benefits of backward Euler but reduce the damping error. One elegant solution is the implicit midpoint method (also known as Crank–Nicolson), which we will study next.

Implicit Midpoint Method¶

The implicit midpoint method is another implicit integrator. It offers second-order accuracy and is symplectic (see below). This makes it especially valuable for problems with an underlying Hamiltonian structure, such as oscillators and orbital dynamics, where energy conservation is critical over long simulations. It can be viewed as the midpoint analogue of backward Euler: instead of evaluating the derivative at the end of the interval, it evaluates it at the midpoint between $x_n$ and $x_{n+1}$.

For the ODE

the implicit midpoint update is

The derivative is evaluated at the midpoint of the interval $[t_n, t_{n+1}]$. This improves accuracy to second order, compared with the first-order backward Euler scheme.

Apply to the linear test equation $dx/dt = \lambda x$:

Rearranging:

Thus,

with amplification factor

• If $\text{Re}(z) \leq 0$, then $|R(z)| \leq 1$.

• Therefore, the Implicit Midpoint method is A-stable, like backward Euler.

• But unlike backward Euler, it does not introduce artificial damping in oscillatory systems. Instead, it preserves oscillatory behavior.

This makes it both stable for stiff problems and faithful to Hamiltonian dynamics over long time integrations.

Let’s apply the implicit midpoint method to the same stiff ODE as before:

From the above plots,

• Forward Euler: unstable unless $\Delta t$ is tiny.

• Backward Euler: stable but introduces extra damping.

• Implicit midpoint: stable and preserves oscillatory features.

This balance of stability and geometric fidelity makes the implicit midpoint method (also known as Crank-Nicolson) a popular scheme for both stiff problems and Hamiltonian systems.

However, the implicit nature makes it less ideal for non-stiff Hamiltonian systems. There there methods that take advantage of Hamiltonian?

Symplectic Algorithms¶

When simulating physical systems, especially those derived from Hamiltonian mechanics, it is not enough for a numerical method to be accurate and stable in the short term. We also want to preserve the geometric structure of the system to ensure that invariants such as energy and momentum remain bounded and that phase-space volume is preserved.

Symplectic integrators are designed precisely for this purpose. They preserve the symplectic structure of Hamiltonian flows, which guarantees long-term qualitative fidelity, i.e., bounded oscillations, stable orbits, and correct statistical properties.

Hamiltonian Systems¶

For a Hamiltonian

with

• $q$: generalized coordinates (e.g. position).

• $p$: conjugate momenta.

• $T(p)$: kinetic energy.

• $V(q)$: potential energy.

The system’s dynamics are governed by Hamilton’s equations:

Examples of Hamiltonian system include:

• Simple harmonic oscillator:

  $$\begin{align} H(q,p) = \frac{p^2}{2m} + \frac{1}{2}k q^2 \end{align}$$

  (39)

• Kepler problem (planetary motion):

  $$\begin{align} H(q,p) = \frac{p^2}{2m} - \frac{GMm}{|q|} \end{align}$$

  (40)

Conservation Laws¶

Hamiltonian systems obey key invariants:

• Energy conservation: if $H$ is not explicitly time-dependent, $dH/dt = 0$.

• Momentum conservation: if $H$ does not depend on a coordinate $q_i$, the corresponding $p_i$ is conserved (Noether’s theorem).

• Phase-space volume preservation: Liouville’s theorem states that Hamiltonian flows preserve volume in $(q,p)$ space.

Why Symplectic Integrators?¶

• Standard methods (e.g. Forward Euler, RK4) often cause systematic drift in energy over long runs. In some sense, unstable schemes artificially inject energy to the system, while stable schemes (including RK4) artificially remove energy from the system.

• Symplectic integrators keep the energy error bounded and oscillatory, never drifting away permanently.

• This makes them especially valuable in:

  • Molecular dynamics (bounded energy in large ensembles).

  • Celestial mechanics (stable planetary orbits).

  • Long-term oscillator simulations.

Symplectic Euler Method¶

The symplectic Euler method updates the coordinates and momenta in a staggered fashion. There are two variants of the Symplectic Euler method:

1. Implicit Symplectic Euler:

   $$\begin{align} p_{n+1} &= p_n - \Delta t \cdot \frac{\partial V}{\partial q}(q_{n}) \\ q_{n+1} &= q_n + \Delta t \cdot \frac{p_{n+1}}{m} \end{align}$$

   (41)

2. Explicit Symplectic Euler:

   $$\begin{align} q_{n+1} &= q_n + \Delta t \cdot \frac{p_n}{m} \\ p_{n+1} &= p_n - \Delta t \cdot \frac{\partial V}{\partial q}(q_{n+1}) \end{align}$$

   (42)

As the formulation suggests, they are both first-order accurate in time. Nevertheless, compared to forward or even backward Euler methods, they preserve the symplectic structure and ensure bounded energy errors. They also offers better stability compared to non-symplectic explicit methods.

From the above plot, we observe:

• Forward Euler Method:

  • Behavior: Deviates significantly from the exact solution over time, exhibiting growing oscillations.

  • Reason: Non-symplectic nature causes energy to accumulate or dissipate artificially.

• Symplectic Euler Method:

  • Behavior: Maintains bounded energy errors, resulting in stable, oscillatory trajectories that closely follow the exact solution.

  • Reason: Preservation of the symplectic structure ensures that energy errors do not systematically drift.

Leapfrog (Verlet) Integrator¶

The Leapfrog (or Verlet) method is one of the most widely used symplectic integrators. It improves upon symplectic Euler by achieving second-order accuracy while still being explicit, symplectic, and time-reversible.

The key idea: positions $q$ and momenta $p$ are evaluated at staggered time points, so that they “leapfrog” over each other. This staggering ensures stability and bounded energy error, even over very long simulations.

Leapfrog is especially important in molecular dynamics, celestial mechanics, and $n$-body simulations, where millions of steps are required to track conservative systems.

Start with Hamilton’s equations:

Momenta are updated on half steps:

then positions are updated on integer steps:

Thus $p$ “leapfrogs” over $q$, each using the other’s most up-to-date value.

Synchronized (velocity-Verlet) Form¶

To avoid half-step bookkeeping, we can write Leapfrog in a synchronized form:

This form is known as velocity-Verlet, and is the most widely used in practice.

Kick-Drift-Kick Form¶

Another equivalent rearrangement is the kick-drift-kick scheme:

This is especially convenient if time steps are adaptive, because only one new force evaluation is needed per step.

The Leapfrog integrator is:

• Accuracy: Second-order in $\Delta t$ (vs. first-order for Symplectic Euler).

• Symplectic: Preserves the Hamiltonian structure, ensuring bounded energy error.

• Time-reversible: Running with $+\Delta t$ then $-\Delta t$ recovers the initial state exactly.

• Stability: Stable for oscillatory motion provided $\Delta t < 2/\omega$ (where $\omega$ is the oscillator frequency).

Derivation and Preservation of Symplectic Structure¶

So far, we have observed that symplectic integrators behave differently: they keep energy bounded, preserve phase-space volume, and remain stable over long times. But why does this happen?

The answer lies in their ability to preserve the symplectic structure of Hamiltonian systems. This structure encodes both the conservation laws (via Noether’s theorem) and Liouville’s theorem (phase-space volume preservation). Understanding how an algorithm preserves this structure explains why symplectic methods are uniquely powerful.

For a separable Hamiltonian

Hamilton’s equations are

A one-step method can be written as a map

For symplectic Euler:

The Jacobian of this transformation is

A method is symplectic if this Jacobian preserves the canonical 2-form:

Carrying out the matrix multiplication confirms that the symplectic condition holds. Thus symplectic Euler is indeed symplectic.

The geometric interpretations here are:

• Symplecticity means the numerical flow preserves the geometry of Hamiltonian mechanics. Phase-space volume is preserved exactly, and invariants such as energy remain bounded.

• Contrast with non-symplectic methods: Forward Euler and RK4 distort the phase-space geometry, leading to secular energy drift.

• Discrete Hamiltonian Flow: Symplectic integrators can be interpreted as approximating a nearby Hamiltonian system rather than the exact one. This explains why energy error oscillates but does not drift.

Some practical considerations for using symplectic intergrators include:

• Choosing $\Delta t$: Symplectic integrators allow larger steps than non-symplectic ones, but accuracy still requires $\Delta t$ small enough to resolve the dynamics.

• Nonlinear forces: For nonlinear $V(q)$, implicit symplectic schemes (e.g. implicit midpoint) may require solving nonlinear equations each step.

• Adaptive time-stepping: Standard symplectic schemes assume constant $\Delta t$. Adaptive step sizes break symplecticity, though specialized adaptive symplectic schemes exist.

Summary¶

• Symplectic integrators preserve the symplectic two-form, ensuring phase-space volume conservation and bounded energy error.

• This explains their excellent long-term stability, especially for Hamiltonian systems like oscillators and planetary motion.

• Examples:

  • Symplectic Euler: first-order, simple, symplectic but not time-symmetric.

  • Leapfrog/Verlet: second-order, symplectic and time-symmetric, widely used in practice.

• Unlike standard explicit methods, symplectic schemes are not just about minimizing local error. They preserve the very geometry of Hamiltonian dynamics. Interesting connection to discrete differential geometry.