Simple pendulum modeled as a differential-algebraic equation (DAE) instead of a ODE in the angle

pendulum_dae.py

# -*- coding: utf-8 -*-
"""
Solve the equations for a pendulum modeled as
a DAE

    x' = u
    y' = v
    u' = λ x
    v' = λ y - gy
    λ' = 3g/L² v

The jacobian of this system is given by

        ⎡0  0  1   0   0⎤
        ⎢               ⎥
        ⎢0  0  0   1   0⎥
        ⎢               ⎥
        ⎢λ  0  0   0   x⎥
     J =⎢               ⎥
        ⎢0  λ  0   0   y⎥
        ⎢               ⎥
        ⎢         3⋅g   ⎥
        ⎢0  0  0  ───  0⎥
        ⎢           2   ⎥
        ⎣          L    ⎦
    

@author: Nicolás Guarín-Zapata
@date: July 2025
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp


def fun(t, X):
    g = 1.0
    L = 1.0
    x, y, u, v, λ = X
    dXdt = [u, v, λ*x, λ*y - g, 3*g/L**2 * v]
    return dXdt


def jac(t, X):
    g = 1.0
    L = 1.0
    x, y, u, v, λ = X
    J = np.array([
        [0, 0, 1, 0, 0],
        [0, 0, 0, 1, 0],
        [λ, 0, 0, 0, x],
        [0, λ, 0, 0, y],
        [0, 0, 0, 3*g/L**2, 0]])
    return J


def euler(fun, t, X0):
    h = t[1] - t[0]
    df0 = fun(0, X0)
    nvars = np.asarray(df0).shape[0]
    ntimes = t.shape[0]
    X = np.zeros((nvars, ntimes))
    X[:, 0] = X0
    for cont in range(1, ntimes):
        X[:, cont] = X[:, cont - 1] + h * np.asarray(fun(0, X[:, cont - 1]))
    return X

    

X0 = [np.sqrt(2)/2, np.sqrt(2)/2, 0, 0, np.sqrt(2)/2]
t0 = 0.0
tf = 10.0
ntimes = 1000000
t = np.linspace(t0, tf, ntimes)
# sol = solve_ivp(fun, [t0, tf], X0, t_eval=t)

X = euler(fun, t, X0)
x, y, u, v, λ = X
T = 0.5*(u**2 + v**2)
V = y
E = T + V


#%%
plt.figure()
# plt.plot(t, x, label="Vertical position")
# plt.plot(x, y)
plt.plot(t, T, label="Kinetic Energy")
plt.plot(t, V, label="Potential Energy")
plt.plot(t, T - V, label="Lagrangian")
plt.plot(t, E, label="Total Energy")
plt.xlabel("Time")
plt.legend()