# Numerical instability of spherical pendulum

## Problem statement

I am trying to simulate a spherical pendulum, with

• rod length $r$
• south-polar angle $\theta$ and azimuthal angle $\phi$
• initial values $(\theta_0,\phi_0)= (0,0)$

My particular pendulum is hanging in a construction, which can be moved around by external forces (say, an experimenter). This construction has well-defined coordinate axes, which therefore also serve as the coordinate axes for the pendulum.

Because of this:

• Pendulum movement takes place in a non-inertial reference frame, with:

• non-zero angular velocity $\mathbf{\Omega}$ and angular acceleration $\mathbf{\dot{\Omega}}$
• non-zero linear velocity $\mathbf{v}_{\text{rel}}$ and linear acceleration $\mathbf{a}_{\text{rel}}$
• Because of the nature of my problem, the gravitational acceleration $\mathbf{g}$ is always pointing in the negative $z$-axis of the pendulum frame.

I am experiencing severe numerical instabilities when integrating the equations of motion, when I include the non-inertial forces and my initial values are such that $|\theta| \lesssim 1^\circ$. The type of the numerical integrator just seems to have an effect on the moment of divergence, but divergence is always there.

Strangely enough, I can integrate everything perfectly fine for any set of initial values, when I set $\mathbf{\Omega}=\dot{\mathbf{\Omega}}=\mathbf{v}_{\text{rel}}=\mathbf{a}_{\text{rel}}=\mathbf{0}$.

## Mathematical description

The equations of motion for a "normal" spherical pendulum are

\begin{align} \ddot{\theta} &= \left(\dot{\phi}^2\cos\theta - g/r \right)\sin\theta\\ \ddot{\phi} &= \frac{-2\dot\theta\dot\phi\cos\theta}{\sin\theta} \end{align}

If I manually set $\ddot\phi=0$ when non-finite values are encountered (precaution for division-by-zero when $\theta=0$), I have no problem whatsoever to numerically integrate these equations.

The total force on this particular pendulum is an addition of the pendulum weight and non-inertial forces:

$$\mathbf{F}_{\text{tot}} = m\left( \mathbf{g} - \dot{\mathbf{\Omega}}\times\mathbf{r} - \mathbf{\Omega}\times\left(\mathbf{\Omega}\times\mathbf{r} \right) - 2\mathbf{\Omega}\times\mathbf{v}\right)$$

which are the gravitational force, Euler force, centripetal force and Coriolis force, respectively.

These cross products are much easier to do in Cartesian coordinates. Therefore, in my software implementation of the differential equation, I convert $r$, $\theta$ and $\phi$ to Cartesian coordinates $x$, $y$, $z$, carry out the cross products, add the forces to the final acceleration, and convert back to the angular accelerations in spherical coordinates $\ddot\theta$ and $\ddot\phi$.

Equations to convert to spherical coordinates (using $\dot{r} = \ddot{r} = 0$):

\begin{align} \mathbf{r} &= r\mathbf{e}_r \\ \mathbf{v} &= r\dot{\phi}\sin\theta \mathbf{e}_\phi + r\dot{\theta}\mathbf{e}_\theta \end{align}

where the unit vectors of the spherical coordinates have been used:

\begin{align} \mathbf{e}_{r} &= \mathbf{r}/|\mathbf{r}| \\ \mathbf{e}_{\theta} &= \left[ \begin{matrix}\cos\theta\cos\phi\\ \cos\theta\sin\phi\\ \sin\theta\end{matrix} \right]\\ \mathbf{e}_{\phi} &= \left[ \begin{matrix}-\sin\phi\\ \cos\phi\\ 0\end{matrix} \right] \end{align}

as well as their time derivatives:

\begin{align} \dot{\mathbf{e}_{r}} &= \dot{\theta}\mathbf{e}_{\theta} + \dot\phi\sin\theta\mathbf{e}_{\phi} \\ \dot{\mathbf{e}_{\theta}} &= -\dot\theta\mathbf{e}_{r} + \dot\phi\cos\theta\mathbf{e}_{\phi} \\ \dot{\mathbf{e}_{\phi}} &= -\dot\phi\sin\theta\mathbf{e}_{r} - \dot\phi\cos\theta\mathbf{e}_{\theta} \end{align}

Deriving the velocity $\mathbf{v}$ once more gives:

