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];
% deThetadth*theta_dot + deThetadph*phi_dot % Should equal eTheta_dot
% 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
