Simulation of a Pendulum hanging on a spinning Disk

Can anybody get this code to run? I know, that it is very long and maybe not easy to understand, but what I am trying to do is to write a simulation for a problem, that I have already posted here: https://math.stackexchange.com/questions/4876146/pendulum-hanging-on-a-spinning-disk

I try to make a nice simulation, that would look like the one someone answered the linked question with. The picture in the answer is written in mathematica and I had no idea how to translate it. Hope you can help me finish this up.

There are two elements of code. One calculating the ODE second degree and one plotting it 3 times. When you plot the ODE, you can see, that the graph line is not doing, what it is supposed to. I don't know, where the mistake is but hopefully you can help.

Here are the two snippets:

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import display
from sympy.vector import CoordSys3D
from scipy.integrate import solve_ivp

def FindDGL():
    N = CoordSys3D('N')
    e = N.i + N.j + N.k
    t = sp.symbols('t')

    x = sp.symbols('x')
    y = sp.symbols('y')
    z = sp.symbols('z')

    x = sp.Function('x')(t)
    y = sp.Function('y')(t)
    z = sp.Function('z')(t)

    p = x*N.i + y*N.j + z*N.k

    m = sp.symbols('m')
    g = sp.symbols('g')
    r = sp.symbols('r')
    omega = sp.symbols('omega')
    q0 = sp.symbols('q0')
    A = sp.symbols('A')
    l = sp.symbols('l')
    
    
    xl = sp.symbols('xl')
    yl = sp.symbols('yl')
    zl = sp.symbols('zl')

    dpdt = sp.diff(x,t)*N.i + sp.diff(y,t)*N.j + sp.diff(z,t)*N.k
    
    #Zwang = ((p-q)-l)
    Zwang = (p.dot(N.i)**2*N.i +p.dot(N.j)**2*N.j +p.dot(N.k)**2*N.k - 2*r*(p.dot(N.i)*N.i*sp.cos(omega*t)+p.dot(N.j)*N.j*sp.sin(omega*t))-2*q0*(p.dot(N.k)*N.k) + r**2*(N.i*sp.cos(omega*t)**2+N.j*sp.sin(omega*t)**2)+q0**2*N.k) - l**2*N.i - l**2*N.j -l**2*N.k
   
    display(Zwang)

    dpdtsq = dpdt.dot(N.i)**2*N.i + dpdt.dot(N.j)**2*N.j + dpdt.dot(N.k)**2*N.k
    
    #La = 0.5 * m * dpdtsq - m * g * (p.dot(N.k)*N.k) + (ZwangA*A)

    L = 0.5 * m * dpdtsq + m * g * (p.dot(N.k)*N.k) - Zwang*A

    #display(La)
    display(L)
    
    Lx = L.dot(N.i)
    Ly = L.dot(N.j)
    Lz = L.dot(N.k)

    Elx = sp.diff(sp.diff(Lx,sp.Derivative(x,t)), t) + sp.diff(Lx,x)
    Ely = sp.diff(sp.diff(Ly,sp.Derivative(y,t)), t) + sp.diff(Ly,y)
    Elz = sp.diff(sp.diff(Lz,sp.Derivative(z,t)), t) + sp.diff(Lz,z)

    display(Elx)
    display(Ely)
    display(Elz)

    ZwangAV = (sp.diff(Zwang, t, 2))/2
    display(ZwangAV)
    ZwangA = ZwangAV.dot(N.i)+ZwangAV.dot(N.j)+ZwangAV.dot(N.k)
    display(ZwangA)

    Eq1 = sp.Eq(Elx,0)
    Eq2 = sp.Eq(Ely,0)
    Eq3 = sp.Eq(Elz,0)
    Eq4 = sp.Eq(ZwangA,0)

    LGS = sp.solve((Eq1,Eq2,Eq3,Eq4),(sp.Derivative(x,t,2),sp.Derivative(y,t,2),sp.Derivative(z,t,2),A))
    #display(LGS)
    #display(LGS[sp.Derivative(x,t,2)].free_symbols)
    #display(LGS[sp.Derivative(y,t,2)].free_symbols)
    #display(LGS[sp.Derivative(z,t,2)])
    XS = LGS[sp.Derivative(x,t,2)]
    YS = LGS[sp.Derivative(y,t,2)]
    ZS = LGS[sp.Derivative(z,t,2)]
    
    #t_span = (0, 10)
    dxdt = sp.symbols('dxdt')
    dydt = sp.symbols('dydt')
    dzdt = sp.symbols('dzdt')
    #t_eval = np.linspace(0, 10, 100)
    XSL = XS.subs({ sp.Derivative(y,t):dydt, sp.Derivative(z,t):dzdt, sp.Derivative(x,t):dxdt, x:xl , y:yl , z:zl})
    YSL = YS.subs({ sp.Derivative(y,t):dydt, sp.Derivative(z,t):dzdt, sp.Derivative(x,t):dxdt, x:xl , y:yl , z:zl})
    ZSL = ZS.subs({ sp.Derivative(y,t):dydt, sp.Derivative(z,t):dzdt, sp.Derivative(x,t):dxdt, x:xl , y:yl , z:zl})
    #display(ZSL.free_symbols)
    XSLS = str(XSL)
    YSLS = str(YSL)
    ZSLS = str(ZSL)
    replace = {"xl":"x","yl":"y","zl":"z","cos":"np.cos", "sin":"np.sin",}
    for old, new in replace.items():
        XSLS = XSLS.replace(old, new)
    for old, new in replace.items():
        YSLS = YSLS.replace(old, new)
    for old, new in replace.items():
        ZSLS = ZSLS.replace(old, new)
    return[XSLS,YSLS,ZSLS]