\begin{align} \mathbf{a} = &\left(-r\left(\dot\theta^2 + \dot\phi^2\sin^2\theta\right) \right) \mathbf{e}_{r} + \\ &\left( r\left(\ddot\theta - \dot\phi^2\sin\theta\cos\theta\right) \right) \mathbf{e}_{\theta} + \\ &\left( r\left(2\dot\phi\dot\theta\cos\theta + \ddot\phi\sin\theta\right) \right) \mathbf{e}_{\phi} \end{align}

so that

\begin{align} \ddot\theta &= \frac{\mathbf{a_{\text{tot}}}\cdot\mathbf{e}_{\theta} + r\dot\phi^2\sin\theta\cos\theta}{r}\\ \ddot\phi &= \frac{\mathbf{a_{\text{tot}}}\cdot\mathbf{e}_{\phi} - 2r\dot\phi\dot\theta\cos\theta }{r\sin\theta} \end{align}

where $\mathbf{a_{\text{tot}}}$ is the acceleration of the pendulum mass, including gravity and all the non-inertial forces.

I am pretty sure these equations are correct, as I have made many assertions, plots, etc. for each and every step of the way. I can also simulate a "regular" spherical pendulum, as well as one in a rotating/accelerating frame, and all the outcomes are very realistic.

Unfortunately, all just for initial values well outside the south-polar region :)

For some strange reason, it is numerically OK to have extremely large values for $\ddot\phi$ when there is only gravity in $\mathbf{a}_{\text{tot}}$, but not if there are additional forces in there...

## What I have tried

Usually I try to avoid spherical coordinates like the plague for precisely this sort of reason, but I fear in the case of a spherical pendulum there is no real way around it.

As the problem is not physical (there is nothing "physically special" about the direct-down position of the pendulums) but just due to the choice of the coordinate system, I thought I'd try to dump the spherical coordinates and integrate everything directly in Cartesian coordinates. This resulted in a much simpler formulation, and much more elegant code, but there was another problem: I could not keep the length of the rod constant :)

When integrating the pendulum like this, you have to include these constraint equations:

• length of the rod is constant
• velocity is perpendicular to the rod

as well as calculate additional terms, like the centripetal acceleration of the pendulum mass towards the hinge point.

I know there are ways to do that, but these usually place demands on the type of integrator used to solve the problem. In my final application, I can only use one of a small subset of pre-defined integrators, so keeping the constraints satisfied is virtually impossible. For any of the integrators I can work with, roundoff error quickly accumulates and starts violating the length constraints, leading to divergence.

Another thing I tried (although I didn't like it much) is to define two "zones of avoidance" at the poles. That is, if the angle $\theta$ is within a small region around the pole, the value for $\ddot\phi$ is defined to be zero. The problem with this approach is:

• how big should the region be in order to guarantee stability?
• it leads to non-physical behavior

That is: while it does seem to be more stable this way, there are "strange discontinuities" visible in the solutions. Moreover, the results completely disagree with results obtained without this "zone of avoidance"...

## What I am looking for

My final application will need to be space-graded. That means a lot of things, but the most important ones are:

• it has to produce meaningful results under all circumstances
• it has to still produce meaningful results after 15 days of continuous operation

Therefore, I am looking for ways to "get rid" of this singularity at the poles somehow, preferably preserving physical correctness :)

I am also looking for a good explanation of why this occurs, as I have been staring at it for three days now and I am totally blinded by my own derivations :)

## Minimum working example

A (not-so-)minimum working example in MATLAB:

function MWE

%% Pendulum parameters

clc

g = [0; 0; -9.8];
L = 2;

tspan = [0 50];

%% Plot Omega forcing term

tt = 0:0.1:25;
Omega = zeros(numel(tt),3);
for ii = 1:numel(tt)
Omega(ii,:) = Omega_fcn(tt(ii)); end
plot(...
tt, Omega(:,1), 'r',...
tt, Omega(:,2), 'g',...
tt, Omega(:,3), 'b');
xlabel('t'), ylabel('Angular rate')
t = title('$\mathbf{\Omega}(t)$');
set(t, 'interpreter', 'LaTeX');
l = legend('$\mathbf{\Omega}_x$', '$\mathbf{\Omega}_y$', '$\mathbf{\Omega}_z$');
set(l, 'interpreter', 'LaTeX');

%% Normal spherical pendulum, "good" initial values

Omega     = zeros(3,1);
Omega_dot = zeros(3,1);

y0 = [pi/4 0 0 0.4];
f  = @(t,y)equationsOfMotion(y, L,g,Omega,Omega_dot);

[~, y] = ode113(f, tspan, y0);

theta = y(:,1);
phi   = y(:,2);

X = +L*cos(phi).*sin(theta);
Y = +L*sin(phi).*sin(theta);
Z = -L*cos(theta);

