Dynamic PINNs (4): Solving Coupled ODEs with PINNs

This notebook extends the hybrid PINN approach to a double pendulum, a system of two coupled nonlinear ODEs. A single network with two outputs learns both angles simultaneously, with per-angle normalisation and two ODE residuals enforced in the physics loss. Curriculum training over five progressive time windows and a reduced physics loss weight help the network learn the oscillatory dynamics accurately. The equations of motion, the first-order conversion for numerical integration, and the training setup are walked through step by step.

pytorch
SciML
PINN
Author

Mei-Chin Pang

Published

May 29, 2026

Coupled ODE for a Double Pendulum

The Equations of Motion (letting \delta = \theta_2 - \theta_1)

\boxed{(m_1+m_2)L_1\frac{d^2\theta_1}{dt^2} + m_2 L_2 \frac{d^2\theta_2}{dt^2}\cos\delta - m_2 L_2 \left(\frac{d\theta_2}{dt}\right)^2\sin\delta + (m_1+m_2)g\sin\theta_1 = 0}

\boxed{L_2\frac{d^2\theta_2}{dt^2} + L_1\frac{d^2\theta_1}{dt^2}\cos\delta + L_1\left(\frac{d\theta_1}{dt}\right)^2\sin\delta + g\sin\theta_2 = 0}

These are the Lagrangian equations of motion for a double pendulum: two rigid rods of lengths L_1, L_2 and point masses m_1, m_2 connected end-to-end, swinging under gravity g. The angles \theta_1 and \theta_2 are measured from the vertical for the upper and lower pendulum respectively, and the shorthand \delta = \theta_2 - \theta_1 is the relative angle between the two rods.

Equation 1 (upper pendulum)

Term Meaning
(m_1+m_2)L_1(d^2\theta_1/dt^2) Angular acceleration of the upper rod, weighted by the total mass because the lower pendulum hangs from its tip.
m_2 L_2 (d^2\theta_2/dt^2)\cos\delta Coupling from the lower pendulum’s angular acceleration, projected along the upper rod’s tangent via \cos\delta.
-m_2 L_2 (d\theta_2/dt)^2\sin\delta Centripetal force that the swinging lower mass exerts on the pivot, directed perpendicular to the upper rod via \sin\delta.
(m_1+m_2)g\sin\theta_1 Gravitational restoring torque on the upper rod.

Equation 2 (lower pendulum)

Term Meaning
L_2(d^2\theta_2/dt^2) Angular acceleration of the lower rod (mass m_2 divides out).
L_1(d^2\theta_1/dt^2)\cos\delta Coupling from the upper pendulum’s angular acceleration, projected along the lower rod.
L_1(d\theta_1/dt)^2\sin\delta Centripetal force from the upper rod’s rotation acting on the lower pivot. Note the sign is positive (opposite to Equation 1) because the geometry reverses: here the upper rod pulls the lower pivot outward.
g\sin\theta_2 Gravitational restoring torque on the lower rod.

Why the system is coupled

Each equation contains d^2\theta/dt^2 of the other pendulum, so neither angle can be solved independently. The \cos\delta and \sin\delta terms also make the coupling nonlinear, since \delta itself depends on both angles. This is what makes the double pendulum famously chaotic for large amplitudes.

NoteUnderstanding the coupled ODE system

(a) What does “coupled” mean in this context?

“Coupled” means the equation of motion for each arm contains terms involving the motion of the other arm. We cannot solve for \theta_1(t) independently of \theta_2(t), the two ODEs must be solved simultaneously. Physically, this reflects the fact that the force the top arm exerts on the pivot depends on what the bottom arm is doing, and vice versa. The coupling is what makes the double pendulum dramatically more complex than two independent single pendulums.

(b) Why does the acceleration of the bottom arm appear in the top arm’s ODE?

The bottom arm is attached to the tip of the top arm, so any acceleration of the bottom arm creates an inertial reaction force on the top arm (Newton’s third law). Specifically, m_2 L_2 (d^2\theta_2/dt^2) \cos\delta is the component of the bottom arm’s linear acceleration that acts along the top arm’s radial direction. Even if the top arm were held fixed, the swinging bottom arm would exert a torque on it.

(c) With \theta_1(0) = 45°, \theta_2(0) = 0°, do you expect regular or chaotic motion?

At these relatively small initial angles (45° and 0°), the motion is expected to be regular (quasi-periodic) rather than chaotic. Double pendulum chaos typically requires larger initial angles – generally when the bottom arm has enough energy to flip over the top pivot. With \theta_1(0) = 45° and \theta_2(0) = 0°, the system has moderate energy and should exhibit complex but non-chaotic oscillatory motion over the 5-second window.

Initial conditions

\theta_1(0) = \pi/4, \quad \theta_2(0) = 0, \quad \frac{d\theta_1}{dt}(0) = \frac{d\theta_2}{dt}(0) = 0

The upper pendulum starts displaced at 45 degrees while the lower pendulum hangs vertically, and both start from rest. Because \theta_1(0) is moderately large (\pi/4 \approx 0.785 rad), the \sin\theta nonlinearity is significant from the outset, and the motion will not be well approximated by the linearised (small-angle) equations.

Loss Terms for Coupled-ODEs

\mathcal{L}_\text{IC} = \bigl(\hat\theta_1(0) - \tfrac{\pi}{4}\bigr)^2 + \bigl(\hat\theta_2(0)\bigr)^2 + \left(\frac{d\hat\theta_1}{dt}(0)\right)^2 + \left(\frac{d\hat\theta_2}{dt}(0)\right)^2

\mathcal{L}_\text{phys} = \frac{1}{N_\text{phys}}\sum_i\bigl[r_1(t_i)^2 + r_2(t_i)^2\bigr]

\qquad \mathcal{L}_\text{data} = \frac{1}{N_\text{data}}\sum_j\bigl[(\tilde{\hat\theta}_1 - \tilde\theta_{1}^\text{exp})^2 + (\tilde{\hat\theta}_2 - \tilde\theta_{2}^\text{exp})^2\bigr]

