Double Pendulum Simulation

A modular, HLS-ready C++ simulation of the double pendulum using the Störmer–Verlet (second-order finite difference) integrator. Includes a serial and parallel ensemble runner and a Python visualiser.

---

Results

Phase Space (chaotic attractor)

Butterfly Effect (exponential divergence)

Performance: Serial vs Parallel

Phase Space Animation

---

Physics

System Description

The double pendulum consists of two point masses $m_1$, $m_2$ suspended by rigid rods of lengths $L_1$, $L_2$. The generalised coordinates are the angles $\theta_1$, $\theta_2$ measured from the vertical.

Equations of Motion

Applying the Euler–Lagrange equations to the Lagrangian of the system yields a coupled nonlinear ODE. At each instant, the angular accelerations $\ddot{\theta}_1$ and $\ddot{\theta}_2$ satisfy the $2 \times 2$ linear system:

$$ (m_1 + m_2),L_1,\ddot{\theta}_1 + m_2,L_2\cos(\Delta\theta),\ddot{\theta}_2 = \varepsilon $$

$$ m_2,L_1\cos(\Delta\theta),\ddot{\theta}_1 + m_2,L_2,\ddot{\theta}_2 = \zeta $$

where $\Delta\theta = \theta_1 - \theta_2$ and the forcing terms are:

$$ \varepsilon = -m_2,L_2,\dot{\theta}_2^2\sin(\Delta\theta) - m_2,g\sin(\theta_2) $$

$$ \zeta = m_2,L_1,\dot{\theta}_1^2\sin(\Delta\theta) - m_2,g\sin(\theta_2) $$

Solved analytically via Cramer's rule. Letting:

$$ D = (m_1 + m_2),L_1 \cdot m_2,L_2 ;-; m_2^2,L_1,L_2\cos^2(\Delta\theta) $$

the accelerations are:

$$ \ddot{\theta}_1 = \frac{\varepsilon,m_2 L_2 - m_2 L_2\cos(\Delta\theta),\zeta}{D} $$

$$ \ddot{\theta}_2 = \frac{(m_1+m_2)L_1,\zeta - m_2 L_1\cos(\Delta\theta),\varepsilon}{D} $$

Chaotic Behaviour

The system is chaotic for most initial conditions near the inverted configuration ($\theta_1 \approx \pi$). Nearby trajectories diverge exponentially — the butterfly effect visible in the divergence plot above.

---

Discretization

Two schemes are applied per timestep, one for each variable type.

Angular Velocity — Forward Euler (first-order)

$$ \dot{\theta}_i^{(n+1)} = \dot{\theta}_i^{(n)} + \ddot{\theta}_i^{(n)} \cdot \Delta t $$

Angle — Störmer–Verlet (second-order central difference)

The central difference approximation of the second derivative,

$$ \ddot{\theta}_i \approx \frac{\theta_i^{(n+1)} - 2,\theta_i^{(n)} + \theta_i^{(n-1)}}{\Delta t^2} $$

rearranged to advance the position:

$$ \theta_i^{(n+1)} = 2,\theta_i^{(n)} - \theta_i^{(n-1)} + \Delta t^2,\ddot{\theta}_i^{(n)} $$

Bootstrapping

The Verlet recurrence requires $\theta^{(-1)}$ (position at $t = -\Delta t$). This is initialised from the initial velocity via a first-order backward step:

$$ \theta_i^{(-1)} = \theta_i^{(0)} - \Delta t \cdot \dot{\theta}_i^{(0)} $$

With zero initial velocities ($\dot{\theta}_i^{(0)} = 0$) this gives $\theta_i^{(-1)} = \theta_i^{(0)}$, which is a clean symmetric start.

Note on Symplecticity

Störmer–Verlet applied purely to positions is symplectic — it conserves a shadow Hamiltonian, bounding long-term energy drift. The forward Euler step on $\dot{\theta}$ breaks this guarantee weakly. A fully symplectic treatment would use the velocity Verlet form. For chaotic systems where short-time accuracy matters more than long-time energy conservation, the mixed scheme used here is standard.

---

Ensemble Parallelism

A single pendulum trajectory is inherently sequential (step $n+1$ depends on step $n$), so parallelism is achieved across an ensemble of $N = 8$ trajectories with perturbed initial conditions:

$$ \theta_1^{(0)}(i) = 3.14 + i \times 0.01, \quad i = 0, 1, \ldots, 7 $$