figure(2)
plot3(X,Y,Z, 'r')
axis tight equal
t = title('$\mathbf{\Omega} = \dot{\mathbf{\Omega}} = \mathbf{0}$, "good" initial values');
set(t, 'interpreter', 'LaTeX');

%% Spherical pendulum in rotating frame, "good" initial values

Omega_dot = zeros(3,1);

y0 = [pi/4 0 0 0.4];
f  = @(t,y)equationsOfMotion(y, L,g,Omega_fcn(t),Omega_dot);

[~, y] = ode113(f, tspan, y0);

theta = y(:,1);
phi   = y(:,2);

X = +L*cos(phi).*sin(theta);
Y = +L*sin(phi).*sin(theta);
Z = -L*cos(theta);

figure(3)
plot3(X,Y,Z, 'r')
axis tight equal
t = title('$\mathbf{\Omega} = \mathbf{\Omega}(t),\ \dot{\mathbf{\Omega}} = \mathbf{0}$, "good" initial values');
set(t, 'interpreter', 'LaTeX');

%% Spherical pendulum in rotating frame, "BAD" initial values

Omega_dot = zeros(3,1);

y0 = [0 0 0 0];
f  = @(t,y)equationsOfMotion(y, L,g,Omega_fcn(t),Omega_dot);

% FAILURE!
[~, y] = ode113(f, tspan, y0);
% FAILURE!

theta = y(:,1);
phi   = y(:,2);

X = L*cos(phi).*sin(theta);
Y = L*sin(phi).*sin(theta);
Z = -L*cos(theta);

figure(4)
plot3(X,Y,Z, 'r')
axis tight equal
t = title('$\mathbf{\Omega} = \mathbf{\Omega}(t),\ \dot{\mathbf{\Omega}} = \mathbf{0}$, "BAD" initial values');
set(t, 'interpreter', 'LaTeX');

end

% Forcing function
function omega = Omega_fcn(t)

omega = [0; 0; 0];

t1Start = 1;
t2Start = 10;
Duration  = 1;

if (t >= t1Start) && (t < t1Start + Duration)
omega = [
+0.0175
-3.2e-4
-3.35e-4]*(t-t1Start);
end

if (t >= t1Start + Duration) && (t < t2Start)
omega = [
+0.0175
-3.20e-4 + 7.500e-006*(t-t1Start-Duration)
-3.35e-4 - 8.125e-006*(t-t1Start-Duration) ];
end

if (t >= t2Start) && (t < t2Start + Duration)
omega= [...
+0.0175 - 0.01750*(t-t2Start)
-2.6e-4 + 3.10e-4*(t-t2Start)
-4e-4   + 3.75e-4*(t-t2Start) ];
end

if (t >= t2Start + Duration)
omega= [...
0
0.5e-004
-0.25e-004];
end

end

function dy = equationsOfMotion(y, L, g, Omega, Omega_dot)

% Extract spherical angles and angular rates
theta = y(1);    theta_dot = y(3);
phi   = y(2);    phi_dot   = y(4);

% pre-compute sines/cosines
sT = sin(theta);   sP = sin(phi);
cT = cos(theta);   cP = cos(phi);

% Unit vectors of the local spherical coordinates
eTheta = [cT.*cP; cT.*sP; sT];
ePhi   = [-sP;  cP;  0];
eR     = [sT*cP;  sT*sP;  -cT];
% ASSERT:
% eR = r/L;
% |eR| = |ePhi| = |eTheta| = 1
% cross(eR,     ePhi)   = eTheta
% cross(ePhi,   eTheta) = eR       (orthonormal basis)
% cross(eTheta, eR)     = ePhi

% Pendulum vector: vector to bob mass, with hinge point as origin
r = L*eR;

% Time derivatives of unit vectors
% eR_dot     =       0         +  theta_dot *eTheta  +  phi_dot*sT*ePhi;
% eTheta_dot = -theta_dot *eR  +         0           +  phi_dot*cT*ePhi;
% ePhi_dot   = -phi_dot*sT*eR  -  phi_dot*cT*eTheta  +         0       ;
%
% ASSERT:
%
% deRdth = [cT*cP; cT*sP; sT];
% deRdph = [-sT*sP; sT*cP; 0];
% deRdth*theta_dot + deRdph*phi_dot % Should equal eR_dot
%
% dePhidph = [-cP; -sP; 0];
% dePhidph*phi_dot % Should equal ePhi_dot
%
% deThetadth = [-sT*cP;  -sT*sP;  cT];
% deThetadph = [-cT*sP;  cT*cP;  0];