The total loss is \mathcal{L_\text{total}} = \lambda_\text{IC}\,\mathcal{L}_\text{IC} + \lambda_\text{phys}\,\mathcal{L}_\text{phys} + \lambda_\text{data}\,\mathcal{L}_\text{data}, combining three terms that each constrain a different aspect of the solution. The hat notation (\hat\theta) denotes the network’s prediction, and the tilde (\tilde\theta) denotes a normalised quantity.

  • IC loss (\mathcal{L}_\text{IC}): penalises deviation from the four known initial values:

    • \theta_1(0) = \pi/4,
    • \theta_2(0) = 0,
    • d\theta_1/dt(0) = 0,
    • d\theta_2/dt(0) = 0.

  • Unlike the single-pendulum case (two conditions), the coupled system requires four IC terms because there are now two second-order ODEs, each needing an initial displacement and an initial velocity.

  • Physics loss (\mathcal{L}_\text{phys}): the mean squared ODE residual evaluated at N_\text{phys} collocation points. Here r_1 and r_2 are the left-hand sides of the two boxed equations above. Both residuals are summed at each collocation point, so the network is simultaneously forced to satisfy both coupled ODEs.

  • Data loss (\mathcal{L}_\text{data}): the mean squared error between the network’s normalised predictions and the normalised experimental observations, summed over both angles at N_\text{data} data points. Each angle is normalised with its own statistics (\mu_1, \sigma_1 for \theta_1 and \mu_2, \sigma_2 for \theta_2) to ensure both contributions are on the same scale.

NoteUnderstanding the loss design for coupled ODEs

(a) Why use one network with two outputs rather than two separate networks?

A single shared network allows the hidden layers to learn shared representations of the time-domain dynamics that are useful for predicting both angles. Since \theta_1 and \theta_2 are physically coupled – their dynamics are intertwined – a shared network can exploit the correlation between the two outputs, potentially learning more efficiently. Two completely separate networks would have no mechanism to share information, and any physical coupling between the outputs would have to be captured indirectly through the physics loss alone. A single network with two outputs is also more parameter-efficient and easier to train as one optimisation problem.

(b) Why must both residuals r_1 and r_2 be included in the physics loss?

The two ODEs are coupled:

r_1 contains \frac{d^2\theta_2}{dt^2} terms and r_2 contains \frac{d^2\theta_1}{dt^2} terms.

If we only enforced r_1 = 0, the network would be constrained to satisfy the first equation but would be completely free to predict any \theta_2(t) that happens to make r_1 small, even if it grossly violates the second equation. The resulting \theta_2 prediction could be physically meaningless. Both residuals must be simultaneously driven to zero for the coupled physics to be correctly enforced.

(c) What could go wrong if we normalised both angles with the same constants?

\theta_1 (top arm, starting at 45°) and \theta_2 (bottom arm, starting at 0° but responding to the top arm’s motion) generally have different amplitude ranges and different statistical distributions. If we applied \mu_1, \sigma_1 to both, the normalised \theta_2 values would not have zero mean and unit variance – they could be systematically biased or have a very different scale. This would make the data loss contribution from \theta_2 disproportionately large or small compared to \theta_1, potentially causing the network to focus on one angle at the expense of the other.

import time
import warnings

import matplotlib.pyplot as plt
import numpy as np
import torch
import torch.nn as nn
from scipy.integrate import odeint
from tqdm import tqdm

# Suppress Intel Iris Xe "not officially supported" warning
warnings.filterwarnings(
    "ignore",
    message=".*detected GPU.*not officially supported by PyTorch XPU.*")

# Define the GPU type ("cuda" for NVIDIA GPUs, "xpu" for Intel GPUs)
gpu_type = "xpu"

if gpu_type == "cuda":
    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
elif gpu_type == "xpu":
    device = torch.device("xpu" if torch.xpu.is_available() else "cpu")
else:
    raise ValueError("Invalid GPU type. Use 'cuda' or 'xpu'.")

print(f"Using device: {device}")
Using device: cpu

Numerical Ground Truth

Physical parameters

Both pendulum arms have equal length (L_1 = L_2 = 1 m) and equal mass (m_1 = m_2 = 1 kg). The time window is t \in [0, 5] s with 1,000 evaluation points. The initial state vector is [\theta_1, \theta_2, \omega_1, \omega_2] = [\pi/4, 0, 0, 0], so the upper arm starts at 45° and the lower arm hangs vertically, both from rest.

Original second-order equations

\boxed{(m_1+m_2)L_1\frac{d^2\theta_1}{dt^2} + m_2 L_2 \frac{d^2\theta_2}{dt^2}\cos\delta - m_2 L_2 \left(\frac{d\theta_2}{dt}\right)^2\sin\delta + (m_1+m_2)g\sin\theta_1 = 0}

\boxed{L_2\frac{d^2\theta_2}{dt^2} + L_1\frac{d^2\theta_1}{dt^2}\cos\delta + L_1\left(\frac{d\theta_1}{dt}\right)^2\sin\delta + g\sin\theta_2 = 0}

Converting to a first-order system

odeint requires first-order ODEs, but the double pendulum equations are second-order. We introduce the angular velocities \omega_1 = d\theta_1/dt and \omega_2 = d\theta_2/dt as new variables, turning two second-order equations into four first-order equations.

Step 1: Define two new variables:

\omega_1 \equiv \frac{d\theta_1}{dt}, \qquad \omega_2 \equiv \frac{d\theta_2}{dt}

This immediately gives two of the four first-order equations:

\frac{d\theta_1}{dt} = \omega_1, \qquad \frac{d\theta_2}{dt} = \omega_2

Step 2: Since d\omega_1/dt = d^2\theta_1/dt^2 and d\omega_2/dt = d^2\theta_2/dt^2, the two original second-order equations become equations for d\omega_1/dt and d\omega_2/dt. Substituting (d^2\theta/dt^2) \to (d\omega/dt) and (d\theta/dt) \to \omega into the boxed equations:

(m_1+m_2)L_1\,\frac{d\omega_1}{dt} + m_2 L_2\,\frac{d\omega_2}{dt}\cos\delta - m_2 L_2\,\omega_2^2\sin\delta + (m_1+m_2)g\sin\theta_1 = 0

L_2\,\frac{d\omega_2}{dt} + L_1\,\frac{d\omega_1}{dt}\cos\delta + L_1\,\omega_1^2\sin\delta + g\sin\theta_2 = 0

Step 3: These two equations are linear in the unknowns d\omega_1/dt and d\omega_2/dt (everything else such as \theta_1, \theta_2, \omega_1, \omega_2 is known at the current time step). Writing them as a matrix system:

\begin{bmatrix} (m_1+m_2)L_1 & m_2 L_2\cos\delta \\ L_1\cos\delta & L_2 \end{bmatrix} \begin{bmatrix} \frac{d\omega_1}{dt} \\ \frac{d\omega_2}{dt} \end{bmatrix} = \begin{bmatrix} m_2 L_2\,\omega_2^2\sin\delta - (m_1+m_2)g\sin\theta_1 \\ -L_1\,\omega_1^2\sin\delta - g\sin\theta_2 \end{bmatrix}

Solving this 2 \times 2 system (e.g. by substitution or Cramer’s rule) gives explicit expressions for d\omega_1/dt and d\omega_2/dt, which is what the code computes via den1, den2 and the numerator expressions.

Result: The full first-order system passed to odeint is:

\frac{d}{dt}\begin{bmatrix}\theta_1 \\ \theta_2 \\ \omega_1 \\ \omega_2\end{bmatrix} = \begin{bmatrix}\omega_1 \\ \omega_2 \\ \frac{d\omega_1}{dt}(\theta_1,\theta_2,\omega_1,\omega_2) \\ \frac{d\omega_2}{dt}(\theta_1,\theta_2,\omega_1,\omega_2)\end{bmatrix}

The state vector is [\theta_1, \theta_2, \omega_1, \omega_2] and pendulum_system returns [d\theta_1/dt,\; d\theta_2/dt,\; d\omega_1/dt,\; d\omega_2/dt].

Solving for \frac{d^2\theta_1}{dt^2} and \frac{d^2\theta_2}{dt^2}

The two boxed equations each contain both (d^2\theta_1/dt^2) and (d^2\theta_2/dt^2), so they must be solved simultaneously as a 2 \times 2 linear system in the accelerations. The code does this algebraically:

  • delta = \theta_2 - \theta_1: the relative angle between the two arms.
  • den1 = (m_1 + m_2)L_1 - m_2 L_1 \cos^2\delta: the determinant-like denominator that arises from eliminating {d^2\theta_2}/{dt^2} to solve for {d^2\theta_1}/{dt^2}. It can never be zero for physical parameters since \cos^2\delta \leq 1 and m_1 > 0.
  • den2 = (L_2 / L_1) \cdot den1: the corresponding denominator for {d^2\theta_2}/{dt^2}, obtained by eliminating {d^2\theta_1}/{dt^2}.
  • The numerators collect the gravitational (g\sin\theta), centripetal (\omega^2\sin\delta), and coupling (\cos\delta) terms after the elimination.

Output

The returned solution array has shape (1000, 4), where columns 0-3 hold \theta_1, \theta_2, \omega_1, \omega_2 at each time point.

g       = 9.81
L1, L2  = 1.0, 1.0
m1, m2  = 1.0, 1.0

t_start, t_end = 0, 5
t = np.linspace(t_start, t_end, 1000)

# Initial conditions: [theta1, theta2, omega1, omega2]
initial_state = [45 * np.pi / 180, 0.0, 0.0, 0.0]

# Define the system of ODEs for the double pendulum
def pendulum_system(states, t, g, L1, L2, m1, m2):
    t1, t2, w1, w2 = states
    dt1 = w1
    dt2 = w2
    delta = t2 - t1
    den1  = (m1 + m2) * L1 - m2 * L1 * np.cos(delta)**2
    dw1   = (m2 * L1 * w1**2 * np.sin(delta) * np.cos(delta) +
             m2 * g * np.sin(t2) * np.cos(delta) +
             m2 * L2 * w2**2 * np.sin(delta) -
             (m1 + m2) * g * np.sin(t1)) / den1
    den2  = (L2 / L1) * den1
    dw2   = (-m2 * L2 * w2**2 * np.sin(delta) * np.cos(delta) +
             (m1 + m2) * g * np.sin(t1) * np.cos(delta) -
             (m1 + m2) * L1 * w1**2 * np.sin(delta) -
             (m1 + m2) * g * np.sin(t2)) / den2
    return [dt1, dt2, dw1, dw2]

solution = odeint(
    pendulum_system,
    initial_state,
    t,
    args=(g, L1, L2, m1, m2))
theta1 = solution[:, 0]
theta2 = solution[:, 1]
omega1 = solution[:, 2]
omega2 = solution[:, 3]

print(f"Numerical solution computed over t ∈ [{t_start}, {t_end}] s")
print(f"theta1 range: [{theta1.min():.3f}, {theta1.max():.3f}] rad")
print(f"theta2 range: [{theta2.min():.3f}, {theta2.max():.3f}] rad")
Numerical solution computed over t ∈ [0, 5] s
theta1 range: [-0.427, 0.785] rad
theta2 range: [-1.144, 0.622] rad

Displacement and Phase Portrait Plots

  • Displacement plot (\theta_1, \theta_2 vs t): shows both angular displacements over 5 seconds. Unlike the single damped pendulum (which decays monotonically), there is no damping here, so neither amplitude decays. The two curves oscillate with different frequencies and amplitudes because energy transfers back and forth between the arms through the coupling terms. \theta_1 starts at \pi/4 and \theta_2 starts at 0, but the lower arm quickly picks up motion from the upper arm. The resulting pattern is quasi-periodic – it looks irregular but repeats (approximately) over a longer timescale.

  • Phase portrait (\frac{d\theta}{dt} vs \theta for each arm): plots angular velocity against displacement with time removed. For a single undamped pendulum this would be a closed ellipse; here, the coupling between the two arms distorts the trajectories. Each arm traces a complex, non-repeating loop because its motion is modulated by the other arm. The trajectories remain bounded (no runaway growth), confirming that the system is in the regular (non-chaotic) regime at these initial conditions.

NoteInterpreting the displacement and phase portrait plots

(a) Why does \theta_2 have larger amplitude than \theta_1?

\theta_2 (bottom arm) typically has larger amplitude than \theta_1 (top arm). This makes physical sense, the bottom arm is free to swing independently and can amplify the motion it receives from the top arm, similar to how the tip of a whip moves faster than the handle. Energy is transferred from the top arm to the bottom arm through the coupling, and the smaller effective restoring torque at the bottom pivot (which is itself moving) allows larger angular excursions.

