Strange behavior of Method ->{"EquationSimplification" -> "Residual"}

Question

I'm experiencing some strange behavior when I try to use ParametricNDSolve and Method -> {"EquationSimplification" -> "Residual"}. (This is a follow up to my previous question.)

Essentially, I'm observing wild behavior for finely tuned values of the parameter L.

f[r_] := 1 - 1/r;

soln = ParametricNDSolve[{r'[\[Lambda]]^2 == 
    1 - L^2/r[\[Lambda]]^2*f[r[\[Lambda]]], r[0] == 1000}, 
  r, {\[Lambda], 0, 2000}, {L}, 
  Method -> {"EquationSimplification" -> "Residual"}]

Plot[Evaluate[
  Table[r[L][\[Lambda]], {L, 10, 100, 10}] /. soln], {\[Lambda], 900, 
  1100}]

Plot[r[59.9][\[Lambda]] /. soln, {\[Lambda], 900, 1100}]
Plot[r[60][\[Lambda]] /. soln, {\[Lambda], 900, 1100}]
Plot[r[60.1][\[Lambda]] /. soln, {\[Lambda], 900, 1100}]
Plot[r[49.8][\[Lambda]] /. soln, {\[Lambda], 900, 1100}]
Plot[r[50][\[Lambda]] /. soln, {\[Lambda], 900, 1100}]
Plot[r[50.2][\[Lambda]] /. soln, {\[Lambda], 900, 1100}]

The solutions behave as expected for most values of L, but for L=50,L=60, and possibly some other isolated values, the plot behaves badly. I'm very confused since the result should vary smoothly with $L$. I've solved this same differential equation in a different manner (Michael E2's solution in the question linked) and did not observe any erratic behavior, so I suspect it has to do with Method -> {"EquationSimplification" -> "Residual"}. I was wondering why this is happening, and how to fix it if possible.

Another strange occurrence is when I replace "Residual" with "Solve" and run the code, abort it, and re-replace "Solve" with "Residual", the code doesn't run. But after I save and run, it works fine. What's going on here?

Answer 1

In DAEs, you can give explicit initial conditions for all variables, including the highest order derivative. Since the IDA method uses a solver such as Newton's method, the solver can be inconsistent in the root it finds.

Clear[f, ivs];
f[r_] := 1 - 1/r;

ivs[r0_, L_] = {r[0], r'[0]} /.   (* use Set[] to pick an explicit formula *)
   First@Solve[
     {r'[λ]^2 == 1 - L^2/r[λ]^2*f[r[λ]], 
       r[0] == r0} /. {λ -> 0}, {r[0], r'[0]}];
soln = ParametricNDSolve[
  {r'[λ]^2 == 1 - L^2/r[λ]^2*f[r[λ]],
   {r[0], r'[0]} == ivs[1000, L]},   (* Beware: ParametericNDSolve[] evaluates ivs[]
                                         It needs to return a formula that depends on L *)
  r, {λ, 0, 2000}, {L}, 
  Method -> {"EquationSimplification" -> "Residual"}]

Plot[Evaluate[Table[r[L][λ], {L, 10, 100, 10}] /. soln], {λ, 900, 1100}]

Answer 2

Notice the problem actually has 2 solutions. When using the default option (in this case it's "EquationSimplification" -> "Solve" as far as I know), Mathematica will notice both solutions, though it has difficulty in finishing calculating one of them:

NDSolve[{r'[λ]^2 == 1 - L^2/r[λ]^2*f[r[λ]], r[0] == 1000} /. 
  L -> 10, r, {λ, 0, 2000}]

NDSolve::mxst

When you set "EquationSimplification" -> "Residual", you're choosing a DAE solver, which will find only one possible solution for the problem. You found the solution sometimes behaves "badly" just because the solution being found is not the desired one.

To make NDSolve find the desired solution, we need to lead NDSolve a bit. Micheal E2 has already shown you one of them. Here's another solution by adjusting the "DefaultStartingValue" option:

f[r_] = 1 - 1/r;

soln = ParametricNDSolve[{r'[λ]^2 == 1 - L^2/r[λ]^2*f[r[λ]], 
    r[0] == 1000}, r, {λ, 0, 2000}, {L}, 
   Method -> {"EquationSimplification" -> "Residual", 
     "DAEInitialization" -> {"Collocation", "DefaultStartingValue" -> -1}}];

Plot[Evaluate[Table[r[L][λ], {L, 10, 100, 10}] /. soln], {λ, 900, 1100}]

"DefaultStartingValue" -> -1 is amount to using $r'(0)=-1$ as the starting value. For more information, try struggling through the obscure tutorial linked by Michael E2 in the comment above.