# Runtime too long for NDSolveValue, FindRoot breaks down at sharp turns

I have two equations:

$$\frac{2}{3}\thinspace\left(\frac{\varepsilon}{\Delta}\right)^{\frac{3}{2}}=\sum_{i=1,2} \sigma_i\lambda_{i}\left(1+\lambda_{i}^{2}\right)^{\frac{1}{4}}E\left(\frac{\pi}{2},k_{i}\right)+\sum_{i=1,2} \sigma_i \frac{(1+\lambda_i^2)^{1/4} }{ 2 ( \lambda_{i}+\sqrt{1+\lambda_{i}^{2} })} F(\frac{\pi}{2}, k_i )$$

and

$$\sqrt{\Delta}\sum_{i=1,2}\biggl[-\left(1+\lambda_{i}^{2}\right)^{\frac{1}{4}}E\left(\frac{\pi}{2},k_{i}\right)+\frac{F\left(\frac{\pi}{2},k_{i}\right)}{2\left(1+\lambda_{i}^{2}\right)^{\frac{1}{4}}}\left( \lambda_{i} + \frac{1}{\sqrt{1+\lambda_{i}^{2}}+\lambda_{i}} \right)\biggr]=\sum_{i=1,2}\biggl[-\left(1+\widetilde{\lambda}_{i}^{2}\right)^{\frac{1}{4}}E\left(\frac{\pi}{2},\widetilde{k}_{i}\right)+\frac{F\left(\frac{\pi}{2},\widetilde{k}_{i}\right)}{2\left(1+\widetilde{\lambda}_{i}^{2}\right)^{\frac{1}{4}}}\left(\widetilde{\lambda}_{i}+\frac{1}{\sqrt{1+\widetilde{\lambda}_{i}^{2}}+\widetilde{\lambda}_{i}} \right)\biggr]$$

where $$\sigma_1=1$$, $$\sigma_2=-1$$, $$k_{i}^{2}=\frac{\sqrt{1+\lambda_{i}^{2}}+\lambda_{i}}{2\sqrt{\left(1+\lambda_{i}^{2}\right)}}$$ for $$i=1,2$$ with $$\lambda_{1}\equiv \frac{\mu-\Delta_{0}}{\Delta}$$, $$\lambda_{2}\equiv \frac{-\Delta_{0}-\mu}{\Delta}$$, $$\widetilde{k}_{i}^{2}=\frac{\sqrt{1+\widetilde{\lambda}_{i}^{2}}+\widetilde{\lambda}_{i}}{2\sqrt{\left(1+\widetilde{\lambda}_{i}^{2}\right)}}$$ for $$i=1,2$$ with $$\widetilde{\lambda}_{1}\equiv \mu$$, $$\widetilde{\lambda}_{2}\equiv -\mu$$. $$F(\frac{\pi}{2},k)$$ and $$E(\frac{\pi}{2},k)$$ are the complete elliptic integral of the first and the second kind respectively which in Mathematica are EllipticK[] and EllipticE[]. In the code below $$\mu$$, $$\Delta_0$$, $$\Delta$$, $$\varepsilon$$ are represented by u, o, d, e.

By first converting those two equations in the form of differential equations and then using NDSolveValue to obtain plots of $$\Delta$$ and $$\mu$$ as a function of $$\Delta_0$$ at $$\varepsilon=0.001$$, I encounter the error:

lambda1[u_, o_, d_] := (u - o)/d;
lambda2[u_, o_, d_] := (-u - o)/d;
k1sq[u_, o_, d_] := (Sqrt[1 + ((lambda1[u, o, d])^2)] + lambda1[u, o, d])/(2*(Sqrt[1 + ((lambda1[u, o, d])^2)]));
k2sq[u_, o_, d_] := (Sqrt[1 + ((lambda2[u, o, d])^2)] + lambda2[u, o, d])/(2*(Sqrt[1 + ((lambda2[u, o, d])^2)]));
a1[u_, o_, d_] := (1 + (lambda1[u, o, d])^2)^(1/4);
a2[u_, o_, d_] := (1 + (lambda2[u, o, d])^2)^(1/4);
b1[u_, o_, d_] := Sqrt[1 + (lambda1[u, o, d])^2] + lambda1[u, o, d];
b2[u_, o_, d_] := Sqrt[1 + (lambda2[u, o, d])^2] + lambda2[u, o, d];

equat1[u_, o_, d_, e_] =
(2/3)*(((e)/d)^(3/2)) == ((lambda1[u, o, d]*a1[u, o, d]*EllipticE[k1sq[u, o, d]])
-((lambda1[u, o, d]*EllipticK[k1sq[u, o, d]])/(2*a1[u, o, d]*b1[u, o, d]))+
((EllipticK[k1sq[u, o, d]])/(2*a1[u, o, d])) -
(lambda2[u, o, d]*a2[u, o, d]* EllipticE[k2sq[u, o, d]])
+ ((lambda2[u, o, d]*EllipticK[k2sq[u, o, d]])/(2*a2[u, o, d]*b2[u, o, d]))-
((EllipticK[k2sq[u, o, d]])/(2*a2[u, o, d])));

