Linearise about the hanging position: the eigenproblem (K − ω²M0)v = 0 gives three frequencies. A 1 mrad nonlinear run (200 s, no damping) is analysed with an FFT; peaks should agree within 1 %.
Undamped run of 60 s from the current start angles with the selected integrator. Relative drift is |ΔE| divided by max(|E0|, 0.1 gΣμiiLi). RK4 at 0.25 ms should stay below 1e-8.
Method A: Benettin renormalisation (150 s). Method B: 48 independent two-trajectory pairs, no renormalisation, slope of ln d(t) in the exponential window. They should agree within 15 % when the motion is chaotic (λ > 0.3 1/s).
A planar triple pendulum has three angles θi measured from the downward vertical, point masses mi on massless rigid rods of length Li. With μij = Σk≥max(i,j) mk:
T = ½ ΣΣ μij Li Lj cos(θi−θj) θ'i θ'j V = −g Σ μii Li cos θiThe Euler–Lagrange equations give the mass-matrix form solved here at every step:
Σj μij Li Lj cos(θi−θj) θ''j = −Σj μij Li Lj sin(θi−θj) θ'j² − g μii Li sin θi − c θ'iIntegration is classical RK4 with a fixed step (default 0.25 ms). Energy is conserved only up to the integrator error; the semi-implicit Euler option is first order and drifts visibly. Joint friction is a simple viscous torque cθ'; real bearings also have dry friction and the rods have inertia, both ignored here.
Two orbits that start ε apart separate as d(t) ≈ ε eλt until the separation saturates at the size of the attractor. λ is the maximal Lyapunov exponent; 1/λ is the horizon after which a measurement error has grown by a factor e. The ghosts are copies displaced by ε (default 10−9 rad) in a random direction of angle space.
Method A (Benettin et al., 1980): follow a shadow orbit, every 0.05 s add ln(d/d0) and rescale the shadow back to d0; λ is the mean of those logs. Method B: fit the slope of ln d(t) over the growth window (100ε to 10−4) for many independent pairs and average. Distances use the 6-dimensional (θ, θ') Euclidean norm. Finite-time estimates fluctuate along an orbit, so the two numbers agree only statistically. Float64 round-off acts like an extra random ε of about 10−16, so prediction beyond roughly ln(1/10−16)/λ is meaningless even for the exact equations.
The Poincaré section records (θ2, θ'2) each time θ3 crosses zero upward: a closed curve means quasi-periodic motion, a scattered cloud means chaos. The flip-time map starts the pendulum at rest from a grid of (θ1, θ2) and colours each pixel by the time until the chosen arm first passes over the top (|θ|>π); the fractal boundaries are the signature of chaotic scattering. At low energy the motion is regular; the regime scan shows λ rising with energy, not strictly monotonic because regular islands persist at high energy (KAM theory).
Small oscillations decouple into three normal modes (the eigenproblem in Experiment 1). The frequencies f = ω/2π for three equal masses and lengths are 0.321, 0.755 and 1.250 Hz at L = 1 m on Earth.