# NDSolveValue stiff system error

I have two equations involving elliptical integrals (equat1 and equat2 in the code below). My goal is to obtain plots of u[o] and d[o]. I turn my equations into ODE with initial conditions found using FindRoot. While trying to solve my ODE with NDSolveValue,I encounter the error:

NDSolveValue::ndsz: At o == 0.259506989709951, step size is effectively zero; singularity or stiff system suspected.

However, I expect a smooth curve as my solution. I tried to do a series expansion of the elliptical integrals and it worked. The solutions were smooth curve. However, as I include higher orders in the expansion, I encounter the same error again. I tried increasing the WorkingPrecesionand changing the Methodoption. Probably StiffnessSwitching will help, but I did not manage to make it work. Any ideas on how I should proceed?

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])));

tlambda1[u_, d_] := u/d;
tlambda2[u_, d_] := -u/d;

tk1sq[u_, d_] := (Sqrt[1 + ((tlambda1[u, d])^2)] +
tlambda1[u, d])/(2*(Sqrt[1 + ((tlambda1[u, d])^2)]));

tk2sq[u_, d_] := (Sqrt[1 + ((tlambda2[u, d])^2)] +
tlambda2[u, d])/(2*(Sqrt[1 + ((tlambda2[u, d])^2)]));

ta1[u_, d_] := (1 + (tlambda1[u, d])^2)^(1/4);
ta2[u_, d_] := (1 + (tlambda2[u, d])^2)^(1/4);
tb1[u_, d_] := Sqrt[1 + (tlambda1[u, d])^2] + tlambda1[u, d];
tb2[u_, d_] := Sqrt[1 + (tlambda2[u, d])^2] + tlambda2[u, 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]))) == (Sqrt[1]*(((tlambda1[u, d])*
EllipticK[tk1sq[u, d]]/(2*ta1[u, d])) - (ta1[u, d]*
EllipticE[tk1sq[u, d]]) + ((EllipticK[tk1sq[u, d]])/(2*ta1[u, d]*
tb1[u, d])) + ((tlambda2[u, d])*EllipticK[tk2sq[u, d]]/(2*ta2[u, d])) -
(ta2[u, d]*EllipticE[tk2sq[u, d]]) + ((EllipticK[tk2sq[u, d]])/(2*ta2[u, d]*
tb2[u, d]))));

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

{dsolj, usolj} = NDSolveValue[{D[equat2[u, o, d] /. {d -> d[o], u -> u[o]}, o],
D[equat1[u, o, d, .1] /. {d -> d[o], u -> u[o]}, o], u[0] == uj, d[0] == dj}, {d, u},
{o, 0, 1}, Method -> {"EquationSimplification" -> "Residual"}]

Plot[{dsolj[t], usolj[t]}, {t, 0, 1}, PlotRange -> {0, 1.5}]

• that appears to simply be the nature of the system, it is singular around .259. You can get it to run margially further with MaxStepSize -> .0001 Dec 29, 2017 at 18:16
• Shouldn't the question be how to determine whether there is stiffness or a singularity? If there's not supposed to be a singularity, but in fact there is one in the coded system, then that seems an important clue. Dec 30, 2017 at 6:09

There's a singularity in the ODE. Perhaps you're thinking d as a function of u is smooth. Here's a picture of the vector field of the ODE. The solution runs right up to where the vector changes direction (turns from red to green). It seems that as a direction field, the field is smooth, but the vector changes orientation because the derivative of o is always 1 and the o-d components rotate past vertical (in the projection onto the o-d plane).

sys = {D[equat2[u, o, d] /. {d -> d[o], u -> u[o]}, o],
D[equat1[u, o, d, .1] /. {d -> d[o], u -> u[o]}, o], u[0] == uj,
d[0] == dj}

phasefield = First@Solve[sys[[{1, 2}]], D[Through[{u, d}[o]], o]];

vplot = VectorPlot3D[
{1, d'[o], u'[o]} /. phasefield /. (f_)[o] :> f // Evaluate,
{o, -0.03, 1}, {d, 0.065, 1.}, {u, 0., 0.15},
VectorScale -> {Small, 0.2, 1 &}, AxesLabel -> {o, d, u},
VectorColorFunctionScaling -> False,
VectorColorFunction ->
Function[{o, d, u, vo, vd, vu, n},
Darker@Hue[1/2 + ArcTan[vd, vu]/(2 Pi)]]];

Show[
vplot,
ParametricPlot3D[{o, dsolj[o], usolj[o]},
{o, 0, 0.2595069897099519},
PlotStyle -> {AbsoluteThickness[4], Blue}]
]


Here is a slice of the vector field showing the reversal of the vectors at the point where NDSolve stopped integrating.

Show[vplot, ParametricPlot3D[{o, dsolj[o], usolj[o]},
{o, 0, 0.2595069897099519},
PlotStyle -> {AbsoluteThickness[4], Blue}],
PlotRange -> {All, All, {0.085, 0.13}}, BoxRatios -> {1, 1, 0.1}]


• Gosh, these images look horrible today. They're from the SE Uploader, Image(pp). They don't look nearly so bad on my screen. Dec 30, 2017 at 6:07

Instead of using NDSolve applied to the total derivatives of equat1 and equat2, as in the question, one can solve the two equations simultaneously with FindRoot. The result is

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


which is indistinguishable from the NDSolve results. Simply stated, FindRoot can find no real roots for o > 0.259507154261615. It is for this reason that NDSolve, as employed in the question, fails at

dsolj["Domain"][[1, 2]]
(* 0.259506989709951 *)


Incidentally, experimentation suggests that there are no real roots for larger o, at least as far as o == 1.