# Implementing a spherical pendulum

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?

• You're talking about accelerations, not forces. There shouldn't be a mass involved anymore. – J. Murray Dec 13 '18 at 14:02
• Have you tried much larger values of initial angular velocity? It could be that the ball is almost, but not quite, hitting the center in these images. – probably_someone Dec 13 '18 at 14:16
• In your code, for $\ddot\phi$, why is there a $4$ instead of a $2$, and why is there a $\sin\theta$ term? – BioPhysicist Dec 13 '18 at 14:43
• Ok but why the $\sin\theta$ term? I would feel like that would make a huge difference if you are dividing by it. When $\theta=0$ you would have an infinite acceleration. – BioPhysicist Dec 13 '18 at 14:46
• Would Computational Science be a better home for this question? – Qmechanic Dec 13 '18 at 17:02

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.

• I upvoted. Your answer was very helpful. I will just take a lot of time to process the information and read up on the new concepts, which I was completely unfamiliar with. Thank you again for your help! – Leander Dec 28 '18 at 10:02

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.1 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.