(b) Why does Arm 1’s phase portrait look smoother than Arm 2’s?

Arm 1’s phase portrait looks more ordered, it traces a relatively smooth, approximately elliptical or rosette-like curve. Arm 2’s phase portrait is more complex and irregular in shape. This is because \theta_1 is constrained by its connection to the fixed pivot and the inertia of the whole system, giving it a smoother, more slowly varying trajectory. \theta_2 responds to the irregular forcing from the moving top arm, producing a more chaotic-looking phase curve even in the non-chaotic regime.

(c) Why is there no amplitude decay to rest?

This model has no damping term, there is no b\frac{d\theta}{dt} dissipation in either equation of motion. Without energy dissipation, the total mechanical energy of the system is conserved, and the pendulum oscillates indefinitely. In the previous single-pendulum example, the b/m damping term continuously removed energy, causing decay to rest. To add damping to the double pendulum, we would need to add -c_1\frac{d\theta_1}{dt} and -c_2\frac{d\theta_2}{dt} terms to the respective equations of motion.

# ---------------------------------------------------------------------------
#  Plot 1: Displacement plot 
plt.figure(figsize=(9, 4))

# Plot the first arm's angle with a solid line
plt.plot(
    t,
    theta1,
    label=r'$\theta_1$ — Top Arm')

# Plot the second arm's angle with a dashed line
plt.plot(
    t,
    theta2,
    label=r'$\theta_2$ — Bottom Arm',
    linestyle='--')
plt.xlabel('Time (s)', fontsize=12)
plt.ylabel('Angular Displacement (rad)', fontsize=12)
plt.title('Double Pendulum — Numerical Solution', fontsize=13)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# ---------------------------------------------------------------------------
#  Plot 2: Phase portrait 
plt.figure(figsize=(9, 4))
plt.plot(
    theta1,
    omega1,
    label=r'Arm 1 ($\theta_1$ vs ${d\theta_1}/{dt}$)')
plt.plot(
    theta2,
    omega2,
    label=r'Arm 2 ($\theta_2$ vs ${d\theta_2}/{dt}$)',
    alpha=0.6)
plt.xlabel(
    'Angular Displacement (rad)', fontsize=12)
plt.ylabel('Angular Velocity (rad/s)', fontsize=12)
plt.title('2-DOF Phase Portrait', fontsize=13)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Synthetic Experimental Data

Sampling strategy

250 time points are drawn uniformly at random from [0, 5] s and sorted into chronological order. The ODE is solved at exactly these times to obtain the “true” displacement values for both arms. Using random (non-uniform) spacing mimics real sensor data, where measurements may arrive at irregular intervals.

Adding noise

Each displacement value is perturbed independently with Gaussian noise of \sigma = 0.005 rad (\approx 0.29°). Crucially, \theta_1 and \theta_2 receive separate noise draws – the measurement errors on the two sensors are independent, which is physically realistic (two separate angle sensors would have uncorrelated noise).

Overlay plots

Two subplots compare the 1,000-point numerical solution (coloured curve) with the 250 noisy data points (black scatter) for each arm separately:

  • Top subplot (\theta_1): the noisy data points should scatter tightly around the blue curve. Since \theta_1 has moderate amplitude (\sim 0.5 rad peak), \sigma = 0.005 rad is small relative to the signal – the noise is barely visible.

  • Bottom subplot (\theta_2): the noisy data points scatter around the red curve. \theta_2 has larger amplitude than \theta_1, so the same absolute noise level is an even smaller fraction of the signal here.

Plotting the two arms separately (rather than on one axis) makes it easier to verify that the noise level is realistic for each angle independently, since their amplitude ranges differ.

# For reproducibility
np.random.seed(42)

 # number of experimental observation points
N_exp_points = 250  

# noise standard deviation (rad)
exp_noise    = 0.005

# Generate synthetic experimental data by sampling the numerical 
# solution at random time points and adding noise
# sort to ensure time points are in order for plotting
t_exp_np = np.sort(np.random.rand(N_exp_points) * (t_end - t_start) + t_start)

# Integrate the system at the experimental time points to get the 
# "true" angles, then add noise
solution_exp  = odeint(
    pendulum_system,
    initial_state,
    t_exp_np,
    args=(g, L1, L2, m1, m2))
theta1_exp_np = solution_exp[:, 0] + exp_noise * np.random.randn(N_exp_points)
theta2_exp_np = solution_exp[:, 1] + exp_noise * np.random.randn(N_exp_points)

# ---------------------------------------------------------------------------
fig, axes = plt.subplots(2, 1, figsize=(9, 8))
axes[0].plot(
    t,
    theta1,
    color='blue',
    label=r'$\theta_1$ Numerical',
    alpha=0.7)
axes[0].scatter(
    t_exp_np,
    theta1_exp_np,
    color='black',
    s=10,
    label=r'$\theta_1$ Experimental')
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Angle (rad)')
axes[0].set_title('Top Pendulum')
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

axes[1].plot(
    t,
    theta2,
    color='red',
    label=r'$\theta_2$ Numerical',
    alpha=0.7)
axes[1].scatter(
    t_exp_np,
    theta2_exp_np,
    color='black',
    s=10,
    label=r'$\theta_2$ Experimental')
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Angle (rad)')
axes[1].set_title('Bottom Pendulum')
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

plt.suptitle(
    r'Synthetic Experimental Data vs Numerical Solution',
    fontsize=16)
plt.tight_layout()
plt.show()

Data Normalisation

This cell applies the same normalisation strategy used in the single-pendulum example, extended to handle two output variables.

Time normalisation

The raw experimental times t \in [0, 5] are mapped to [-1, +1]:

\tilde{t} = \frac{2(t - t_\text{min})}{t_\text{max} - t_\text{min}} - 1

Displacement normalisation (per-angle)

Each angle is standardised independently using its own experimental statistics:

\tilde{\theta}_1 = \frac{\theta_1 - \mu_{\theta_1}}{\sigma_{\theta_1}}, \qquad \tilde{\theta}_2 = \frac{\theta_2 - \mu_{\theta_2}}{\sigma_{\theta_2}}

Separate statistics are essential because \theta_1 and \theta_2 have different amplitude ranges and means. Using the same constants for both would bias one angle’s normalised values away from zero mean and unit variance, distorting the data loss balance between the two outputs.

