Simulation of a spring loaded pendulum on a spinning disk

I want to write a simulation in Python similar to the one described in Simulation of a Pendulum hanging on a spinning Disk.

But I want the system to be spring loaded. So instead of the mass hanging from a thread, I want the mass to hang from a spring, that rotates. I have tried putting the ODE together but it takes forever to calculate:

import sympy as sp
from IPython.display import display

R = sp.symbols('R')
omega = sp.symbols('omega')
t = sp.symbols('t')
phi = sp.Function('phi')(t)
theta = sp.Function('theta')(t)
s = sp.Function('s')(t)
L = sp.symbols('L')
m = sp.symbols('m')
k = sp.symbols('k')
g = sp.symbols('g')

x = R*sp.cos(omega*t)+(L+s)*(sp.sin(theta)*sp.cos(phi))
y = R*sp.sin(omega*t)+(L+s)*(sp.sin(theta)*sp.sin(phi))
z = -(L+s)*sp.cos(theta)

xs = sp.diff(x,t)
ys = sp.diff(y,t)
zs = sp.diff(z,t)

v = xs**2 + ys**2 + zs**2
Ekin = 0.5*m*v
Epot = g*(L+s)*sp.cos(theta)+0.5*k*s**2

L = Ekin + Epot

#display(L)

ELTheta = sp.diff(sp.diff(L,sp.Derivative(theta,t)), t) + sp.diff(L,theta)
ELPhi = sp.diff(sp.diff(L,sp.Derivative(phi,t)), t) + sp.diff(L,phi)
ELs = sp.diff(sp.diff(L,sp.Derivative(s,t)), t) + sp.diff(L,s)

Eq1 = sp.Eq(ELTheta,0)
Eq2 = sp.Eq(ELPhi,0)
Eq3 = sp.Eq(ELs,0)

LGS = sp.solve((Eq1,Eq2,Eq3),(sp.Derivative(theta,t,2),sp.Derivative(phi,t,2),sp.Derivative(s,t,2)))

thetadd = sp.simplify(LGS[sp.Derivative(theta,t,2)])
phidd = sp.simplify(LGS[sp.Derivative(phi,t,2)])
sdd = sp.simplify(LGS[sp.Derivative(s,t,2)])

I don't know whether I chose the right original condition. Is there a simpler and faster way to compute this or a different formula to the problem, that would simplify it?

Once you remove the constraint of a fixed-length pendulum it is probably easier to solve this in Cartesian coordinates using Newton's 2nd Law (F=ma), rather than via the Lagrangian equations of motion.

Consider a spring of natural length L and stiffness k, one end attached to a point on the edge of a disk of radius R rotating with angular velocity ω, the other holding a bob of mass m.

The tether point on the disk has time-dependent coordinates

The bob has coordinates

The bob is subject to two forces: the elastic force from the spring of magnitude

and its weight

The equation of motion is just

or

where n is a unit vector toward the tether and g is the acceleration due to gravity.

Note that there is no energy-dissipating mechanism, so make sure that you avoid resonance (i.e. the natural circular frequency of the spring, sqrt(k/m), should not be close to the circular frequency, omega, of the rotating point).

For a stable, rather than very "bouncy", orbit, see below the main code.

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                                                             # acceleration due to gravity (m/s^2)
gravity = np.array( [ 0.0, 0.0, -g ] )

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( -2.0, 2.0 );   ax.set_ylim( -2.0, 2.0 );   ax.set_aspect('equal')

   ax.plot( qx, qy, qz, 'k' )                                        # tether point
   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 ):                                                 # update anything that has changed
      a[0].set_data_3d( [qx[i],x[i]], [qy[i],y[i]], [qz[i],z[i]] )   
      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 deriv( t, Y, R, omega, L, k, m ):
    x, v = Y[0:3], Y[3:6]                                                      # bob position and velocity (vectors)
    q = np.array( [ R*math.cos(omega*t), R*math.sin(omega*t), 0.0 ] )          # position vector of tether
    spring = q - x                                                             # vector from bob to tether
    length = np.linalg.norm( spring )                                          # stretched length of spring                                             
    unitVector = spring / length                                               # direction of elastic force
    extension = length - L                                                     # extension of spring
    a = ( k / m ) * extension * unitVector + gravity                           # F combines elastic force and gravity
    return np.hstack( ( v, a ) )                                               # return velocity and acceleration


R = 0.5            # radius of disk (m)
omega = 2.0        # angular velocity of disk (rad/s)
L = 1.0            # natural length of spring (m)
k = 20.0           # stiffness of spring (N/m)             # warning: avoid resonance!
m = 2.0            # mass (kg)
Y0 = np.array( [ R, 0.0, -L, 0.0, 0.0, 0.0 ] )                                 # initial position and velocity

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

t = np.linspace( 0, tmax, 1000 )
Y = solution.sol( t )

# Position of bob
x = Y[0,:]
y = Y[1,:]
z = Y[2,:]

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

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

It is possible to give the pendulum bob an initial position and velocity that will put it in a stable orbit (steady motion around a circle, in synchronisation with the tether on the disk). If you set the centripetal force equal to the horizontal component of the elastic force and the weight equal to the vertical component of the elastic force then you can solve a couple of simultaneous equations for the tangent of the angle to the vertical and the spring extension and thereby provide the relevant initial conditions for Y0. Try replacing the initialisation of Y0 with the following lines. There are also probably synchronous solutions with the bob lagging the tether point, but I haven't investigated them.

# "Equilibrium motion" - moves in a circle, synchronised with the tether
t, e = 0, 0                  # tan(theta) and extension
w2g = omega ** 2 / g
mgk = m * g / k
for _ in range( 100 ):
     c = 1 / math.sqrt( 1 + t*t );   s = t * c        # cosine and sine
     t = w2g * ( R + (L+e) * s )
     e = mgk / c
     # print( t, e )
x = R + (L+e) * s
z =   - (L+e) * c
Y0 = np.array( [ x, 0.0, z, 0.0, x * omega, 0.0 ] )