Implementing a spherical pendulum

Question

I am trying to implement a spherical pendulum. The Lagrangian (which I haven't fully understood so yet) based on $l$, θ and φ taken from this page result in the equations:

\begin{align} \ddot{\theta}&=\dot{\varphi}^2\sin\theta\cos\theta-\frac{g}{l}\sin\theta\\ \ddot{\varphi}&=-2\dot{\theta}\dot{\varphi}\cot\theta \end{align}

I am treating $\ddot{\theta}$ as the acceleration and of $\dot{\theta}$ as the velocity of $\theta$. Is this correct? Now, each step of the motion is implemented as follows (Python):

theta_f = pow(phi_v, 2) * sin(theta) * cos(theta) - G / L * sin(theta)
phi_f = - 2 * theta_v * phi_v / tan(theta)

theta_v += theta_f / timesteps
phi_v += phi_f / timesteps

theta += theta_v
phi += phi_v

This works as long as phi_v ($\dot{\phi}$) is 0 or close to 0.

If $\dot{\phi}_\text{initial} > 0$, the movement is erroneous.

My initial values are

theta = 0.8
phi = 0.5
timesteps = 60
L ~ 2
G = 2.0
theta_v = 0.0
phi_v = 0.1

After a few iterations, the code produces a math range error as phi_v gets too large. I have found this question which could explain the math rounding error.

I am using using 60 samples per second, because there will be real-time interaction. Approximated values will be totally fine, but I can't believe that the current state is simply and approximation error.

How can I correct my code to simulate the spherical pendulum?

Answer 1

The problem seems to be that you are not accounting for the step size in time when integrating. This should be obvious when you are doing things like

phi_v += phi_f

in the code.¹ The acceleration cannot be simply added to the velocity (the units don't match!). The expected relationship is,

$$a_\phi =\frac{\mathrm d v_\phi}{\mathrm dt}\approx\frac{\Delta v_\phi}{\Delta t}\implies v_\phi\approx v_\phi'+ a_\phi \Delta t$$ where $v_\phi'$ is the previous value is the velocity.

What you should be doing is using the 4 equations (2 positions, 2 velocities), $$ \mathbf{v}=\frac{\mathrm d\mathbf x}{\mathrm dt} \text{ and } \mathbf{a}=\frac{1}{m}\frac{\mathrm d\mathbf f}{\mathrm dt} , $$ where $\mathbf v=[\dot{\phi},\,\dot{\theta}]^T$ and similarly for $\mathbf x$, and solve those via a leap-frog integrator, such as velocity Verlet. However, since you have a velocity-dependent acceleration (i.e., $a= f(\mathbf x,\,\mathbf v)$), you need to use a modified version of the integrator (which I talk about in this other answer of mine),

\begin{align} a_1(t)&=\text{compute from $x(t)$ and $v(t)$}\tag{1}\\ x(t+\Delta t)&=x(t)+v(t)\Delta t+\frac{1}{2}a_1(t)\Delta t^2\tag{2}\\ v(t+\Delta t)&=v(t)+a_1(t)\Delta t\tag{3}\\ a_2(t+\Delta t)&=\text{compute from $x(t+\Delta t)$, $v(t+\Delta t)$}\tag{4}\\ v(t+\Delta t)&=v(t+\Delta t)+\frac{1}{2}\left(a_2(t+\Delta t)-a_1(t)\right)\Delta t\tag{5} \end{align} which you then just loop through until satisfied (e.g., $t\geq t_\text{end}$). Provided $\Delta t$ is small enough, this should provide a more stable solver than velocity Verlet integration alone.

^{1. OP added the timesteps term to the equations in version 6; this answer was posted prior to that edit.}

Answer 2

It looks to me like you are trying to solve a system of ordinary differential equations by estimating the solution for small steps. Am I correct? If so then there are a number of things to consider.

First, you need to have a good estimate for the derivative operator to approximate the future state given the initial conditions. This type of problem is called an initial value problem (i.v.p).

Since the equation involves both first and second order time derivatives a common step is to define momentum variables and propagate a larger system. This will help prevent errors as the total time increases. I've seen that problem in ray tracing codes and what I am suggesting is so common that most people never propagate a second order equation. The idea is:

p_theta = d(theta)/dt

p_phi = d(phi)/dt

then your equation is

d(p_theta)/dt = (p_phi)^2*cos(theta)*sin(theta) - g/l*sin(theta)

d(p_phi)/dt = -2*p_theta*p_phi*cot(theta)

plus the two equations that define the p's.

A simple Euler step would be implemented as,

p_theta(t0+dt) = p_theta(t0) + (d(p_theta)/dt)(t0)*dt
p_phi(t0+dt) = p_phi(t0) + (d(p_phi)/dt)(t0)*dt

plus the equations for theta and phi in the definitions of p.

This is the typical setup. Now, the Euler step is very poor and never recommended. You would do better to implement a higher order method like RK4 or RK5(4), etc, with step size control.

Aside from that the use of angles is sometimes an issue as it leads to the possibility that at some step you cannot uniquely determine the next value due to the system being singular. In aerospace simulations they use quarterions to fix this. You can refer to Goldstein Classical Mechanics for details on the math or Zipfel Modeling and Simulation of Aerospace Vehicle Dynamics. I think you don't need this machinery right now. You need to write the equations properly and try a simple step algorithm before getting too sophisticated. I would think that Python has ODE solver package so you don't need to write your own, just set up correctly and call. I hope that helps.