Tensor conversion

  • t_exp_norm has shape (250, 1): a column vector of normalised times, matching the network’s single input.
  • theta_exp_norm has shape (250, 2): the two normalised angles are stacked column-wise using torch.stack(..., dim=1), matching the network’s two outputs. Column 0 holds \tilde{\theta}_1 and column 1 holds \tilde{\theta}_2.
# ----------------------------------------------------------------------------
# Time normalisation constants
t_min = t_start
t_max = t_end
dt    = t_max - t_min
t_exp_norm_np = 2 * (t_exp_np - t_min) / dt - 1.0

t_exp_norm    = torch.from_numpy(
    t_exp_norm_np).float().view(-1, 1).to(device)
print(f"t_exp_norm shape: {t_exp_norm.shape}")
print(f"t_exp_norm range (should be near [-1, +1]): " 
      + f"[{t_exp_norm.min():.3f}, {t_exp_norm.max():.3f}]")
print()
# ----------------------------------------------------------------------------
# Displacement normalisation: separate stats for each angle
theta1_mean = np.mean(theta1_exp_np)
theta1_std  = np.std(theta1_exp_np)
theta2_mean = np.mean(theta2_exp_np)
theta2_std  = np.std(theta2_exp_np)

# Standardise theta1 using its own statistics
theta1_exp_norm_np = (theta1_exp_np - theta1_mean) / theta1_std

# Standardise theta2 using its own statistics
theta2_exp_norm_np = (theta2_exp_np - theta2_mean) / theta2_std

# Convert to tensors, stack into shape (N, 2)
theta_exp_norm = torch.stack([
    torch.from_numpy(theta1_exp_norm_np).float(),
    torch.from_numpy(theta2_exp_norm_np).float()
], dim=1).to(device)

print(f"theta1_exp_norm mean (should be close to 0.0):  " 
      + f"{theta_exp_norm[:, 0].mean():.4f}")
print(f"theta2_exp_norm mean (should be close to 0.0):  " 
      + f"{theta_exp_norm[:, 1].mean():.4f}")
print(f"theta_exp_norm shape: " 
      f"{theta_exp_norm.shape}")
t_exp_norm shape: torch.Size([250, 1])
t_exp_norm range (should be near [-1, +1]): [-0.990, 0.974]

theta1_exp_norm mean (should be close to 0.0):  0.0000
theta2_exp_norm mean (should be close to 0.0):  0.0000
theta_exp_norm shape: torch.Size([250, 2])

Network Definition and Training Setup

Neural network architecture

FullyConnectedNN is a fully connected feedforward network with Tanh activations and two outputs, one for each pendulum angle. It takes a single normalised time input and returns normalised predictions for both angles simultaneously:

Component Code Description
fcs Linear(1, 64) + Tanh Input layer: maps the 1-D normalised time to 64 hidden units.
fch 3 × Linear(64, 64) + Tanh Hidden layers: three additional layers (N_LAYERS - 1 = 3) for greater representational capacity.
fce Linear(64, 2) Output layer: maps to two normalised displacement values [\tilde{\theta}_1, \tilde{\theta}_2] (no activation).

Normalisation helpers

  • normalise_t(t_tensor): maps raw physical time to [-1, +1] via

    \tilde{t} = \frac{2(t - t_\text{min})}{(t_\text{max} - t_\text{min})} - 1

  • denormalise_theta(theta_tensor): reverses the standardisation per-channel:

    • column 0 is denormalised with (\mu_{\theta_1}, \sigma_{\theta_1}) and
    • column 1 with (\mu_{\theta_2}, \sigma_{\theta_2}).
    • This ensures the physics residual is evaluated in correct physical units for each angle independently.

Curriculum training strategy

PINNs often struggle to learn oscillatory solutions over long time horizons because errors in the physics residual accumulate from early to late times. Curriculum training mitigates this by training on progressively longer time windows:

Stage Time window Epochs Purpose
1 [0, 1] s 5,000 Learn the initial transient accurately.
2 [0, 2] s 5,000 Extend the learned dynamics one second further.
3 [0, 3] s 5,000 Continue building outward from the established solution.
4 [0, 4] s 5,000 Near-full domain coverage.
5 [0, 5] s 5,000 Final stage covers the entire time domain.

At each stage, the collocation points and experimental data are restricted to the current window, but the network retains its weights from the previous stage. This means each stage refines and extends, rather than relearns the solution. The total training budget is 25,000 epochs across all five stages.

Training configuration

  • Architecture: 1 input → 4 hidden layers of 64 Tanh units → 2 outputs.
  • Collocation points: 1,000 uniformly spaced points within each stage’s time window, with requires_grad=True for autograd-based derivative computation. These are regenerated at the start of each curriculum stage.
  • IC point: a single point at t = 0 for the initial condition loss (four conditions: \theta_1, \theta_2, \dot\theta_1, \dot\theta_2).
  • Loss weights: \lambda_\text{IC} = 1, \lambda_\text{phys} = 0.1, \lambda_\text{data} = 1. The reduced physics weight prevents the large-magnitude ODE residual (which involves terms scaled by g \approx 9.81) from dominating the loss early in training, allowing the data loss to guide the network toward the correct waveform shape first.
  • Optimiser: Adam with LR = 10^{-3}.
  • LR scheduler: StepLR halves the learning rate every 5,000 epochs (aligning with curriculum stage boundaries), giving progressively finer convergence in later stages.
torch.manual_seed(42)

# -------------------------------------------------------------------
# PINN network architecture: fully connected with Tanh activations

class FullyConnectedNN(nn.Module):
    def __init__(self, N_INPUT, N_OUTPUT, N_HIDDEN, N_LAYERS):
        super().__init__()

        activation = nn.Tanh

        # Input layer + first hidden layer
        self.fcs = nn.Sequential(
            nn.Linear(N_INPUT, N_HIDDEN), activation())

        # Hidden layers (N_LAYERS - 1 because the first hidden
        # layer is already defined above)
        self.fch = nn.Sequential(*[
            nn.Sequential(
                nn.Linear(N_HIDDEN, N_HIDDEN), activation())
            for _ in range(N_LAYERS - 1)])

        # Output layer (linear, no activation)
        self.fce = nn.Linear(N_HIDDEN, N_OUTPUT)

    def forward(self, x):
        x = self.fcs(x)
        x = self.fch(x)
        x = self.fce(x)
        return x

