Question: 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?

Answer:  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 ] )