Why "Repeated convergence test failure" disappears in rare cases while occurring in most runs?

Question

I am trying to solve huge ordinary differential equation of second order using ParametricNDSolve. My code worked fine until I moved to a new model of plasma pressure. Now it exits with the following message:

ParametricNDSolveValue::ndcf: Repeated convergence test failure at z == 0.00007180340719427522`; unable to continue.

The most frustrating thing is that in rare cases it works without errors with the same set of input parameters.

At first, I thought that such random behavior of the code could be caused by lack of memory. However, testing the code on three Windows 10 (Professional or Corporate) computers with different memory sizes (from 16 GB to 64 GB) did not reveal any differences in the behavior of the code.

Sometimes it was possible to bring the calculation to the end without errors with a clean start of Wolfram Mathematica. But a second or third run of the calculation, even in this case, ended with the error. Clearing the system cache does not help.

I am attaching the most concise version of the code that reproduces the error.

Preparing equation for ParametricNDSolve

(* Coefficient A1 for ODE *)
A1[bv_, p_] =
  1/(2 bv) (1/3 + 1/(6 p) + ((1 - I Sqrt[3]) (1 + 2 p)^2)/(
     12 p (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
        p (-6 + 6 I Sqrt[3] bv Sqrt[
            1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(1/3)) + 
     1/(12 p) (1 + I Sqrt[3]) (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
        p (-6 + 6 I Sqrt[3] bv Sqrt[
            1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(
      1/3)) (2/(1/3 + 1/(
        6 p) + ((1 - I Sqrt[3]) (1 + 2 p)^2)/(12 p (-1 + 
             6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
             p (-6 + 6 I Sqrt[3] bv Sqrt[
                 1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(1/3)) + 
        1/(12 p) (1 + I Sqrt[3]) (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
           p (-6 + 6 I Sqrt[3] bv Sqrt[
               1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(
         1/3)) - (bv (18 bv p - (27 I Sqrt[3] bv^2 p^2)/Sqrt[
          1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3] + 
          I Sqrt[3] Sqrt[
           1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 
            8 p^3]) (I (I + Sqrt[3]) (1 + 2 p)^2 + (1 + 
             I Sqrt[3]) (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
             p (-6 + 6 I Sqrt[3] bv Sqrt[
                 1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(
           2/3)))/(3 (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
          p (-6 + 6 I Sqrt[3] bv Sqrt[
              1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(
        4/3) (1/3 + 1/(
          6 p) + ((1 - I Sqrt[3]) (1 + 2 p)^2)/(12 p (-1 + 
               6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
               p (-6 + 6 I Sqrt[3] bv Sqrt[
                   1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(1/3)) + 
          1/(12 p) (1 + I Sqrt[3]) (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 +
               p (-6 + 
                 6 I Sqrt[3] bv Sqrt[
                  1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(1/3))^2));

====

(* Coefficient B1 for ODE *)
B1[bv_, p_] =
  1/(4 bv^2)
    p (2/3 - 1/(6 p) + (I (I + Sqrt[3]) (1 + 2 p)^2)/(
     12 p (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
        p (-6 + 6 I Sqrt[3] bv Sqrt[
            1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(1/3)) - 
     1/(12 p) (1 + I Sqrt[3]) (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
        p (-6 + 6 I Sqrt[3] bv Sqrt[
            1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(1/3)) (1/3 + 1/(
     6 p) + ((1 - I Sqrt[3]) (1 + 2 p)^2)/(
     12 p (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
        p (-6 + 6 I Sqrt[3] bv Sqrt[
            1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(1/3)) + 
     1/(12 p) (1 + I Sqrt[3]) (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
        p (-6 + 6 I Sqrt[3] bv Sqrt[
            1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(1/3)) (4/3 + 1/(
     6 p) + ((1 - I Sqrt[3]) (1 + 2 p)^2)/(
     12 p (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
        p (-6 + 6 I Sqrt[3] bv Sqrt[
            1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(1/3)) + 
     1/(12 p) (1 + I Sqrt[3]) (-1 + 6 (-2 + 9 bv^2) p^2 - 8 p^3 + 
        p (-6 + 6 I Sqrt[3] bv Sqrt[
            1 + 6 p - 3 (-4 + 9 bv^2) p^2 + 8 p^3]))^(1/3));

======

(* Coefficient C1 for ODE *)
C1[z_] = (
  3 Derivative[1][bv][z]^2 - 2 bv[z] (bv^\[Prime]\[Prime])[z])/(
  4 bv[z]^2);

=======

(* Rule to put actual magnetic field profile *)
BvRule2[q_, M_, R_] = {bv -> ((1 + (M - 1) Sin[(\[Pi] #1)/2]^q)/R &)};
(* Right-End point of the interval where we solve ODE, bv[zR]==1 *)
zR[R_, {q_, M_}] = (2 ArcSin[((-1 + R)/(-1 + M))^(1/q)])/\[Pi];
(* derivative bv'[z] at z==zR;bv[zR]==1 
*)
bv1R[R_, {q_, M_}] = (\[Pi] q Sqrt[-1 + ((-1 + M)/(-1 + R))^(2/q)] (-1 + R))/(
  2 R);

========

(* ODE in symbolic form *)
ODE = D[(\[CapitalLambda] + 1 - 2 B1[bv[z], p]) D[\[Phi][z], z], 
    z] + \[Phi][
     z] (-D[Derivative[1][bv][z] A1[bv[z], p] (1 - B1[bv[z], p]), 
       z]) + \[Phi][
     z] (-2 B1[bv[z], p]*
       C1[z] - (1/2) Derivative[1][bv][z]^2 A1[bv[z], 
        p]^2 (1 - B1[bv[z], p]));

======

(* Take real part of the ODE coefficients since they can have small imaginary \
parts due to inaccuracy in float arithmetics *)
ODE2 = Coefficient[ODE, (\[Phi]^\[Prime]\[Prime])[z]];
ODE1 = Coefficient[ODE, Derivative[1][\[Phi]][z]];
ODE0 = Coefficient[ODE, \[Phi][z]];
ode2[z_, p_, \[CapitalLambda]_] = 
  Re@ODE0 *\[Phi][z] + Re@ODE1*Derivative[1][\[Phi]][z] + 
   Re@ODE2*(\[Phi]^\[Prime]\[Prime])[z];
Short[ode2[z, p, \[CapitalLambda]], 20];
"Coefficients of ODE are singular at bv[z]=p=1"
Echo[ODE0 /. {bv[z] -> 1, p -> 1} // Quiet, "ODE0="];
Echo[ODE1 /. {bv[z] -> 1, p -> 1} // Quiet, "ODE1="];
Echo[ODE2 /. {bv[z] -> 1, p -> 1} // Quiet // Simplify, "ODE2="];
ClearAll[ODE0, ODE1, ODE2];

=========

(* Compose ODE in the final form *)
ode2$[z_, p_, {q_, M_, R_}, \[CapitalLambda]_] = 
  ode2[z, p, \[CapitalLambda]] /. BvRule2[q, M, R];

==============

(* input parameters *)
qm = 2; Km = 8; Rm = 3; 
\[CapitalLambda]s = {1, 1.1, 2, 5, 10};
ps = {p00, pmin, pmax} = {0.25, 0.1, 0.99999};

debug = 2; xRootDebug = 0;

========

(* ODE is singular at bv==1 (i.e. at z=zR) and p==1,  *)
(* ****** add numerical factor <1 to exclude singularity at the right-end \
point ***** *)
zR$ = zR[Rm, {qm, Km}];
zR$ = (1 - 10^-6) zR[Rm, {qm, Km}

]

============

(* compose coefficient for BC at the right end *)
If[debug > 1, Print["\!\(\*
StyleBox[\"BC1$\",\nFontColor->RGBColor[1, 0, 0]]\)"]];
BC1$[p_, \[CapitalLambda]_] = (1/(\[CapitalLambda] + 1)) Re@A1[1, p]*
    bv1R[Rm, {qm, Km}] /. BvRule2[qm, Km, Rm];
If[debug > 1, Print[Short[BC1$[p, \[CapitalLambda]], 10]]];

Solve ParametricNDSolve

(* some chunk of main code, borrowed from a Do cycle over q, M, and R *)
aborted =.;
Check[
 Echo[tmp0 = Timing[pf2 = ParametricNDSolveValue[
       {ode2$[z, p, {qm, Km, Rm}, \[CapitalLambda]] == 0, \[Phi][0] ==
          1, \[Phi]'[0] == 0}
       , {\[Phi]'[zR$], \[Phi][zR$]}
       , {z, 0, (1. + 10^-7) zR$}
       , {{p, 0, 1/2}, {\[CapitalLambda], 1, 500}}
       , Method -> {"EquationSimplification" -> "Residual"}
       ]]
   ];
 ,
 Print["Error at evaluating ParametricNDSolveValue."]; aborted = True;
  Abort[];
 ]

=============

{0.015625,ParametricFunction[Expression: {([Phi]^[Prime])[0.359],[Phi][0.359]} Parameters: {p,[CapitalLambda]}

]}

=============

Print["Going to check BC1$ and pf2 for ps=", ps];
If[debug > 0,
  Do[
    Check[
      Print[
        StringForm[
         "pm=``, \[CapitalLambda]m=``, zR$=``, \
BC1[pm,\[CapitalLambda]m]=``, pf2[pm,\[CapitalLambda]m]=``"
         , pm, \[CapitalLambda]m, zR$, BC1$[pm, \[CapitalLambda]m], 
         pf2[pm, \[CapitalLambda]m]]];
      ,
      Print["Error at evaluating BC1$ and pf2 for p=", pm]; 
      aborted = True; Abort[];
      ];
    , {pm, ps}, {\[CapitalLambda]m, \[CapitalLambda]s}
    ];
  ];

=================

During evaluation of In[168]:= Going to check BC1$ and pf2 for ps={0.25,0.1,0.99999}
>     
>     During evaluation of In[168]:= ParametricNDSolveValue::ndcf: Repeated convergence test failure at z == 0.00007180334257120158`;
> unable to continue.
>     
>     During evaluation of In[168]:= InterpolatingFunction::dmval: Input value {0.359017} lies outside the range of data in the interpolating
> function. Extrapolation will be used.
>     
>     During evaluation of In[168]:= InterpolatingFunction::dmval: Input value {0.359017} lies outside the range of data in the interpolating
> function. Extrapolation will be used.
>     
>     During evaluation of In[168]:= pm=0.25`, \[CapitalLambda]m=1, zR$=0.3590166769543402`, BC1[pm,\[CapitalLambda]m]=-0.551922,