# -------------------------------------------------------------------
# Normalisation helpers
def normalise_t(t_tensor):
    return 2.0 * (t_tensor - t_min) / dt - 1.0

def denormalise_theta(theta_tensor):
    # Reverses normalisation independently for each output channel
    th1 = theta_tensor[:, 0:1] * theta1_std + theta1_mean
    th2 = theta_tensor[:, 1:2] * theta2_std + theta2_mean
    return torch.cat([th1, th2], dim=1)

# -------------------------------------------------------------------
# Hyperparameters
n_input  = 1
n_output = 2       # predicts [theta1, theta2] simultaneously
n_hidden = 64
n_layers = 4
learning_rate = 1e-3

lambda_ic = 1
lambda_physics  = 0.1
lambda_data     = 1

# -------------------------------------------------------------------
# Curriculum training stages
# Train on progressively longer time windows so the network learns
# early-time dynamics first, then extends to the full domain.
curriculum_stages = [
    {"t_end": 1.0, "epochs": 5000},
    {"t_end": 2.0, "epochs": 5000},
    {"t_end": 3.0, "epochs": 5000},
    {"t_end": 4.0, "epochs": 5000},
    {"t_end": 5.0, "epochs": 5000},
]

# -------------------------------------------------------------------
# Instantiate the PINN model and set up training components
pinn = FullyConnectedNN(
    n_input, n_output, n_hidden, n_layers).to(device)

# Initial condition: t=0 -> theta1=45 deg, theta2=0 deg
t_ic = torch.tensor(
    [[0.0]], device=device, requires_grad=True)
N_phys  = 1000

optimiser = torch.optim.Adam(
    pinn.parameters(),
    lr=learning_rate)


scheduler = torch.optim.lr_scheduler.StepLR(
    optimiser, step_size=5000, gamma=0.5)


total_epochs = sum(s["epochs"] for s in curriculum_stages)
n_params = sum(p.numel() for p in pinn.parameters())
stages_summary = [
    (s['t_end'], s['epochs']) for s in curriculum_stages]

print(f"Network: {n_input} input -> "
      f"{n_layers}x{n_hidden} Tanh -> {n_output} outputs")
print(f"Parameters: {n_params}")
print(f"lambda_physics = {lambda_physics}  "
      f"(reduced to let data guide early training)")
print(f"Curriculum stages: {stages_summary}")
print(f"Total training epochs: {total_epochs}")
Network: 1 input -> 4x64 Tanh -> 2 outputs
Parameters: 12738
lambda_physics = 0.1  (reduced to let data guide early training)
Curriculum stages: [(1.0, 5000), (2.0, 5000), (3.0, 5000), (4.0, 5000), (5.0, 5000)]
Total training epochs: 25000

PINN Training Loop

The training loop iterates over five curriculum stages. At the start of each stage, 1,000 collocation points are regenerated over the current time window [0, t_\text{end}^\text{stage}], and only the experimental data falling within that window is selected. Within each stage, every epoch executes four parts then updates the network weights.

Part 1: Initial Condition (IC) Loss

  1. Normalise t = 0 and forward-pass through the network to get [\tilde{\theta}_1(0),\, \tilde{\theta}_2(0)].
  2. Denormalise both outputs to physical units.
  3. Compute [d\theta_1/dt,\, d\theta_2/dt] at t = 0 via torch.autograd.grad. Because the network has two outputs, torch.ones_like(thetas_ic) sums the gradients across both columns (each is a scalar at the single IC point), giving a (1, 2) Jacobian row.
  4. Penalise deviations from all four known initial values:

\mathcal{L}_\text{IC} = \bigl(\hat\theta_1(0) - \tfrac{\pi}{4}\bigr)^2 + \bigl(\hat\theta_2(0) - 0\bigr)^2 + \left(\frac{d\hat\theta_1}{dt}(0)\right)^2 + \left(\frac{d\hat\theta_2}{dt}(0)\right)^2

Part 2: Physics (Coupled ODE Residual) Loss

  1. Normalise the 1,000 collocation times (within the current stage’s window) and forward-pass to get normalised predictions.
  2. Denormalise to physical [\theta_1, \theta_2].
  3. Compute first and second derivatives for each angle separately via four calls to torch.autograd.grad (one per derivative per angle). All use create_graph=True so the computational graph extends through the derivatives for backpropagation.
  4. Evaluate both ODE residuals in physical space:

r_1 = (m_1+m_2)L_1\frac{d^2\theta_1}{dt^2} + m_2 L_2 \frac{d^2\theta_2}{dt^2}\cos\delta - m_2 L_2 \left(\frac{d\theta_2}{dt}\right)^2\sin\delta + (m_1+m_2)g\sin\theta_1

r_2 = L_2\frac{d^2\theta_2}{dt^2} + L_1\frac{d^2\theta_1}{dt^2}\cos\delta + L_1\left(\frac{d\theta_1}{dt}\right)^2\sin\delta + g\sin\theta_2

where \delta = \theta_2 - \theta_1.

\mathcal{L}_\text{phys} = \text{mean}(r_1^2) + \text{mean}(r_2^2)

Both residuals are included so the network is forced to satisfy both coupled ODEs simultaneously.

Part 3: Data Loss (in normalised space)

Pass the pre-normalised experimental times within the current stage’s window through the network. The MSE is computed across both angles:

\mathcal{L}_\text{data} = \text{mean}\!\bigl[(\tilde{\hat\theta}_1 - \tilde\theta_1^\text{exp})^2 + (\tilde{\hat\theta}_2 - \tilde\theta_2^\text{exp})^2\bigr]

Because each angle was normalised with its own statistics, the two columns contribute equally to the loss. As the curriculum window grows, more experimental data points are included.

Part 4: Total Loss and Weight Update

\mathcal{L_\text{total}} = \lambda_\text{IC}\,\mathcal{L}_\text{IC} + \lambda_\text{phys}\,\mathcal{L}_\text{phys} + \lambda_\text{data}\,\mathcal{L}_\text{data}

total_loss.backward() computes gradients, optimiser.step() updates weights via Adam, and scheduler.step() halves the learning rate every 5,000 epochs (coinciding with curriculum stage boundaries). All four loss components (total, IC, physics, data) are recorded at every epoch for the training loss plot.

