Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Consider a second order differential equation with a potential that diverges at some generic value in the variable. For example:

$$-y^{\prime\prime}(s)+\frac1{\mathrm{cn}{(s\mid k^2)}}y(s)=0$$

where $\mathrm{cn}(s\mid k^2)$ is the JacobiCN[s, k^2] function. Now, suppose, the differential equation has to be solved on the interval of a full period, namely $s\in[0,2\mathbb{K}]$, where $\mathbb{K}$ is the complete elliptic integral of the first kind EllipticK[k^2], with the initial conditions $y(0)=0$ and $y^\prime(0)=1$.

If I just naïvely start the NDSolve function for, say, $k^2=0.7$:

NDSolve[{-y''[s] + y[s]/JacobiCN[s, 0.7] == 0, y[0] == 0, y'[0] == 1},
  y[s], {s, 0, 2 EllipticK[0.7]}]

Mathematica rightfully complains that a singularity is detected at half of the evaluation range, namely at the value s = EllipticK[k^2] = 2.07536:

NDSolve::ndsz: At s == 2.0753631352632516`, step size is effectively zero; singularity or stiff system suspected.

And also the resulting interpolating function terminates at this value:

{{y[s] -> InterpolatingFunction[{{0., 2.07536}},<>][s]}}

So, evidently, the solution does not overcome the point of singularity and is given only for the half of the required interval. I heard that there are numerical techniques to overcome singularities when solving differential equations. It makes me wonder that Mathematica does not activate these by default. What can I do to solve a differential equation beyond a singularity numerically in Mathematica?

share|improve this question
This documentation about events might be helpful. –  Silvia May 9 '13 at 9:11
(I'm not posting an answer since this is rather specific to the DE concerned.) You have a coefficient that is a Jacobi elliptic function, which makes Picard's theory applicable (in complete analogy to the use of Floquet theory for differential equations with trigonometric coefficients, e.g. Mathieu). You might have better luck with an analytic approach instead of piggishly persisting with crossing through a singularity. –  J. M. is back. May 9 '13 at 10:17
One other possibility you can do if you're unable or unwilling to follow through with Picard analysis is to do a bit of "pole-vaulting"; that is, integrate along a path in the complex plane that avoids the poles of $\mathrm{nc}(u\mid m)$. –  J. M. is back. May 9 '13 at 12:36
"...involves a much more complicated, but still strictly elliptic potential term." - but you at least know where the poles might possibly be, yes? If so, your only job is to find a nice contour that avoids those nasty poles. In any event, if you want worked examples of Picard analysis, you might want to look up treatments of the well-studied Lamé equation. –  J. M. is back. May 9 '13 at 13:15
I will just quote "A powerful way of computing the twelve Jacobian elliptic functions for real or complex values of both the argument z and the modulus k is to use the definitions in terms of theta functions given in ..." as found in NIST Handbook of Mathematical Functions pg. 566 –  Spawn1701D May 9 '13 at 13:34

1 Answer 1

up vote 5 down vote accepted

The the function y is the integral of a logarithmic singularity, so it is relatively easy to understand. The function JacobiCN[s, 7/10] is relatively flat near s == EllipticK[7/10], so that the integral of its reciprocal can be approximated by a logarithm:

y[s]/JacobiCN[s, 7/10] ~= y0 / (jp (s - EllipticK[7/10]))

where y0 ~= y[s] at s == EllipticK[7/10] and jp = D[JacobiCN[s, 7/10], s] /. {s -> EllipticK[7/10]}.

Mathematica graphics *Plots of JacobiCN[s, 7/10] and the variation in its derivative. *

GraphicsRow@{Plot[JacobiCN[s, 7/10], {s, s1, s2}, 
   GridLines -> {{EllipticK[7/10]}, None}],
  Plot[D[JacobiCN[s, 7/10], s] - (D[JacobiCN[s, 7/10], s] /. s -> EllipticK[7/10]) // 
    Evaluate, {s, s1, s2}, WorkingPrecision -> $wp, 
   GridLines -> {{EllipticK[7/10]}, None}]}

We integrate once to find y[s] when the integration stops just before the singularity at s == s0. This is adjusted by integrating y''[s] ~= y0 / (jp (s - EllipticK[7/10])) to estimate y[s] at s == EllipticK[7/10].

The threshold for switching to a linear approximate of y[s]/JacobiCN[s, 7/10] is set by solving JacobiCN[s, 7/10] == eps for a small positive eps. This will yield two points s1, s2, before and after EllipticK[7/10] respectively. We then integrate to estimate the change in p[s] == y'[s] over the interval s1 < s < s2.

The displacements of the points s1 < s0 < s2 from EllipticK[7/10], 6] turn out to be

N[{s1, s0, s2} - EllipticK[7/10], 6]
(*  {-1.82574*10^-10, -1.78664*10^-12, 1.82574*10^-10}  *)

Here is the whole code:

$wp = 30;                             (* working precision ($)*)
  (* estimate y at singularity *)
{sol0} = NDSolve[
   {-y''[s] + y[s]/JacobiCN[s, 7/10] == 0, y[0] == 0, y'[0] == 1},
   y, {s, 0, EllipticK[7/10]},
   Method -> "StiffnessSwitching", 
   "ExtrapolationHandler" -> {Indeterminate &, 
     "WarningMessage" -> False}, WorkingPrecision -> MachinePrecision];
  (* data for slipping past singularity *)
eps = 1*^-10;                                                (* threshold for switching *)
jp = D[JacobiCN[s, 7/10], s] /. {s -> EllipticK[7/10]};      (* for lin. approx. *)
s0 = SetPrecision[y["Grid"] /. sol0 // Last // First, $wp];  (* end data of sol0 ($)*)
y0 = SetPrecision[y["ValuesOnGrid"] /. sol0 // Last, $wp];   (*   "  "   ($)*)
p0 = SetPrecision[y'["ValuesOnGrid"] /. sol0 // Last, $wp];  (*   "  "   ($)*)
s1 = InverseJacobiCN[eps, 7/10];                             (* interval for switching *)
s2 = InverseJacobiCN[-eps, 7/10];                            (*   "  "   *)
dy = Chop@NIntegrate @@                        (* estimate change in y to singularity *)
   N[{p0 + y0 (x Log[x] - x /. x -> EllipticK[7/10] - s)/jp, {s, s0, 
      EllipticK[7/10]}, PrecisionGoal -> 10, 
     WorkingPrecision -> $wp}, $wp]; 
dp = NIntegrate @@                             (* estimate mean p = y' to singularity *)
    N[{(y0 + dy) Log[Abs[s - EllipticK[7/10]]]/jp, {s, s1, 
       EllipticK[7/10], s2}, WorkingPrecision -> $wp}, $wp]/(s2 - s1);

  (* final integration *)
{sol1} = NDSolve[{
    y''[s] == Piecewise[{             (* switch p' = y'' at singularity *)
       {dp, s1 < s < s2}},              (* use average rate of change near singularity *)
      y[s]/JacobiCN[s, 7/10]],          (* default *)
    y[0] == 0, y'[0] == 1},
   y, {s, 0, 2 EllipticK[7/10]},
   Method -> "StiffnessSwitching", 
   "ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False},
   WorkingPrecision -> $wp];


Plot[y[s] /. {sol1} // Evaluate, {s, 0, 2 EllipticK[0.7]}]

Mathematica graphics

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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