pf2[pm,[CapitalLambda]m]={-2.1293110^13+0. I,-1.5286910^12+0. I}

During evaluation of In[168]:= Error at evaluating BC1$ and pf2 for p=0.25
>     
>     Out[169]= $Aborted

====================

Print["Going to evaluate bcMode."];
Check[
  eqn2[p_?NumericQ, \[CapitalLambda]_?NumericQ] = 
    pf2[Min[p, 
       pmax], \[CapitalLambda]] . {1, -BC1$[
        Min[p, pmax], \[CapitalLambda]]};
  ,
  Print["Error at evaluating bcMode."]; aborted = True; Abort[]
  ];

If[debug > 0,
  Check[
    \[CapitalLambda]1 = First[\[CapitalLambda]s]; \[CapitalLambda]2 = 
     Last[\[CapitalLambda]s];
    Print["eqn2[p,\[CapitalLambda]]/.{p\[Rule]", p00, 
     ", \[CapitalLambda]\[Rule]", \[CapitalLambda]1, "}=", 
     eqn2[p, \[CapitalLambda]] /. {p -> 
        p00, \[CapitalLambda] -> \[CapitalLambda]1},
     "\neqn2[p,\[CapitalLambda]]/.{p\[Rule]", p00, 
     ", \[CapitalLambda]\[Rule]", \[CapitalLambda]2, "}=", 
     eqn2[p, \[CapitalLambda]] /. {p -> 
        p00, \[CapitalLambda] -> \[CapitalLambda]2}
     ];
    , Print["Error at evaluating eqn2."]; aborted = True; Abort[];
    ];
  ];

