Double Pendulum Simulator: Equations of Motion, Energy, and Chaos

Simulate the double pendulum with four integrators and read the energy error each one commits. The Lagrangian is non-separable, so plain Verlet does not apply. The derivation is on the page.

At a glance

Computes
Double pendulum motion under four methods, with each method's energy error
You supply
Two angles, angular speeds, masses, rod lengths, gravity, and timestep
Use when
You want a defensible trajectory or want to compare integrator error
Not for
Predicting the exact position far ahead in time, because chaos forbids it
Initial state
Rig
Methods

The frame rate never reaches an integrator.

Relative energy deviation

Press play to start the clock

Energy is measured against t = 0 on a non-zero system scale.

The reach circle is an equal-aspect witness. Stretching would turn it into an ellipse.
Simulated time
0.00 s
Total energy
-1.801452e-15 J
Upper angle
90.000 degrees
Lower angle
90.000 degrees
Mass-matrix condition number
6.854
Driver
Implicit midpoint

The canonical solve is well conditioned. An extreme legal ratio triggers a digit-loss warning and warns instead of refusing.

Paused at the build-time reference state.

Signed relative energy deviation against Simulated time (s); reference at Initial energy

0246810−0.0003−0.0002−0.00010.00000.00010.00020.0003Initial energyImplicit midpoint
Moderate worked example, 10 seconds at dt = 1/240 s
IntegratorOrderSymplecticStepsFunction evaluationsSolver iterations meanSubstepsAcceptedRejectedMax constraint residualMax deviationDrift rateSecular fractionEnergy shapeStatus
RK44 no2,4009,60009.727e-10measured from tracemeasured from tracemeasured from tracepaused
Implicit midpoint2 yes, to solver tolerance2,400measured livemeasured live05.406e-6measured from tracemeasured from tracemeasured from tracepaused
Dormand-Prince 5(4)5(4) no2,400adaptive0measured livemeasured live1.021e-11measured from tracemeasured from tracemeasured from tracepaused
RATTLE2 yes, to solver tolerance2,400measured live0measured live3.795e-5measured from tracemeasured from tracemeasured from tracepaused

Integrator energy comparison

Signed relative energy deviation against Simulated time (s); reference at Initial energy

0246810−0.000000−0.0000000.0000000.0000000.000000Initial energyRK4

Signed relative energy deviation against Simulated time (s); reference at Initial energy

0246810−0.000006−0.000004−0.000002−0.0000000.0000020.0000040.000006Initial energyImplicit midpoint

Signed relative energy deviation against Simulated time (s); reference at Initial energy

0246810−0.000000−0.0000000.0000000.0000000.000000Initial energyDormand-Prince 5(4)

Signed relative energy deviation against Simulated time (s); reference at Initial energy

0246810−0.00004−0.000020.000000.000020.00004Initial energyRATTLE

Finite-time sensitivity

State separation against Simulated time (s), logarithmic vertical axis

051015200.0000000.0000000.0000000.0000000.0000010.0000100.0001000.0010000.010000Scaled state separation

Finite-time Lyapunov estimate is computed from the scaled separation trace; it is not an asymptotic exponent.

Phase portrait and Poincare section

p1 (kg m squared per s) against theta1 (radians)

−10−50510−2−1012
Points drawn at 0.45 opacity because 1000 overlap.

p1 at crossing against theta1 at upward theta2 crossing

7.58.08.5−0.20.00.20.40.6
Export

This page shows the reference configuration and computed energy trace from build time. Running the animation needs JavaScript.

Why plain Verlet does not apply here

The kinetic term depends on the momenta and the angle difference. The Hamiltonian is not separable, so kick-drift-kick Verlet in angle coordinates is a category error.

What a symplectic method actually promises

A fixed-step symplectic method follows a nearby modified Hamiltonian. Its energy error can occupy a bounded band without being exactly zero.

How the drift rate is measured

The signed relative energy deviation is fitted against simulated time. The secular fraction distinguishes a ramp from an oscillating band and is this page's convention.

Chaos and sensitivity

The sensitivity figure is a finite-time estimate, not an asymptotic exponent. A right-angle start loses predictive digits rapidly even when energy behavior is excellent.

The double slit experiment is analytic. The kinetic theory simulator advances exact collision events.

H(theta,p) = 0.5 p^T M(theta)^(-1) p + V(theta) How?

How this is calculated

The mass matrix contains cos(theta1 - theta2). RK4, implicit midpoint and Dormand-Prince step the same canonical state. RATTLE uses Cartesian coordinates and constrains both rods.

Formula: H(theta,p) = 0.5 p^T M(theta)^(-1) p + V(theta)

Assumptions and model
  • Two point masses on rigid massless rods.
  • A fixed frictionless pivot, no drag, and no damping.
  • Angles are absolute and measured from the downward vertical.
  • Chaos limits trajectory prediction even when numerical error is small.

Sources

  1. Geometric numerical integration illustrated by the Stormer-Verlet method. Acta Numerica. Retrieved .
  2. Andersen, RATTLE. Journal of Computational Physics. Retrieved .
  3. Leimkuhler and Skeel, Symplectic numerical integrators in constrained Hamiltonian systems. Journal of Computational Physics. Retrieved .
  4. Dormand and Prince, A family of embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics. Retrieved .
  5. Benettin et al., Lyapunov characteristic exponents for smooth dynamical systems. Meccanica. Retrieved .