total_loss_history = []
ic_loss_history = []
physics_loss_history = []
data_loss_history = []

start_training_time = time.time()

for stage_idx, stage in enumerate(curriculum_stages):
    t_end_stage   = stage["t_end"]
    stage_epochs  = stage["epochs"]

    # Build collocation points for the current time window
    t_phys = torch.linspace(
        t_start, t_end_stage, N_phys, device=device).view(-1, 1)
    t_phys.requires_grad_(True)

    # Select only experimental data within the current window
    mask       = t_exp_np <= t_end_stage
    t_exp_stage = t_exp_norm[mask]
    theta_exp_stage = theta_exp_norm[mask]

    print(f"Stage {stage_idx+1}/{len(curriculum_stages)} | "
          f"t ∈ [0, {t_end_stage}] s | {stage_epochs} epochs")

    for i in range(stage_epochs):
        
        optimiser.zero_grad()

        # --------------------------------------------------------------------
        # Part 1: Initial Condition Loss
        t_ic_norm      = normalise_t(t_ic)
        thetas_ic_norm = pinn(t_ic_norm)
        thetas_ic      = denormalise_theta(thetas_ic_norm)

        dthetas_dt_ic  = torch.autograd.grad(
            outputs=thetas_ic,
            inputs=t_ic,
            grad_outputs=torch.ones_like(thetas_ic),
            create_graph=True)[0]

        loss_theta_ic  = (
            torch.mean((thetas_ic[:, 0] - (45*torch.pi/180))**2)
            + torch.mean((thetas_ic[:, 1] - 0.0)**2))

        loss_dtheta_dt_ic = torch.mean(dthetas_dt_ic**2)
        total_loss_ic        = loss_theta_ic + loss_dtheta_dt_ic

        # --------------------------------------------------------------------
        # Part 2: Physics (Coupled ODE Residual) Loss
        t_phys_norm = normalise_t(t_phys)
        thetas_norm    = pinn(t_phys_norm)
        thetas         = denormalise_theta(thetas_norm)

        th1 = thetas[:, 0:1]
        th2 = thetas[:, 1:2]

        dth1_dt  = torch.autograd.grad(
            outputs=th1,
            inputs=t_phys,
            grad_outputs=torch.ones_like(th1),
            create_graph=True)[0]

        dth2_dt  = torch.autograd.grad(
            outputs=th2,
            inputs=t_phys,
            grad_outputs=torch.ones_like(th2),
            create_graph=True)[0]

        d2th1_dt2 = torch.autograd.grad(
            outputs=dth1_dt,
            inputs=t_phys,
            grad_outputs=torch.ones_like(dth1_dt),
            create_graph=True)[0]
        
        d2th2_dt2 = torch.autograd.grad(
            outputs=dth2_dt,
            inputs=t_phys,
            grad_outputs=torch.ones_like(dth2_dt),
            create_graph=True)[0]

        delta = th2 - th1

        # Compute the residuals of the ODEs (r1 and r2 should ideally be zero)
        r1 = (
            (m1+m2)*L1*d2th1_dt2
            + m2*L2*d2th2_dt2*torch.cos(delta)
            - m2*L2*dth2_dt**2*torch.sin(delta)
            + (m1+m2)*g*torch.sin(th1))

        r2 = (
            L2*d2th2_dt2
            + L1*d2th1_dt2*torch.cos(delta)
            + L1*dth1_dt**2*torch.sin(delta)
            + g*torch.sin(th2))

        loss_physics = torch.mean(r1**2) + torch.mean(r2**2)

        # --------------------------------------------------------------------
        # Part 3: Data Loss (normalised space)
        thetas_data_norm = pinn(t_exp_stage)
        loss_data        = torch.mean(
            (thetas_data_norm - theta_exp_stage)**2)

        # --------------------------------------------------------------------
        # Part 4: Total Loss
        total_loss = (
            lambda_ic * total_loss_ic
            + lambda_physics * loss_physics
            + lambda_data * loss_data)

        total_loss.backward()
        optimiser.step()
        scheduler.step()

        total_loss_history.append(total_loss.item())
        ic_loss_history.append(total_loss_ic.item())
        physics_loss_history.append(loss_physics.item())
        data_loss_history.append(loss_data.item())

        if i % 1000 == 0 or i == stage_epochs - 1:
            print(f"  Epoch {i:>5d}/{stage_epochs}  " 
                  f"Loss: {total_loss.item():.6f}")

end_training_time = time.time()
training_duration = end_training_time - start_training_time

print(f"\nTraining complete. Final loss: {total_loss_history[-1]:.6f}")
print(f"Total training time: {training_duration:.2f} seconds")
Stage 1/5 | t ∈ [0, 1.0] s | 5000 epochs
  Epoch     0/5000  Loss: 1.937241
  Epoch  1000/5000  Loss: 0.000584
  Epoch  2000/5000  Loss: 0.000591
  Epoch  3000/5000  Loss: 0.000609
  Epoch  4000/5000  Loss: 0.000634
  Epoch  4999/5000  Loss: 0.000471
Stage 2/5 | t ∈ [0, 2.0] s | 5000 epochs
  Epoch     0/5000  Loss: 3.019290
  Epoch  1000/5000  Loss: 0.000696
  Epoch  2000/5000  Loss: 0.000476
  Epoch  3000/5000  Loss: 0.000428
  Epoch  4000/5000  Loss: 0.000890
  Epoch  4999/5000  Loss: 0.000414
Stage 3/5 | t ∈ [0, 3.0] s | 5000 epochs
  Epoch     0/5000  Loss: 2.086826
  Epoch  1000/5000  Loss: 0.000753
  Epoch  2000/5000  Loss: 0.000500
  Epoch  3000/5000  Loss: 0.000434
  Epoch  4000/5000  Loss: 0.000397
  Epoch  4999/5000  Loss: 0.000375
Stage 4/5 | t ∈ [0, 4.0] s | 5000 epochs
  Epoch     0/5000  Loss: 2.054655
  Epoch  1000/5000  Loss: 0.033731
  Epoch  2000/5000  Loss: 0.001052
  Epoch  3000/5000  Loss: 0.000682
  Epoch  4000/5000  Loss: 0.000568
  Epoch  4999/5000  Loss: 0.000470
