Problem statementI am trying to simulate a spherical pendulum, withrod 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 descriptionThe 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 triedUsually 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 constantvelocity is perpendicular to the rodas 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 behaviorThat 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 forMy 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 circumstancesit has to still produce meaningful results after 15 days of continuous operationTherefore, 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 exampleA (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 functionfunction 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];    endendfunction 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
Numerical instability of spherical pendulum
stability;time integration
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 andthat 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.