% Velocities
v = L*phi_dot*sT*ePhi + L*theta_dot*eTheta;
% ASSERT: v must equal:
%  v = R*eR_dot

% Pendulum weight
arel = g;

% d(Omega)/dt × r
% The Euler force
Euler = cross(Omega_dot, r);

% 2·(Omega × v)
% The Coriolis force
Coriolis = 2*cross(Omega, v);

% Omega × (Omega × r)
% The local centripetal force
Centripetal = cross(Omega, cross(Omega, r));

% Ficticious forces due to non-intertialness of the frame:
a = arel - Centripetal - Coriolis - Euler;

% All components parallel to the rod are absorbed by the tension/compression of
% the rod. Moreover, they are irrelevant for theta_dot_dot and phi_dot_dot.
% Thus, just to be sure, subtract those components:
a = a - dot(a,eR)*eR;

% Now convert the acceleration to spherical coordinates
theta_dot_dot = ( dot(a,eTheta) + L*phi_dot^2*sT*cT )/L;
phi_dot_dot   = ( dot(a,ePhi)   - 2*L*phi_dot*theta_dot*cT ) / (L*sT);

% division by zero
if ~isfinite(phi_dot_dot)
phi_dot_dot = 0; end

% Thus, the derivative of the state vector:
dy = [
theta_dot
phi_dot
theta_dot_dot
phi_dot_dot
];

end


Interesting problem. But this will be only a hint to a possible solution.

I am not so much in the details, but I assume, you can write your system in the coordinates of the endpoint of the pendulum $p=[p_1,p_2,p_3]$ and velocities $v=[v_1,v_2,v_3]$ in cartesian coordinates subject to the system of equations

\begin{align} \dot p &= v + f_1 \\ \dot v + g_p(p)^T\lambda &= f_2, \\ g(p) := p_1^2 + p_2^2 + p_3^2 - r^2 &= 0 \end{align}

where $f_*$ contain all external forces, $\lambda$ is a multiplier, and $r$ is the length of the rod.

This is a DAE and reduces to a multibody system in an inertial coordinate system. Standard integration will inevitably lead to deviations from the constraint.

For multibody system there is the approach of minimal extension, that reduces the index of the DAE. From a practical point of view, this means finding an equivalent reformulation of the system, for which standard integrators, e.g. BDF methods work.

In the paper by Kunkel&Mehrmann, in which minimal extension was introduced, there is a concrete example of how to simulate the movement induced by gravity of a masspoint inside a paraboloid (Eqn. 62).

I think this is quite similar to what you are trying, except from your additional difficulties that come with the

• accelerated frame and
• that the paraboloid is easier to parametrize than the sphere .

The authors of the mentioned paper still use a specialized DAE solver. However, if you are fine with lower accuracy - Backward Euler should work well for the reformulated system, as we have shown for semi-discretized Euler equations in a recent preprint.

Disclaimer: Volker Mehrmann is the supervisor of my PhD thesis.

• Thanks! I will certainly try this. To be honest though, it seems a bit complicated for this type of problem, especially since it is already solvable everywhere except at two small patches... – Rody Oldenhuis Sep 25 '13 at 14:08
• Yes, since the approach is formulated for a general systems, there is a lot of abstract considerations and computations involved. If you have the equations ready (in the multibody formulation), I can do the derivation of the minimal extended system for you. – Jan Sep 25 '13 at 14:32

Have you tried using a pseudo-Cartesian (not sure of the proper terminology) coordinate system like $u = \theta \cos \phi$, $v = \theta \sin \phi$? This will avoid the instability at the south pole but is not continuous at the north pole - I don't know whether that is OK for your purposes.

• Thanks! Yes, I thought of something similar to this, and I suppose for the current version it will indeed be sufficient. Thing is -- I won't be surprised if after a few weeks my customer will come back to me complaining that he cannot swing his pendulum all the way around without blowing up the simulation :) Your solution essentially shifts the problem from both equilibrium positions to just the unstable one (North pole), so it is a significant improvement, but only a partial solution. – Rody Oldenhuis Oct 14 '13 at 5:08
• Moreover, the North pole is actually the stable equilibrium position in some of the scenarios I need to simulate. That's why I ended up re-defining the spherical coordinate system, such that the poles are in the XY plane at those locations where the probability that the pendulum would ever come there is minimized (for my specific application, of course). Also not very elegant, and also not a complete solution, but sufficient for my case. Thanks again! – Rody Oldenhuis Oct 14 '13 at 5:10