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}]
  • $\begingroup$ 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 $\endgroup$
    – george2079
    Commented Dec 29, 2017 at 18:16
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Dec 30, 2017 at 6:09

2 Answers 2


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

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

Mathematica graphics

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}]

Mathematica graphics

  • $\begingroup$ Gosh, these images look horrible today. They're from the SE Uploader, Image(pp). They don't look nearly so bad on my screen. $\endgroup$
    – Michael E2
    Commented 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]

enter image description here

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.