Stage 5/5 | t ∈ [0, 5.0] s | 5000 epochs
  Epoch     0/5000  Loss: 1.595923
  Epoch  1000/5000  Loss: 0.081953
  Epoch  2000/5000  Loss: 0.015572
  Epoch  3000/5000  Loss: 0.001536
  Epoch  4000/5000  Loss: 0.001006
  Epoch  4999/5000  Loss: 0.000679

Training complete. Final loss: 0.000679
Total training time: 316.82 seconds

Evaluation and Plots

PINN displacement prediction

500 evenly spaced test points are generated over [0, 5] s, normalised, and passed through the trained network inside torch.no_grad() (no gradient tracking needed at inference). The two-column normalised output is denormalised back to physical angles for plotting.

Plot 1: Displacement comparison (stacked)

Two vertically stacked subplots, one per pendulum arm, each overlaying three elements:

  • PINN prediction (black solid line): the trained network’s angle \hat\theta_i(t).
  • Numerical solution (dashed line): the ground-truth odeint result. If training succeeded, the two curves should be nearly indistinguishable.
  • Experimental data (coloured scatter): the 250 noisy observations used for training.

Plot 2: Training loss breakdown

Four loss curves on a log scale over 25,000 epochs (five curriculum stages of 5,000 epochs each):

  • Total loss (black): the weighted sum \lambda_\text{IC}\,\mathcal{L}_\text{IC} + \lambda_\text{phys}\,\mathcal{L}_\text{phys} + \lambda_\text{data}\,\mathcal{L}_\text{data}.
  • IC loss (orange): deviation from the four initial conditions.
  • Physics loss (red): mean squared ODE residual across collocation points.
  • Data loss (blue): MSE between predicted and observed normalised angles.

Orange dashed vertical lines mark the LR reduction points (every 5,000 epochs, coinciding with curriculum stage boundaries). Each reduction should produce a visible change in convergence behaviour as the optimiser enters a finer learning regime.

# Evaluate PINN
t_test      = torch.linspace(t_start, t_end, 500, device=device).view(-1, 1)
t_test_norm = normalise_t(t_test)

with torch.no_grad():
    thetas_pred_norm = pinn(t_test_norm)
    thetas_pred      = denormalise_theta(thetas_pred_norm)

t_plot   = t_test[:, 0].detach().cpu().numpy()
th1_pinn = thetas_pred[:, 0].detach().cpu().numpy()
th2_pinn = thetas_pred[:, 1].detach().cpu().numpy()

# ----------------------------------------------------------------------------
# Plot 1: Displacement — both arms
fig, axes = plt.subplots(
    2, 1,
    figsize=(9, 10))

axes[0].plot(
    t_plot,
    th1_pinn,
    color='black',
    linewidth=2,
    label=r'PINN $\hat\theta_1$')

axes[0].plot(
    t,
    theta1,
    color='blue',
    linestyle='--',
    alpha=0.7,
    label=r'Numerical $\theta_1$')

axes[0].scatter(
    t_exp_np,
    theta1_exp_np,
    s=15,
    color='blue',
    alpha=0.8,
    label='Exp Data')

axes[0].set_xlabel(
    'Time (s)',
    fontsize=12)
axes[0].set_ylabel(
    'Angle (rad)',
    fontsize=12)
axes[0].set_title(
    r'Top Pendulum Dynamics ($\theta_1$)',
    fontsize=13)
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

axes[1].plot(
    t_plot,
    th2_pinn,
    color='black',
    linewidth=2,
    label=r'PINN $\hat\theta_2$')
axes[1].plot(
    t,
    theta2,
    color='red',
    linestyle='--',
    alpha=0.7,
    label=r'Numerical $\theta_2$')
axes[1].scatter(
    t_exp_np,
    theta2_exp_np,
    s=15,
    color='red',
    alpha=0.8,
    label='Exp Data')
axes[1].set_xlabel(
    'Time (s)', fontsize=12)
axes[1].set_ylabel('Angle (rad)', fontsize=12)
axes[1].set_title(r'Bottom Pendulum Dynamics ($\theta_2$)', fontsize=13)
axes[1].legend(fontsize=10); axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# ----------------------------------------------------------------------------
# Plot 2: Training loss
fig, ax = plt.subplots(figsize=(9, 5))

ax.plot(
    total_loss_history,
    marker='o',
    markersize=2,
    color='black',
    linewidth=3,
    label="Total Loss")
ax.plot(
    ic_loss_history,
    color='orange',
    linewidth=0.8,
    label="IC Loss")
ax.plot(
    physics_loss_history,
    color='red',
    linewidth=1.2,
    alpha=0.8,
    label="Physics Loss")
ax.plot(
    data_loss_history,
    color='blue',
    linewidth=1.2,
    alpha=0.8,
    label="Data Loss")


for ep in [5000, 10000, 15000]:
    plt.axvline(
        x=ep,
        color='orange',
        linestyle='--',
        alpha=0.5)

ax.axvline(
    x=5000,
    color='orange',
    linestyle='--',
    alpha=0.5,
    label='LR reductions')

ax.set_yscale('log')
ax.set_xlabel('Epoch', fontsize=12)
ax.set_ylabel('Total Loss (log scale)', fontsize=12)
ax.set_title('2-DOF Hybrid PINN Training Loss', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, which='both', linestyle='-', alpha=0.2)

plt.tight_layout()
plt.show()

NoteSummary
  1. A 2-DOF PINN uses a single network with two output neurons to simultaneously predict the motion of both pendulum arms from a single time input.
  2. The physics loss enforces two coupled ODE residuals (r_1 and r_2) at every collocation point — both must be satisfied for the solution to be physically valid.
  3. Per-output normalisation ensures each angle is scaled independently, preventing one variable from dominating the data loss.
  4. A deeper network (4 layers) and longer training (20,000 epochs) are required to capture the more complex, coupled dynamics.
  5. The PINN framework scales naturally from 1-DOF to multi-DOF — the key extensions are the output dimension, the number of ODE residuals, and the per-output normalisation.

References

  1. Francis Fernandes (2026). Mastering Dynamic PINNs.
  2. Murad, J. ‘The Double Pendulum: Equations of Motion & Lagrangian Mechanics’, The Engineered Mind. Available at: https://www.engineered-mind.com/engineering/double-pendulum-1/ (Accessed: 29 May 2026).