Mfcttrf

I have two equations:

$$
frac23thinspaceleft(fracvarepsilonDeltaright)^frac32=sum_i=1,2 sigma_ilambda_ileft(1+lambda_i^2right)^frac14Eleft(fracpi2,k_iright)+sum_i=1,2 sigma_i frac(1+lambda_i^2)^1/4 2 ( lambda_i+sqrt1+lambda_i^2 ) F(fracpi2, k_i )
$$

and

$$
sqrtDeltasum_i=1,2biggl[-left(1+lambda_i^2right)^frac14Eleft(fracpi2,k_iright)+fracFleft(fracpi2,k_iright)2left(1+lambda_i^2right)^frac14left( lambda_i +
frac1sqrt1+lambda_i^2+lambda_i right)biggr]=sum_i=1,2biggl[-left(1+widetildelambda_i^2right)^frac14Eleft(fracpi2,widetildek_iright)+fracFleft(fracpi2,widetildek_iright)2left(1+widetildelambda_i^2right)^frac14left(widetildelambda_i+frac1sqrt1+widetildelambda_i^2+widetildelambda_i right)biggr]
$$

where $sigma_1=1$, $sigma_2=-1$,
$k_i^2=fracsqrt1+lambda_i^2+lambda_i2sqrtleft(1+lambda_i^2right)$ for $i=1,2$ with $lambda_1equiv fracmu-Delta_0Delta$, $lambda_2equiv frac-Delta_0-muDelta $, $widetildek_i^2=fracsqrt1+widetildelambda_i^2+widetildelambda_i2sqrtleft(1+widetildelambda_i^2right)$ for $i=1,2$ with $widetildelambda_1equiv mu$, $widetildelambda_2equiv -mu $. $F(fracpi2,k)$ and $E(fracpi2,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.