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.)