Prev
Unit 9 · Ordinary Differential Equations: Initial Value Problems
Systems, Higher-Order ODEs and Geometric Integration
60m100 XP
Next
🛰️Ordinary Differential Equations: Initial Value Problems

Systems, Higher-Order ODEs and Geometric Integration

First-order systems, events, Störmer–Verlet and symplectic integration, rigid-body attitude dynamics

Kristiania, 1907
The aurora borealis appears in a ring around the magnetic pole, not over it. Kristian Birkeland, in Norway, had suggested that it is caused by charged particles from the Sun guided by the Earth's magnetic field. His colleague, the mathematician Carl Størmer, set out to compute where such particles actually go. The equation of motion of a charged particle in a dipole field has no closed-form solution. Størmer integrated it numerically — by hand, with the help of students — computing hundreds of trajectories over years of work. For the integration he used a simple scheme that steps positions directly from the second-order equation using accelerations at the current position, with no velocities at all. In 1967 the French physicist Loup Verlet, simulating the motion of argon atoms in a computer, used the same formula and found it remarkably well behaved over very long runs. Molecular dynamics adopted it under his name. Only in the 1980s and 1990s did numerical analysts explain why: the method preserves the geometric structure of Hamiltonian mechanics, which almost every other integrator destroys.
Everything is a first-order system. Solvers accept y′=f(t,y)\mathbf y' = \mathbf f(t, \mathbf y). A second-order equation θ′′=−gLsin⁡θ\theta'' = -\tfrac gL\sin\theta becomes two first-order ones by introducing the velocity: θ′=ω\theta' = \omega, ω′=−gLsin⁡θ\omega' = -\tfrac gL\sin\theta. A spacecraft's translational state is six numbers; with attitude (a quaternion) and angular velocity it is thirteen. Write the state vector, write its derivative, and every method of this unit applies unchanged.
📐The large-amplitude pendulum
The small-angle period 2πL/g=2.00642\pi\sqrt{L/g} = 2.0064 s holds for small swings only. Integrating the full pendulum with solve_ivp and an event at each upward zero crossing gives periods of 2.0074 s at 5°, 2.0413 s at 30°, 2.3683 s at 90° and 4.8944 s at 170° — 2.44 times the small-angle value. Events let the solver locate such crossings accurately without a fine output grid.
💡Störmer–Verlet
For q′′=a(q)\mathbf q'' = \mathbf a(\mathbf q), with p\mathbf p the velocity, one step is a kick–drift–kick: pn+1/2=pn+h2a(qn),qn+1=qn+h pn+1/2,pn+1=pn+1/2+h2a(qn+1).\mathbf p_{n+1/2} = \mathbf p_n + \tfrac h2\mathbf a(\mathbf q_n), \qquad \mathbf q_{n+1} = \mathbf q_n + h\,\mathbf p_{n+1/2}, \qquad \mathbf p_{n+1} = \mathbf p_{n+1/2} + \tfrac h2\mathbf a(\mathbf q_{n+1}). Second order, one force evaluation per step (the last one is reused), time-reversible — and symplectic: the map it applies to (q,p)(\mathbf q, \mathbf p) preserves phase-space area exactly, as the true flow of any Hamiltonian system does. A symplectic method exactly conserves a 'shadow' energy that differs from the true one by O(h2)O(h^2). So its energy error oscillates but never drifts, over any number of steps.
Derivation
Verlet is the central second difference, and a conserved 'shadow' energy you can check
Störmer's form. Replace q′′\mathbf q'' by the second difference of lesson 3-3 at step nn:
qn+1−2qn+qn−1h2=a(qn)⟹qn+1=2qn−qn−1+h2a(qn).\frac{\mathbf q_{n+1} - 2\mathbf q_n + \mathbf q_{n-1}}{h^2} = \mathbf a(\mathbf q_n) \qquad\Longrightarrow\qquad \mathbf q_{n+1} = 2\mathbf q_n - \mathbf q_{n-1} + h^2\mathbf a(\mathbf q_n).
The kick–drift–kick form computes exactly this. Its drifts give qn+1−qn=hpn+1/2\mathbf q_{n+1} - \mathbf q_n = h\mathbf p_{n+1/2} and qn−qn−1=hpn−1/2\mathbf q_n - \mathbf q_{n-1} = h\mathbf p_{n-1/2}, and the two half-kicks around qn\mathbf q_n give pn+1/2=pn−1/2+ha(qn)\mathbf p_{n+1/2} = \mathbf p_{n-1/2} + h\mathbf a(\mathbf q_n). Subtract the drifts and substitute: the same formula. It is symmetric in time — swap n+1n + 1 and n−1n - 1 and nothing changes — which is part of why it does not drift.
A conserved quantity you can verify. Take the oscillator x′′=−xx'' = -x and the symplectic Euler of lesson 9-1: vn+1=vn−hxnv_{n+1} = v_n - hx_n, then xn+1=xn+hvn+1x_{n+1} = x_n + hv_{n+1}. As a map on (x,v)(x, v) it is a matrix,
[xn+1vn+1]=[1−h2h−h1][xnvn],det⁡=(1−h2)+h2=1,\begin{bmatrix} x_{n+1} \\ v_{n+1} \end{bmatrix} = \begin{bmatrix} 1 - h^2 & h \\ -h & 1 \end{bmatrix}\begin{bmatrix} x_n \\ v_n \end{bmatrix}, \qquad \det = (1 - h^2) + h^2 = 1,
so it preserves areas in the (x,v)(x, v) plane exactly — the defining property of a symplectic map in two dimensions. (Explicit Euler's matrix [1h−h1]\begin{bmatrix} 1 & h \\ -h & 1 \end{bmatrix} has determinant 1+h21 + h^2: every step inflates areas, and the energy, by that factor.) Expanding the new values shows that
E~=x2+v2−h xv\tilde E = x^2 + v^2 - h\,xv
takes exactly the same value after every step. It differs from twice the true energy by a term of order hh, so the true energy wobbles within a narrow band and can never drift away. Störmer–Verlet has a shadow energy of the same kind, closer to the true one by a further power of hh.
📐A thousand orbits
An orbit of eccentricity 0.5 with 100 steps per revolution. RK4 — fourth order, four evaluations per step — has relative energy errors of 3.0×10−43.0 \times 10^{-4}, 3.0×10−33.0 \times 10^{-3} and 3.0×10−23.0 \times 10^{-2} after 10, 100 and 1000 orbits: it grows linearly, and the orbit slowly decays. Störmer–Verlet — second order, one evaluation per step — has 1.1×10−21.1 \times 10^{-2} after 10 orbits and still 1.1×10−21.1 \times 10^{-2} after 1000. Over short times RK4 is far better; over long times the symplectic method wins, at a quarter of the cost. The orbit's orientation precesses slightly under Verlet, but its size and shape — set by energy and angular momentum — do not decay.
The mathematical difficulties of the theory of rotation arise chiefly from the want of geometrical illustrations and sensible images, by which we might fix the results of analysis in our minds.
– James Clerk Maxwell, On a Dynamical Top (1857)
💡Rigid-body attitude as a system
A spacecraft with principal moments of inertia I1<I2<I3I_1 < I_2 < I_3 and no external torque obeys Euler's equations: I1ω˙1=(I2−I3)ω2ω3,I2ω˙2=(I3−I1)ω3ω1,I3ω˙3=(I1−I2)ω1ω2.I_1\dot\omega_1 = (I_2 - I_3)\omega_2\omega_3, \quad I_2\dot\omega_2 = (I_3 - I_1)\omega_3\omega_1, \quad I_3\dot\omega_3 = (I_1 - I_2)\omega_1\omega_2. Spin about the axis of least or greatest inertia is stable. Spin about the intermediate axis is not: integrated for two minutes from a 1 rad/s spin with a 10−310^{-3} wobble, the spin component falls to almost zero and reverses sign three times — the body flips end over end, while kinetic energy and angular momentum stay constant to 10−1010^{-10}. Cosmonaut Vladimir Dzhanibekov filmed a wing nut doing exactly this aboard Salyut 7 in 1985. Explorer 1, the first US satellite, showed a related instability in 1958: spinning about its long, minimum-inertia axis, it dissipated energy through its flexible antennas and ended up tumbling about the maximum-inertia axis.
Derivation
Euler's equations, and why the middle axis is unstable
Where the equations come from. Without external torque the angular momentum is constant in inertial space. In axes fixed to the body (principal axes), L=(I1ω1,I2ω2,I3ω3)\mathbf L = (I_1\omega_1, I_2\omega_2, I_3\omega_3), and a vector constant in space appears to change in rotating axes at the rate −ω×L-\boldsymbol\omega \times \mathbf L:
L˙+ω×L=0.\dot{\mathbf L} + \boldsymbol\omega \times \mathbf L = \mathbf 0.
The first component of ω×L\boldsymbol\omega \times \mathbf L is ω2I3ω3−ω3I2ω2\omega_2I_3\omega_3 - \omega_3I_2\omega_2, so I1ω˙1=(I2−I3)ω2ω3I_1\dot\omega_1 = (I_2 - I_3)\omega_2\omega_3; the other two follow cyclically.
Spin about the middle axis. Let the body spin at rate Ω\Omega about axis 2, with a small wobble ω1,ω3\omega_1, \omega_3. To first order in the wobble, ω2≈Ω\omega_2 \approx \Omega stays constant and
ω˙1=I2−I3I1 Ω ω3,ω˙3=I1−I2I3 Ω ω1.\dot\omega_1 = \frac{I_2 - I_3}{I_1}\,\Omega\,\omega_3, \qquad \dot\omega_3 = \frac{I_1 - I_2}{I_3}\,\Omega\,\omega_1.
Differentiate the first and substitute the second:
ω¨1=(I2−I3)(I1−I2)I1I3 Ω2 ω1.\ddot\omega_1 = \frac{(I_2 - I_3)(I_1 - I_2)}{I_1I_3}\,\Omega^2\,\omega_1.
With I1<I2<I3I_1 < I_2 < I_3 both brackets are negative, so the coefficient is positive: ω1\omega_1 grows exponentially, as ekte^{kt} with k=Ω(I3−I2)(I2−I1)/(I1I3)k = \Omega\sqrt{(I_3 - I_2)(I_2 - I_1)/(I_1I_3)} — an eigenvalue with positive real part (lesson 6-1). The same calculation for spin about axis 1 or axis 3 gives a negative coefficient: ω′′=−k2ω\omega'' = -k^2\omega, a bounded wobble.
💡Where you will use this
Orbital Mechanics: long-term propagation of planetary systems and debris uses symplectic integrators (Wisdom–Holman). Spacecraft Attitude Control: Euler's equations plus quaternion kinematics, with renormalisation or Lie-group integrators to keep the quaternion unit-length. Structural dynamics: the Newmark average-acceleration method is closely related to Störmer–Verlet and conserves energy for linear systems. Molecular and particle simulations: Verlet everywhere.
📚Further reading
  • Ernst Hairer, Christian Lubich and Gerhard Wanner, 'Geometric numerical integration illustrated by the Störmer–Verlet method', Acta Numerica 12, 399–450 (2003) — the whole theory of this lesson through one method.
  • Benedict Leimkuhler and Sebastian Reich, Simulating Hamiltonian Dynamics (Cambridge University Press, 2004) — symplectic integrators, shadow energies and rigid-body integration.
  • Ernst Hairer, Christian Lubich and Gerhard Wanner, Geometric Numerical Integration (2nd ed., Springer, 2006) — the full reference.
↓But there is a problem
Many of the systems in this unit had something in common once linearised: a spacecraft near its target attitude, a structure vibrating by small amounts, an electrical network. Their equations are linear with constant coefficients, x′=Ax+Bu\mathbf x' = A\mathbf x + B\mathbf u — and for those, stepping forward approximately is unnecessary. The exact solution is known: an exponential, of a matrix. Computing that exponential sounds trivial. In 1978 two numerical analysts published a survey of the ways people did it, and called it Nineteen Dubious Ways.