equat2[u_, o_, d_] =
Sqrt[d]*(((lambda1[u, o, d])*
EllipticK[k1sq[u, o, d]]/(2*a1[u, o, d])) - (a1[u, o, d]*
EllipticE[
k1sq[u, o, d]]) + ((EllipticK[k1sq[u, o, d]])/(2*a1[u, o, d]*
b1[u, o, d])) + ((lambda2[u, o, d])*
EllipticK[k2sq[u, o, d]]/(2*a2[u, o, d])) - (a2[u, o, d]*
EllipticE[
k2sq[u, o, d]]) + ((EllipticK[k2sq[u, o, d]])/(2*a2[u, o, d]*
b2[u, o, d]))) == ((lambda1[u, 0, 1])*
EllipticK[k1sq[u, 0, 1]]/(2*a1[u, 0, 1])) - (a1[u, 0, 1]*
EllipticE[
k1sq[u, 0, 1]]) + ((EllipticK[k1sq[u, 0, 1]])/(2*a1[u, 0, 1]*
b1[u, 0, 1])) + ((lambda2[u, 0, 1])*
EllipticK[k2sq[u, 0, 1]]/(2*a2[u, 0, 1])) - (a2[u, 0, 1]*
EllipticE[
k2sq[u, 0, 1]]) + ((EllipticK[k2sq[u, 0, 1]])/(2*a2[u, 0, 1]*
b2[u, 0, 1]));

{dk, uk} = {d, u} /. FindRoot[{equat2[u, o, d], equat1[u, o, d, .001]} /. o -> 0, {{d, .5}, {u, .1}}];

{dsolk, usolk} =
NDSolveValue[{D[equat2[u, o, d] /. {d -> d[o], u -> u[o]}, o],
D[equat1[u, o, d, .001] /. {d -> d[o], u -> u[o]}, o],
u[0] == uk, d[0] == dk}, {d, u}, {o, 0, 1},
Method -> {"EquationSimplification" -> "Residual"}];

Plot[{dsolk[t], usolk[t]}, {t, 0, 1}, PlotRange -> {-.1, 1.1}, Frame -> True]

NDSolveValue::mxst: Maximum number of 100000 steps reached at the point o == 0.6097033678408679`.

I expect the blue curve ($$\Delta$$) to continue levelling off and the yellow curve ($$\mu$$) to continue increasing. Increasing MaxSteps->Infinity probably helps but the run-time took too long to evaluate and I gave up. Eventually, I will be obtaining these plots for different values of $$\varepsilon$$ and so increasing the run-time by increasing MaxSteps does not seem to be the way to go.

Using FindRoot directly gives somwhat the plot I want and the run-time is relatively shorter, but fails around the sharp turn at $$\Delta_0=0.29$$ and throws the error:

Table[Values@
FindRoot[{equat2[u, o, d],
equat1[u, o, d, .001]}, {{d, .5}, {u, .1}}], {o, 0, 1, .001}];
ListPlot[Transpose@%, DataRange -> {0, 1}, Joined -> True,
AxesLabel -> {o, "d, u"},
LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

I think the problem is a 0/0 term somewhere in the expression for equat2[] that arises at the time the integration begins to fail, from around 0.60 to 0.62. This causes numeric noise in computing u'[o], which leads to NDSolve[] thinking there's an error. A common factor can be eliminated for o -> 0 via Simplify. I let Simplify run for a several minutes on general o to no avail. So I resorted to Limit, which is slow. To resolve the issue, I needed higher WorkingPrecision, so I had to avoid Method -> {"EquationSimplification" -> "Residual"}, which can only be used with the machine-precision IDA method. I also replaced the parameter e value 0.001 with 1/1000 to avoid machine precision. I timed all the time-consuming computations.

{dk, uk} = {d, u} /.
FindRoot[{equat2[u, o, d], equat1[u, o, d, 1/1000]} /.
o -> 0, {{d, .5}, {u, .1}}, WorkingPrecision -> 100];

newODE = First@Solve[
{D[equat2[u, o, d] /. {d -> d[o], u -> u[o]}, o],
D[equat1[u, o, d, 1/1000] /. {d -> d[o], u -> u[o]}, o]},
{d'[o], u'[o]}
] /. Rule -> Equal; // AbsoluteTiming
(*  {7.57532, Null}  *)

ClearAll[uprime];
uprime[o1_?NumericQ, d1_, u1_] := Block[{o},
Quiet[Check[
newODE[[2, 2]] /. {d[o] -> d1, u[o] -> u1, o -> o1},
nlim++;
Limit[newODE[[2, 2]] /. {d[o] -> d1, u[o] -> u1}, o -> o1],
{Power::infy, Infinity::indet}],
{Power::infy, Infinity::indet}]
];

newODE2 = {First@newODE, u'[o] == uprime[o, d[o], u[o]]};

nlim = 0;
{dsolk, usolk} =
NDSolveValue[{newODE2, u[0] == uk, d[0] == dk}, {d, u}, {o, 0, 1},
PrecisionGoal -> 8, AccuracyGoal -> 8,
WorkingPrecision -> 50]; // AbsoluteTiming
nlim
(*
{99.1584, Null}
3  <-- number of times Limit[] was called
*)

ListLinePlot@{dsolk, usolk}

In sum, equat2[] is numerically ill-conditioned in places. Perhaps simplification can help, or manual analysis, although the equations are too unwieldy for a quick look-see. However, I don't have any more time to look into it or the FindRoot problem, which is likely to be related, at this point.

• Note also that the total number of steps according to dsolk@"Grid" // Length is only 264. Jul 1, 2019 at 0:34