# Stiffness problem which abruptly stops the evaluation

We have a Lagrangian and we found the EOM. We need to solve the EOM and find the Phase-space plot. However, we are facing a stiffness problem which abruptly stops the evaluation (NDSolve::ndsz). The StiffnessSwitching method seems to not work. We have tried {"DiscontinuityProcessing" -> False, "TimeIntegration" -> "Extrapolation"} too. They are futile too.

How to resolve the issue??

here is our sample code.

lagrangian =
1/2 (1 + r[t]^2/
l^2 - (2 Sqrt[2] (1 + 1/l^2) r[t]^2)/(1 + r[t]^2)^(3/2) K^2 +
1/(1 + r[t]^2/l^2 - (2 Sqrt[2] (1 + 1/l^2) r[t]^2)/(1 + r[t]^2)^(
3/2)) r'[t]^2 + r[t]^2 L^2);

eulerLagrange[lagrangian_, vars_, dvars_] :=
Thread[(Table[D[D[lagrangian, dvar], t], {dvar, dvars}] -
Table[D[lagrangian, var], {var, vars}]) ==
ConstantArray[0, Length@vars]];

equationsOfMotion =
eulerLagrange[lagrangian, {r[t]}, {r'[t]}] /. l -> 1 // Simplify

tmax = 1000; K = 1; L = 20;
sol = NDSolve[{equationsOfMotion, r[0] == 1.5, r'[0] == 1.5}, {r[t],
r'[t]}, {t, 0, tmax}]

ParametricPlot[Evaluate[{r[t], r'[t]}] /. sol, {t, 0, 10},
AspectRatio -> 1, PlotPoints -> 200]

• It's probably not stiffness since r'[t] $\sim$ r[t]^2 near the end of integration. That is, $r(t) \sim K/(t-t_1)$, which would be a singularity instead of stiffness. A singularity is usually a feature, not a problem. If you think you have a reason that there should be no singularities, then the problem might be a mistake in modeling (setting up the DEs). Note the singularity persists if we use WorkingPrecision -> 32. That means it is less likely to be a numerical issue. Commented Aug 31, 2023 at 22:00