1
$\begingroup$

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]
$\endgroup$
1
  • 2
    $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Aug 31, 2023 at 22:00

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.