Why does this non-stiff ode requires a stiff solver? This post is related to another one in physics.stackexchange. But it might be more of a quadrature than a physics problem.
The following set of ordinary differential equations describe the motion of two mass points. They are coupled by a couple of geometric constraints yielding in the following equations in $x_1(t), z_2(t)$ (see Theoretische Physik, Bartelmann et al., 1e, 2015):
$$
\begin{eqnarray}
m_1 \ddot{x}_1 & = & -\frac{\dot{x}_1^2 + \dot{z}_2^2 - z_2g}{2(x_1^2/m_1 + z_2^2/m_2)}x_1, \\
m_2 \ddot{z}_2 & = & -\frac{\dot{x}_1^2 + \dot{z}_2^2 - z_2g}{2(x_1^2/m_1 + z_2^2/m_2)}z_2 - m_2g.
\end{eqnarray}
$$
Integrating with a Runge-Kutta scheme of order 4, and initial conditions chosen consistently with the constriants to be $\mathbf{x}_0 = (x_{10}, z_{20}) = (0, -L)$, violates one of the geometric constraints $x_1^2 + z_2^2 - L^2 = 0$ for time $t > 0$.
If integrated with an implicit Runge-Kutta method of the Radau IIA family of order 5, considered a stiff solver, the error is reduced by an order of magnitude of 3.
If I consider a characteristic of the system, stating that the expression $E := m/2(\dot{x}_1^2 + \dot{z}_2^2 + (x_1 + z_2)g)$ forms an invariant/is preserved, if $m_1 = m_2 = m$, and considering the geometric constraint $x_1^2 + z_2^2 = L^2$, the same system can be expressed as
$$
\begin{eqnarray}
m\ddot{x}_1 & = & -\frac{E - (x_1/2 + z_2)g}{L^2}x_1, \\
m\ddot{z}_2 & = & -\frac{E - (x_1/2 + z_2)g}{L^2}z_2,
\end{eqnarray}
$$
i.e. a system of ordinary differential equations which can be split into a decoupled linear, and a nonlinear part.
Question:
As can be seen from the linear part of the model, the eigenfrequencies are identical. Particularly in the case of $g=0$, how come I need a stiff solver to account for numerical instabilities then?
 A: According to my movement equations, the movement behavior is completely acceptable.
Calling
$$
\cases{
p_1=(x_1,y_1)\\
p_2 = (x_2,y_2)\\
r_1 = \|p_1\|^2-l_1^2/4\\
r_2 = \|p_2-2p_1\|^2-l_2^2/4
}
$$
the lagrangian reads
$$
L(p_1,\dot p_1, p_2,\dot p_2,\mu_1,\mu_2) = \frac 12(m_1\|\dot p_1\|^2+m_2\|\dot p_2\|^2)-m_1 g y_1-m_2 g y_2+\mu_1 r_1+\mu_2 r_2
$$
so the movement equations are obtained from
$$
\cases{
\frac{d}{dt}\frac{\partial L}{(\dot p_1,\dot p_2)}-\frac{\partial L}{( p_1,p_2)}=0\\
\ddot r_1 = 0\\
\ddot r_2 = 0
}
$$
Follows a script in MATHEMATICA which derives symbolically those equations, as well as integration results for your initial conditions.
p1 = {x1[t], y1[t]};
p2 = {x2[t], y2[t]};
T = 1/2 (m1 D[p1, t]. D[p1, t] + m2  D[p2, t]. D[p2, t])
U = m1 g y1[t] + m2 g y2[t]
r1 = p1.p1 - l1^2/4
r2 = (p2 - 2 p1).(p2 - 2 p1) - l2^2/4
L = T - U + mu1 r1 + mu2 r2
vars = Join[p1, p2]
equs = D[Grad[L, D[vars, t]], t] - Grad[L, vars]
equsmov = Join[equs, {D[r1, t, t], D[r2, t, t]}]
solvars = Solve[equsmov == 0, Join[D[vars, t, t], {l1, l2}]][[1]]
parms = {m1 -> 1, m2 -> 1, g -> 1, l1 -> 1, l2 -> 1}
equs0 = D[vars, t, t] /. solvars /. parms
cinits = {x1[0] == x1'[0] == 0, y1[0] == -l1/2, y1'[0] == 0, x2[0] == l2/2, x2'[0] == 0, y2[0] == -l1, y2'[0] == 0} /. parms

tmax = 30;
solmov = NDSolve[Join[Thread[D[vars, t, t] == equs0], cinits], vars, {t, 0, tmax}]
gr1 = ParametricPlot[Evaluate[p1 /. solmov], {t, 0, tmax}];
gr2 = ParametricPlot[Evaluate[p2 /. solmov], {t, 0, tmax}];
Show[gr1, gr2, PlotRange -> All]


Plot[Evaluate[(p2 - 2 p1).(p2 - 2 p1) /. solmov], {t, 0, tmax}]


Note that the distance decays very slowly and with a hamiltonian integrator this decay can be almost eliminated,