Result = FindDGL() 
print(Result[0])
print(Result[1])
print(Result[2])

Here is the second one:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D

def Q(t):
    omega = 1 
    return r * (np.cos(omega * t) * np.array([1, 0, 0]) + np.sin(omega * t) * np.array([0, 1, 0])) + np.array([0, 0, q0])

def equations_of_motion(t, state, r, omega, q0, l):
    x, y, z, xp, yp, zp = state
    dxdt = xp
    dydt = yp
    dzdt = zp
    dxpdt = dxdt**2*r*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dxdt**2*x/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dxdt*omega*r**2*np.sin(omega*t)*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dxdt*omega*r*x*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dydt**2*r*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dydt**2*x/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dydt*omega*r**2*np.cos(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dydt*omega*r*x*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dzdt**2*r*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dzdt**2*x/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + g*q0*r*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*q0*x/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*r*z*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + g*x*z/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*r**2*x*np.cos(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*r**2*y*np.sin(omega*t)*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*x**2*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*x*y*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2)
    dypdt = dxdt**2*r*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dxdt**2*y/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dxdt*omega*r**2*np.sin(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dxdt*omega*r*y*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dydt**2*r*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dydt**2*y/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dydt*omega*r**2*np.sin(omega*t)*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dydt*omega*r*y*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dzdt**2*r*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dzdt**2*y/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + g*q0*r*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*q0*y/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*r*z*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + g*y*z/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*r**2*x*np.sin(omega*t)*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*r**2*y*np.sin(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*x*y*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*y**2*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2)
    dzpdt = dxdt**2*q0/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dxdt**2*z/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dxdt*omega*q0*r*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dxdt*omega*r*z*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dydt**2*q0/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dydt**2*z/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dydt*omega*q0*r*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dydt*omega*r*z*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dzdt**2*q0/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dzdt**2*z/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*r**2*np.sin(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*r**2*np.cos(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*g*r*x*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*g*r*y*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*x**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*y**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*q0*r*x*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*q0*r*y*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*x*z*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*y*z*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2)
    return [dxdt, dydt, dzdt, dxpdt, dypdt, dzpdt]

r = 0.5  
omega = 1.2  
q0 =  1 
l = 1  
g = 9.81

#{x[0] == r, y[0] == x'[0] == y'[0] == z'[0] == 0, z[0] == q0 - l}
initial_conditions = [r, 0, 0, 0, 0, q0-l] 

tmax = 200
solp = solve_ivp(equations_of_motion, [0, tmax], initial_conditions, args=(r, omega, q0, l), dense_output=True)#, method='DOP853')

t_values = np.linspace(0, tmax, 1000)
p_values = solp.sol(t_values)
print(p_values.size)

d =0.5

Qx = [Q(ti)[0] for ti in t_values]
Qy = [Q(ti)[1] for ti in t_values]
Qz = [Q(ti)[2] for ti in t_values]

fig = plt.figure(figsize=(20, 16))
ax = fig.add_subplot(111, projection='3d')
ax.plot(p_values[0], p_values[1], p_values[2], color='blue')
ax.scatter(r, 0, q0-l, color='red')
ax.plot([0, 0], [0, 0], [0, q0], color='green') 

ax.plot(Qx, Qy, Qz, color='purple')
#ax.set_xlim(-d, d)
#ax.set_ylim(-d, d)
#ax.set_zlim(-d, d)
ax.view_init(30, 45)
plt.show()

Solution

Edited to include the direct solution in Cartesians, since this was the original direction taken by the OP. See the bottom of this answer.

Edited again to provide an alternative derivation of the Cartesian equations from Newton's Second Law (F=ma) rather than Lagrangian Mechanics. (This has the side benefit that it also finds the tension in the pendulum rod.)

Lagrangian Equations (i) Angle Coordinates

This method, related to, but not quite the same as, the replies in your Stack Exchange post is quite appealing because it treats your 2-degree-of-freedom problem with two independent coordinates, so not having to deal with the constraint (that the length of pendulum is fixed).

Let ϕ be the angle between the vertical plane containing the string and the x-axis. Let θ be the angle made by the string with the downward vertical. Then the coordinates of the bob (relative to the centre of the ring) are:

Differentiate with respect to time and we get velocity components

Summing, squaring and simplifying using trig formulae:

The Lagrangian (strictly, Lagrangian divided by mass, but mass would cancel in the analysis that follows) is

whence

From the Lagrangian equations for a conservative system

we get (after a LOT of algebra!) the key equations for our degrees of freedom:

The code sample below solves these using solve_ivp. Note that the denominator of the second equation means that you shouldn’t start with the pendulum vertical, as sin⁡θ would be 0.

Then your animation.

I use FuncAnimation from matplotlib.animation. The objective here is to update any elements of your plot (in this case the ends of the pendulum string) in each frame. The code plots an animation by default, but you can remove the subsequent comment to store it as a (rather large) graphics file.

If you want the trajectory rather than the animation then re-comment the commands at the bottom to choose plot_figure instead.

Note that the system is chaotic and very sensitive to the relative sizes of disk and pendulum and the angular velocity ω. If I do a dimensional analysis, the trajectory shape should be a function of g/ω^2 L (or the ratios of the periods of the disk and the stand-alone pendulum), R/L, and the initial angles.

CODE

import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D

g = 9.81

def plot_animation( t, qx, qy, qz, x, y, z ):
   plotInterval = 1
   fig = plt.figure( figsize=(4,4) )
   ax = fig.add_subplot(111, projection='3d' )
   ax.view_init(30, 45)
   ax.set_xlim( -1.0, 1.0 );   ax.set_ylim( -1.0, 1.0 );   ax.set_aspect('equal')

   ax.plot( qx, qy, qz, 'k' )                                            # ring
   a = ax.plot( [qx[0],x[0]], [qy[0],y[0]], [qz[0],z[0]], 'g' )          # pendulum string
   b = ax.plot( [x[0]], [y[0]], [z[0]], 'ro' )                           # pendulum bob
   def animate( i ):
      a[0].set_data_3d( [qx[i],x[i]], [qy[i],y[i]], [qz[i],z[i]] )       # update anything that has changed
      b[0].set_data_3d( [x[i]], [y[i]], [z[i]] )
   ani = animation.FuncAnimation( fig, animate, interval=4, frames=len( t ) )
   plt.show()
#  ani.save( "demo.gif", fps=50 )


def plot_figure( qx, qy, qz, x, y, z ):
   fig = plt.figure(figsize=(8,8))
   ax = fig.add_subplot(111, projection='3d')
   ax.plot( qx, qy, qz, 'k' )                                                  # disk
   ax.plot( x, y, z, 'b' )                                                     # bob trajectory
   ax.plot( ( qx[-1], x[-1] ), ( qy[-1], y[-1] ), ( qz[-1], z[-1] ), 'g' )     # final line
   ax.plot( x[-1], y[-1], z[-1], 'ro' )                                        # final bob
   ax.view_init(30, 45)
   ax.set_xlim( -1.0, 1.0 );   ax.set_ylim( -1.0, 1.0 );   ax.set_aspect('equal')
   plt.show()



def deriv( t, Y, R, omega, L ):
    theta, phi, thetaprime, phiprime = Y
    ct, st = math.cos( theta ), math.sin( theta )
    phase = omega * t - phi
    thetaprime2 = ( phiprime ** 2 * L * st * ct + omega ** 2 * R * ct * math.cos( phase ) - g * st ) / L
    phiprime2   = ( -2 * thetaprime * phiprime * L * ct + omega ** 2 * R * math.sin( phase ) ) / ( L * st )
    return [ thetaprime, phiprime, thetaprime2, phiprime2 ]



R = 0.5
omega = 2.0
L = 1.0

Y0 = [ 0.01, 0.01, 0.0, 0.0 ]
period = 2 * np.pi / omega
tmax = 5 * period
solution = solve_ivp( deriv, [0, tmax], Y0, args=(R,omega,L), rtol=1.0e-6, dense_output=True )

t = np.linspace( 0, tmax, 1000 )
Y = solution.sol( t )
theta = Y[0,:]
phi   = Y[1,:]

# Position on disk
qx = R * np.cos( omega * t )
qy = R * np.sin( omega * t )
qz = np.zeros_like( qx )

# Trajectory
x = qx + L * np.sin( theta ) * np.cos( phi )
y = qy + L * np.sin( theta ) * np.sin( phi )
z =    - L * np.cos( theta )

#plot_figure( qx, qy, qz, x, y, z )
plot_animation( t, qx, qy, qz, x, y, z )

Trajectory:

Lagrangian Equations (ii) Cartesian Coordinates

The difficulty with Cartesian coordinates is that you are using 3 coordinates to solve a 2-degree-of-freedom problem, and so the constraint (length of pendulum being fixed) must be dealt with through a Lagrange multiplier λ.

The following vastly simplifies the solution given in Stack Exchange. The Lagrangian is

where x is the location of the bob and "q"=(R cos⁡ωt,R sin⁡ωt,0). The three Lagrangian equations, rewritten in vector form give

λ must be found from the constraint equation

Differentiating twice, and noting that

together with the expression for the acceleration above, gives the expression for 2λ below. Thus, your equation set is

This is far simpler than the expression given in Stack Exchange, but probably equivalent. It is shown in code below. The only real difference from the first code is the routine deriv(), since you now have different dependent variables.

The other problem with a constrained Lagrangian is that the initial conditions (position and velocity) must satisfy both the constraint equation AND its derivative. i.e. the string length is L and it is not initially lengthening.

Code with Cartesian dependent variables.

import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D

g = 9.81

def plot_animation( t, qx, qy, qz, x, y, z ):
   plotInterval = 1
   fig = plt.figure( figsize=(4,4) )
   ax = fig.add_subplot(111, projection='3d' )
   ax.view_init(30, 45)
   ax.set_xlim( -1.0, 1.0 );   ax.set_ylim( -1.0, 1.0 );   ax.set_aspect('equal')

   ax.plot( qx, qy, qz, 'k' )                                            # ring
   a = ax.plot( [qx[0],x[0]], [qy[0],y[0]], [qz[0],z[0]], 'g' )          # pendulum string
   b = ax.plot( [x[0]], [y[0]], [z[0]], 'ro' )                           # pendulum bob
   def animate( i ):
      a[0].set_data_3d( [qx[i],x[i]], [qy[i],y[i]], [qz[i],z[i]] )       # update anything that has changed
      b[0].set_data_3d( [x[i]], [y[i]], [z[i]] )
   ani = animation.FuncAnimation( fig, animate, interval=4, frames=len( t ) )
   plt.show()
#  ani.save( "demo.gif", fps=50 )


def plot_figure( qx, qy, qz, x, y, z ):
   fig = plt.figure(figsize=(8,8))
   ax = fig.add_subplot(111, projection='3d')
   ax.plot( qx, qy, qz, 'k' )                                                  # disk
   ax.plot( x, y, z, 'b' )                                                     # bob trajectory
   ax.plot( ( qx[-1], x[-1] ), ( qy[-1], y[-1] ), ( qz[-1], z[-1] ), 'g' )     # final line
   ax.plot( x[-1], y[-1], z[-1], 'ro' )                                        # final bob
   ax.view_init(30, 45)
   ax.set_xlim( -1.0, 1.0 );   ax.set_ylim( -1.0, 1.0 );   ax.set_aspect('equal')
   plt.show()



def deriv( t, Y, R, omega, L ):
    x, y, z, xdot, ydot, zdot = Y
    qx, qy = R * math.cos( omega * t ), R * math.sin( omega * t )
    qxdot, qydot = -omega * R * math.sin( omega * t ), omega * R * math.cos( omega * t )
    twoLambda = ( g * z + ( R ** 2 - qx * x - qy * y ) * omega ** 2 - ( xdot - qxdot ) ** 2 - ( ydot -qydot ) ** 2 - zdot ** 2 ) / L ** 2
    xddot = twoLambda * ( x - qx )
    yddot = twoLambda * ( y - qy )
    zddot = twoLambda * z - g
    return [ xdot, ydot, zdot, xddot, yddot, zddot ]



R = 0.5
omega = 2.0
L = 1.0

Y0 = [ R, 0.0, -L, 0.0, 0.0, 0.0 ]
period = 2 * np.pi / omega
tmax = 5 * period
solution = solve_ivp( deriv, [0, tmax], Y0, args=(R,omega,L), rtol=1.0e-6, dense_output=True )

t = np.linspace( 0, tmax, 1000 )
Y = solution.sol( t )
x = Y[0,:]
y = Y[1,:]
z = Y[2,:]

# Position on disk
qx = R * np.cos( omega * t )
qy = R * np.sin( omega * t )
qz = np.zeros_like( qx )

#plot_figure( qx, qy, qz, x, y, z )
plot_animation( t, qx, qy, qz, x, y, z )

Newtonian Mechanics

Your more recent post (on elastic pendulums) has prompted me to ask how the problem would be solved by Newtonian mechanics, F=ma, and not with Lagrangians.

The pendulum bob experiences two forces: an (unknown) tension T along the vector from bob to tether point and the weight of the bob, mg. Dividing by mass we get acceleration

We need to find the tension T. We do this by applying the constraint that the pendulum has a fixed length:

Differentiating twice and dividing by 2 we get

Multiplying out, using the expression for acceleration above and noting that

gives

whence we get an expression for the tension T:

Finally, substituting this back in the equation for the acceleration (and flipping signs to multiply by x-q rather than q-x) gives after a bit of rearrangement:

which is the same as the equation derived by Lagrangian mechanics above.