# Numerically solve 2nd order differential equation with singularity

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, y' == 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?

• 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.'s ennui 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.'s ennui 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.'s ennui 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

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]}. *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, y' == 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, y' == 1}, y, {s, 0, 2 EllipticK[7/10]}, Method -> "StiffnessSwitching", "ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}, WorkingPrecision ->$wp];


Plot:

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