---

During evaluation of In[66]:= Going to evaluate bcMode.

During evaluation of In[66]:= InterpolatingFunction::dmval: Input value {0.359017} lies outside the range of data in the interpolating function. Extrapolation will be used.

During evaluation of In[66]:= InterpolatingFunction::dmval: Input value {0.359017} lies outside the range of data in the interpolating function. Extrapolation will be used.

During evaluation of In[66]:= ParametricNDSolveValue::mconly: For the method NDSolve`IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.

During evaluation of In[66]:= InterpolatingFunction::dmval: Input value {0.359017} lies outside the range of data in the interpolating function. Extrapolation will be used.

During evaluation of In[66]:= General::stop: Further output of InterpolatingFunction::dmval will be suppressed during this calculation.

During evaluation of In[66]:= eqn2[p,[CapitalLambda]]/.{p->0.25, [CapitalLambda]->1}=-2.2136810^13+0. I eqn2[p,[CapitalLambda]]/.{p->0.25, [CapitalLambda]->10}=9.1453410^51+0. I

During evaluation of In[66]:= Error at evaluating eqn2.

Out[68]= $Aborted

Reduced version

The coefficients of the equation are singular at bv==1 and p=1. And this could be the reason for the code to be unstable, although the code bypasses the singularity point, in particular by means of the following lines:

ps = {p00, pmin, pmax} = {0.25, 0.1, 0.99999};
zR$ = (1. - 10^-6) zR[Rm, {qm, Km}]

However, the error disappears if I use the reduced version of the coefficients A1 and B1, namely:

(* Coefficient A1 for ODE *)
A1b[bv_, p_] = -(1/Sqrt[6 (1 - bv) + (1 - p)^2]);
(* Coefficient B1 for ODE *)
B1b[bv_, p_] = 1/6 (Sqrt[6 (1 - bv) + (1 - p)^2] - (1 - p));

They are also singular in the same point. Here is result of successful computation:

> Going to check BC1$ and pf2 for ps={0.25,0.1,0.99999}
> 
> pm=0.25`, \[CapitalLambda]m=1, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-2.20769,
> pf2[pm,\[CapitalLambda]m]={-1.86753,0.683703}
> 
> pm=0.25`, \[CapitalLambda]m=1.1`, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-2.10256,
> pf2[pm,\[CapitalLambda]m]={-1.79862,0.70001}
> 
> pm=0.25`, \[CapitalLambda]m=2, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-1.47179,
> pf2[pm,\[CapitalLambda]m]={-1.33924,0.795023}
> 
> pm=0.25`, \[CapitalLambda]m=5, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-0.735895,
> pf2[pm,\[CapitalLambda]m]={-0.712639,0.900225}
> 
> pm=0.25`, \[CapitalLambda]m=10, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-0.401398,
> pf2[pm,\[CapitalLambda]m]={-0.398732,0.946213}
> 
> pm=0.1`, \[CapitalLambda]m=1, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-1.83974,
> pf2[pm,\[CapitalLambda]m]={-1.62503,0.703965}
> 
> pm=0.1`, \[CapitalLambda]m=1.1`, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-1.75213,
> pf2[pm,\[CapitalLambda]m]={-1.56212,0.719155}
> 
> pm=0.1`, \[CapitalLambda]m=2, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-1.22649,
> pf2[pm,\[CapitalLambda]m]={-1.15156,0.807825}
> 
> pm=0.1`, \[CapitalLambda]m=5, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-0.613246,
> pf2[pm,\[CapitalLambda]m]={-0.60705,0.906316}
> 
> pm=0.1`, \[CapitalLambda]m=10, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-0.334498,
> pf2[pm,\[CapitalLambda]m]={-0.33842,0.949465}
> 
> pm=0.99999`, \[CapitalLambda]m=1, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-165576.,
> pf2[pm,\[CapitalLambda]m]={-91559.2,0.552935}
> 
> pm=0.99999`, \[CapitalLambda]m=1.1`, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-157692.,
> pf2[pm,\[CapitalLambda]m]={-90928.2,0.576608}
> 
> pm=0.99999`, \[CapitalLambda]m=2, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-110384.,
> pf2[pm,\[CapitalLambda]m]={-78695.1,0.712925}
> 
> pm=0.99999`, \[CapitalLambda]m=5, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-55192.2,
> pf2[pm,\[CapitalLambda]m]={-47532.9,0.861227}
> 
> pm=0.99999`, \[CapitalLambda]m=10, zR$=(2 ArcSin[Sqrt[2/7]])/\[Pi],
> BC1[pm,\[CapitalLambda]m]=-30104.8,
> pf2[pm,\[CapitalLambda]m]={-27858.9,0.92538}

The run is successful even if I do not exclude the point bv=1 by changing zR$:

zR$ = (1 - 0 10^-6) zR[Rm, {qm, Km}]

I would appreciate any idea or suggestion how to overcome the error ParametricNDSolveValue::ndcf.

In the Wolfram Mathematica help, when this error occurs, it is suggested to contact the technical support service of this product and send them a code sample that causes this error. But I doubt that they will fix the bug before the release of a new version of the product.

Answer 1

Try Method->StiffnessSwitching.