Simulation of a pendulum on a spinning Disk

Question

I don't know if this is the right place to ask this Question, but I have previously asked a similar question where i asked how to write a simulation on this phenomenon. I got a great answer with a Mathematica code. I started to play around a little with the values and wrote a similar program myself in python. You can find this question here: Pendulum hanging on a spinning Disk

But as I was plugging in different values I started to notice, that the constraint for the pendulum seemed to sometimes not work. The pendulum got to points, that were further, from the radius than the distance l. Also the Z axis behaved very strange often. I often had distances to the rotating disk, that seemed like the ball was jumping. I have some picutres here, so you know, what I am talking about:

For the first one:

r = 100
omega = 2
q0 =  6000
l =6000
g = 9810
tmax = 500

And the second one:

r = 0.001
omega = 0.2  
q0 =  6
l =6
g = 9.81
tmax = 500

In both cases the Result is not satisfying. I also plotted the Z-Axis just to see how the values changes. Here is the Z-Axis to the second picture:

But we can see in the question linked above, that the ODE can produce a very plausible picutre of the motion of such a pendulum.

Now the Question is why is this behaviour? It seems to be something that goes off in the integration, where we try to solve the ODE numerically. I have already plugged in the python method='DOP853' which uses the Runge-Kutta eighth-order. But still these strange simulations are the result. Is this maybe the result Overflow problem, because the numbers are to small. What do you think?

How can I encounter this problem?

I provide my Python code here, so you can try it yourself:

    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 = (g*q0*r*np.cos(omega*t) - g*q0*x - g*r*z*np.cos(omega*t) + g*x*z + omega**2*r**2*x*np.cos(omega*t)**2 + omega**2*r**2*y*np.sin(2*omega*t)/2 - omega**2*r*x**2*np.cos(omega*t) - omega**2*r*x*y*np.sin(omega*t) + omega*r**2*xp*np.sin(2*omega*t) - 2.0*omega*r**2*yp*np.cos(omega*t)**2 - 2.0*omega*r*x*xp*np.sin(omega*t) + 2.0*omega*r*x*yp*np.cos(omega*t) + r*xp**2*np.cos(omega*t) + r*yp**2*np.cos(omega*t) + r*zp**2*np.cos(omega*t) - x*xp**2 - x*yp**2 - x*zp**2)/(1.0*q0**2 - 2.0*q0*z + 1.0*r**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + 1.0*x**2 + 1.0*y**2 + 1.0*z**2)
    dypdt = (g*q0*r*np.sin(omega*t) - g*q0*y - g*r*z*np.sin(omega*t) + g*y*z + omega**2*r**2*x*np.sin(2*omega*t)/2 + omega**2*r**2*y*np.sin(omega*t)**2 - omega**2*r*x*y*np.cos(omega*t) - omega**2*r*y**2*np.sin(omega*t) + 2.0*omega*r**2*xp*np.sin(omega*t)**2 - 1.0*omega*r**2*yp*np.sin(2*omega*t) - 2.0*omega*r*xp*y*np.sin(omega*t) + 2.0*omega*r*y*yp*np.cos(omega*t) + r*xp**2*np.sin(omega*t) + r*yp**2*np.sin(omega*t) + r*zp**2*np.sin(omega*t) - xp**2*y - y*yp**2 - y*zp**2)/(1.0*q0**2 - 2.0*q0*z + 1.0*r**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + 1.0*x**2 + 1.0*y**2 + 1.0*z**2)
    dzpdt =(-1.0*g*r**2 + 2.0*g*r*x*np.cos(omega*t) + 2.0*g*r*y*np.sin(omega*t) - 1.0*g*x**2 - 1.0*g*y**2 + 1.0*omega**2*q0*r*x*np.cos(omega*t) + 1.0*omega**2*q0*r*y*np.sin(omega*t) - 1.0*omega**2*r*x*z*np.cos(omega*t) - 1.0*omega**2*r*y*z*np.sin(omega*t) + 2.0*omega*q0*r*xp*np.sin(omega*t) - 2.0*omega*q0*r*yp*np.cos(omega*t) - 2.0*omega*r*xp*z*np.sin(omega*t) + 2.0*omega*r*yp*z*np.cos(omega*t) + 1.0*q0*xp**2 + 1.0*q0*yp**2 + 1.0*q0*zp**2 - 1.0*xp**2*z - 1.0*yp**2*z - 1.0*z*zp**2)/(1.0*q0**2 - 2.0*q0*z + 1.0*r**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + 1.0*x**2 + 1.0*y**2 + 1.0*z**2)
    return [dxdt, dydt, dzdt, dxpdt, dypdt, dzpdt]

r = 0.001
omega = 0.2  
q0 =  6
l =6
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 = 500

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)

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(Qx, Qy, 0, color='purple')
ax.view_init(30, 45)
plt.show()

plt.plot(t_values, p_values[2], color='blue')
plt.grid(True)
plt.show()

Answer 1

Most numerical integration schemes don't exactly preserve quantities that are conserved by the true ODE. It is often the case that conserved quantities will slightly drift or oscillate over time. There are many methods for fixing this such as using more accurate integration, explicitly constraining the numerical scheme, using symplectic integration, etc. An overview of these methods can be found in this excellent answer by Chris Rackaukas. He is one of the authors of the Julia ODE package he uses in this post, DifferentialEquations.jl, which has Python bindings, found here.