Each thread owns exactly one trajectory. There is no inter-thread communication ($\text{comm_ns} = 0$), so speedup is near-linear in $N$. This contrasts with a grid-based PDE solver (e.g. heat2d) where threads must exchange ghost rows at every timestep, capping practical speedup.

---

Project Structure

double-pendulum-eqn/
  assets/                         committed visualisation outputs
  source/
    include/
      config.hpp                  compile-time constants
      dp_types.hpp                PendulumState, PendulumPrev (plain-old-data)
      dp_solver.hpp               HLS-friendly core: compute_step(), run_simulation()
      dp_state_utils.hpp          SolverMetrics struct, CSV helpers
      dp_serial.hpp               serial ensemble solver declaration
      dp_parallel.hpp             parallel ensemble solver declaration
      runner.hpp                  RunnerResult, run_pendulum2d(), print_summary()
    src/
      dp_solver.cpp               physics kernel (no STL, no dynamic memory -- HLS-ready)
      dp_state_utils.cpp          ensure_dir, save_trajectory_csv, save_metrics_csv
      dp_serial.cpp               N_TRAJ trajectories sequentially, timed
      dp_parallel.cpp             N_TRAJ trajectories in parallel threads, timed
      runner.cpp                  orchestrates serial + parallel, prints summary
    double_pendulum_tb.cpp        main()
    testbench.cpp                 run_testbench() -> run_pendulum2d()
    visualize.py                  phase space, butterfly effect, timing plots
    run_viz.sh                    bootstrap Python venv and launch visualizer
    requirements.txt
    Makefile
    CMakeLists.txt

---

Build and Run

Requires clang++ with C++17 support.

cd source
make          # compile
make run      # compile and run
make clean    # remove binary and output/

Expected output:

[serial]   Running 8 trajectories x 10000 steps...
[parallel] Running 8 trajectories x 10000 steps (8 threads)...

===== double-pendulum ensemble run summary =====
  Trajectories : 8
  Steps/traj   : 10000
  Threads(par) : 8

  Serial     | total=... ms  compute=... ms  ...
  Parallel   | total=... ms  compute=... ms  ...

  Compute speedup (serial/parallel): ~8x
  (comm_ns = 0 -- no inter-thread data exchange needed)

---

Output Files

After running, the output/ directory contains:

output/
  serial/
    traj_0.csv .. traj_7.csv     columns: step, theta1, theta2
  parallel/
    traj_0.csv .. traj_7.csv
  serial_metrics.csv
  parallel_metrics.csv

Serial and parallel files for the same trajectory index are byte-for-byte identical. This is guaranteed: each trajectory is a fully independent sequential chain with no shared mutable state between threads — identical inputs and operation order produce identical floating-point results.

---

Visualisation

./run_viz.sh                    # interactive (requires display)
./run_viz.sh --save             # save PNGs and GIFs to output/plots/
./run_viz.sh --mode serial      # serial trajectories only
./run_viz.sh --no-anim --save   # static plots only

Output file	Description
phase_space.png	$\theta_1$ vs $\theta_2$ for all trajectories
divergence.png	$
timing_metrics.png	serial vs parallel compute time and speedup vs ideal
phase_animation_*.gif	animated phase-space trajectories

---

HLS Notes

The core physics kernel (dp_solver.cpp, dp_types.hpp) is written to be HLS-compatible:

• No STL containers or dynamic memory (std::vector, new, malloc absent)
• No HLS pragmas yet; these will be added in Vivado/Vitis HLS
• All state passed by reference as plain structs, allocated on the stack
• The ensemble runner (dp_serial, dp_parallel, runner) uses STL and is host-side infrastructure only

---

Parameters

Defined in include/config.hpp:

Constant	Value	Description
DP_G	9.8	gravitational acceleration ($m/s^2$)
DP_L1	2.0	length of pendulum 1 ($m$)
DP_L2	1.0	length of pendulum 2 ($m$)
DP_M1	1.0	mass of bob 1 ($kg$)
DP_M2	2.0	mass of bob 2 ($kg$)
DP_DT	0.001	timestep $\Delta t$ ($s$)
DP_N_STEPS	10000	steps per trajectory
DP_N_TRAJ	8	ensemble size $N$
DP_IC_STEP	0.01	per-trajectory $\theta_1$ perturbation