# Unexpected/unphysical behavior of ODE near singularity

I'm trying to solve an ODE that becomes stiff or singular at a point. The ODE and initial conditions are

\begin{align*} \frac{d^2 z}{dt^2} &= f(z) f'(z) \bigg(1-\frac{3 \eta^2}{2} f(z) \bigg)\\ \frac{dz}{dt}\Bigg|_{t=t_0} &= \sqrt{1-\eta^2}\\ z\big|_{t=t_0}&=0. \end{align*}

where $f(z) = 1+\frac{z^2}{l^2} - \frac{\mu^d z^d}{l^{2d}}$. Here, $\mu$, $d$, and $l$ are constants, $\eta$ is an adjustable parameter between $0$ and $1$, and $t_0$ denotes the initial time.

I've implemented the system in the code snippet below:

d = 3;
l = 5.; μ = 26.^(1/d);

f[z_] := 1 + z^2/l^2 - μ^d z^d/l^(2 d);

Subscript[θ, 0] = 2.4;
tSoln = ParametricNDSolve[{ZZ''[t] == f[ZZ[t]] f'[ZZ[t]] (1 - (3 η^2)/2 f[ZZ[t]]),
ZZ'[-l Subscript[θ, 0]] == Sqrt[1 - η^2],
ZZ[-l Subscript[θ, 0]] == 0}, {ZZ},
{t, -l Subscript[θ, 0], l Subscript[θ, 0]}, {η}]
z[t_, η_] := ZZ[η][t] /. tSoln;

Plot[z[t, 0], {t, -5*2.4, 5*2.4}]


From the underlying physics, I know that for sufficiently small $\eta$, the function $z(t)$ should asymptotically approach a constant value (here $z \to 25$). However, I believe the system becomes stiff/singular as $z(t)$ approaches that value, and Mathematica doesn't find the correct solution.

As an example, the above code gives me this plot of $z(t)$: For some reason, the graph dips downward. (With other small, nonzero values of $\eta$, the graph crosses the asymptote of $z=25$.) If I take the $\eta=0$ case and integrate it directly rather than via an ODE solver, I obtain the behavior that I expect: Is there an ODE solution method that is particularly suited for this problem? I've mixed and matched with "BDF", "StiffnessSwitching", and some others without really understanding how they work, and none of them give consistent results. How should I approach this?

In short: How can I get Mathematica to understand the correct asymptotic behavior of a differential equation when it becomes stiff or singular?

Sorry for such a lengthy post, it's kind of hard for me to articulate the exact problem. I'm also not very familiar with the advanced tools of Mathematica, so I very much appreciate any guidance. Thank you so much! (Please let me know if I can clarify anything.)

It's just a matter of precision. WorkingPrecision -> 32 will resolve your problem:

d = 3;
l = 5; μ = 26^(1/d);

f[z_] = 1 + z^2/l^2 - μ^d z^d/l^(2 d);

Subscript[θ, 0] = 24/10;
tSoln = ParametricNDSolve[{ZZ''[t] == f[ZZ[t]] f'[ZZ[t]] (1 - (3 η^2)/2 f[ZZ[t]]),
ZZ'[-l Subscript[θ, 0]] == Sqrt[1 - η^2],
ZZ[-l Subscript[θ, 0]] == 0}, {ZZ}, {t, -l Subscript[θ, 0],
l Subscript[θ, 0]}, {η}, WorkingPrecision -> 32];
z[t_, η_] = ZZ[η][t] /. tSoln;

Plot[z[t, 0], {t, -5*2.4, 5*2.4}] • Cool! That does seem to be the case in this instance. However, I'm a bit hesitant to accept this answer because it doesn't seem to work as well for different values of $d$. For example, if I run the code as shown with $d=3$, the plot outputs immediately and exhibits the correct behavior. However, with $d=4$ or $d=5$, these are no longer true. Would I just need more WorkingPrecision for these cases? Is there a more "universal" method that takes known asymptotic behavior into account and works roughly equally for different cases of $d$? – pianyon Sep 20 '17 at 21:32
• @qm-arv Choosing WorkingPrecision -> 48 together with Method -> ExplicitRungeKutta (or the calculation will be too slow) seems to handle $d=4$ and $5$. I guess this behavior might be similar to the one mentioned in this post, if so, there doesn't exist a more general solution for so long, as far as I know. Anyway, feel free to wait for a moment, maybe someone more knowledgable will come up with a better answer :) . – xzczd Sep 21 '17 at 1:29