How to Solve for the Motion of a Double Pendulum

Question

I want to plot the motion of a double pendulum with a spring in python. I need to plot the theta1, theta2, r, and their first derivatives. I have found my equations for the motion, which are second-order ODEs so I then converted them to first-order ODEs where x1=theta1, x2=theta1-dot, y1=theta2, y2=theta2-dot, z1=r, and z2=r-dot. Here is a picture of the double pendulum problem: enter image description here

Here is my code:

from scipy.integrate import solve_ivp
from numpy import pi, sin, cos, linspace

g = 9.806 #Gravitational acceleration
l0 = 1 #Natural length of spring is 1
k = 2 #K value for spring is 2
OA = 2 #Length OA is 2
m = 1 #Mass of the particles is 1

def pendulumDynamics1(t, x): #Function to solve for theta-1 double-dot
    x1 = x[0]
    x2 = x[1]
    y1 = y[0]
    y2 = y[1]
    z1 = z[0]
    z2 = z[1]

    Fs = -k*(z1-l0)
    T = m*(x2**2)*OA + m*g*cos(x1) + Fs*cos(y1-x1)

    x1dot = x2
    x2dot = (Fs*sin(y1-x1) - m*g*sin(x1))/(m*OA) # angles are in radians
    return [x1dot,x2dot]

def pendulumDynamics2(t, y): #Function to solve for theta-2 double-dot
    x1 = x[0]
    x2 = x[1]
    y1 = y[0]
    y2 = y[1]
    z1 = z[0]
    z2 = z[1]
    Fs = -k*(z1-l0)
    
    y1dot = y2
    y2dot = (-g*sin(y1) - (Fs*cos(y1-x1)*sin(x1))/m + g*cos(y1-x1)*sin(x1) - x2*z1*sin(x1))/z1
    return [y1dot,y2dot]

def pendulumDynamics3(t, z): #Function to solve for r double-dot (The length AB which is the spring)
    x1 = x[0]
    x2 = x[1]
    y1 = y[0]
    y2 = y[1]
    z1 = z[0]
    z2 = z[1]
    Fs = -k*(z1-l0)
    
    z1dot = z2
    z2dot = g*cos(y1) - Fs/m + (y2**2)*z1 + x2*OA*cos(y1-x1) - (Fs*(sin(y1-x1))**2)/m + g*sin(x1)*sin(y1-x1)
    return [z1dot,z2dot]

# Define initial conditions, etc
d2r = pi/180
x0 = [30*d2r, 0] # start from 30 deg, with zero velocity
y0 = [60*d2r, 0] # start from 60 deg, with zero velocity
z0 = [1, 0] #Start from r=1
t0 = 0
tf = 10

#Integrate dynamics, initial value problem
sol1 = solve_ivp(pendulumDynamics1,[t0,tf],x0,dense_output=True) # Save as a continuous solution
sol2 = solve_ivp(pendulumDynamics2,[t0,tf],y0,dense_output=True) # Save as a continuous solution
sol3 = solve_ivp(pendulumDynamics3,[t0,tf],z0,dense_output=True) # Save as a continuous solution

t = linspace(t0,tf,200) # determine solution at these times
dt = t[1]-t[0]
x = sol1.sol(t)
y = sol2.sol(t)
z = sol3.sol(t)

I have 3 functions in my code, each to solve for x, y, and z. I then use solve_ivp function to solve for x, and y, and z. The error in the code is:

`File "C:\Users\omora\OneDrive\Dokument\AERO 211\project.py", line 13, in pendulumDynamics1 y1 = y[0]

NameError: name 'y' is not defined`

I don't understand why it is saying that y is not defined, because I defined it in my functions.

before going in the mathematics stuff line by line...you didn't define any variable named y,x,z or assigned it to value ..instead you defined y0 , x0 ,z0 ... for any one understands ODEs better can tell if the following is true or not because I don't remember too much ...if x0 , y0 , z0 are same as y , x , z in the initial condition of course ...you can rename x0 y0 z0 to x y z for the code to run ....but I don't if it is mathematically true — Hussien Mostafa, Dec 05 '21 at 04:42
Just to see that you have still a lot to study, there are some posts here and elsewhere on the simple double pendulum, like https://stackoverflow.com/questions/51722052/odd-slowing-effect-on-a-double-pendulum-animation, https://stackoverflow.com/questions/65224923/i-want-to-have-the-pendulum-blob-in-my-double-pendulum, https://stackoverflow.com/questions/61044226/solving-a-system-of-2nd-order-differential-equations-from-sympy. With the spring you get even more terms and possibly hard singularities in the solution. — Lutz Lehmann, Dec 05 '21 at 09:59

score 1 · Answer 1 · answered Dec 05 '21 at 17:34

Your system is closed without friction, thus can be captured by the Lagrange or Hamiltonian formalism. You have 3 position variables, thus a 6-dimensional dynamical state, complemented either by the velocities or the impulses.

Let q_k be theta_1, theta_2, r, Dq_k their time derivatives and p_k the impulse variables to q_k, then the dynamics can be realized by

def DoublePendulumSpring(u,t,params):
    m_1, l_1, m_2, l_2, k, g = params
    q_1,q_2,q_3 = u[:3]
    p = u[3:]

    A = [[l_1**2*(m_1 + m_2), l_1*m_2*q_3*cos(q_1 - q_2), -l_1*m_2*sin(q_1 - q_2)], 
         [l_1*m_2*q_3*cos(q_1 - q_2), m_2*q_3**2, 0], 
         [-l_1*m_2*sin(q_1 - q_2), 0, m_2]]

    Dq = np.linalg.solve(A,p)

    Dq_1,Dq_2,Dq_3 = Dq
    T1 = Dq_2*q_3*sin(q_1 - q_2) + Dq_3*cos(q_1 - q_2)
    T3 = Dq_1*l_1*cos(q_1 - q_2) + Dq_2*q_3

    Dp = [-l_1*(m_2*Dq_1*T1 + g*(m_1+m_2)*sin(q_1)), 
          l_1*m_2*Dq_1*T1 - g*m_2*q_3*sin(q_2), 
          m_2*Dq_2*T3 + g*m_2*cos(q_2) + k*(l_2 - q_3) ]
    return [*Dq, *Dp]

For a derivation see the Euler-Lagrange equations and their connection to the Hamilton equations. You might get asked about such a derivation.

This, after suitable defining the parameter tuple and initial conditions, can be fed to odeint and produces a solution that can then be plotted, animated or otherwise examined. The lower bob traces a path like the one below, not periodic and not very deterministic. (The fulcrum and the arc of the upper bob are also inserted, but less interesting.)

score 0 · Answer 2 · answered Dec 05 '21 at 05:52

def pendulumDynamics1(t, x):
    x1 = x[0]
    x2 = x[1]
    y1 = y[0]
    y2 = y[1]
    z1 = z[0]
    z2 = z[1]

You only pass x as a parameter. The code inside the function has no idea what y and z refer to.

You will need to change the function call to also include those variables.

def pendulumDynamics1(t, x, y, z):

How to Solve for the Motion of a Double Pendulum

2 Answers2