# Numerical solution of IVP for linear ODE with variable coefficient runs wild soon

Cross posted in scicomp.SE.

A friend of mine showed me this initial value problem (IVP) for a linear ordinary differential equation (ODE) with variable coefficient:

$$y''(x)=\left(x^2-1\right) y(x)$$$$y(0)=1$$$$y'(0)=0$$

Seems to be a simple one, right? Actually it can be solved analytically by DSolve:

asol = DSolve[{y''[x] == (x^2 - 1) y[x], y[0] == 1, y'[0] == 0}, y[x], x]

{{y[x] -> E^(-(x^2/2))}}


But its numerical solution given by NDSolve runs wild very fast:

l = 10;
nsol = NDSolve[{y''[x] == (x^2 - 1) y[x], y[0] == 1, y'[0] == 0},
y, {x, 0, l}]; // AbsoluteTiming
Manipulate[Plot[y[x] /. nsol // Evaluate, {x, 0, l2},
PlotRange -> {{0, 7}, {-10, 1}}], {l2, 1/10, 7}]

{0.015600, Null}


How to resolve the problem?

Well, actually I found two solution for this problem but one is time-consuming and the other is limited so I'm not quite satisfied.

Solution 1

A higher WorkingPrecision will help:

l = 10;
nsol = NDSolve[{y''[x] == (x^2 - 1) y[x], y[0] == 1, y'[0] == 0},
y, {x, 0, l}, WorkingPrecision -> 50]; // AbsoluteTiming
Plot[y[x] /. nsol, {x, 0, l}, PlotRange -> All]

{2.857000, Null}


Nonetheless, this solution is slow, and will need a higher WorkingPrecision and be even slower when l gets larger.

Solution 2

Noticing the analytic solution involves a Exp, given the experience that Exp often causes trouble in numerical calculation, the transformation $y(x)=e^{z(x)}$ is used:

l = 50;
rule = y -> (Exp[z@#] &);

(* z[0] == 0 is manually substituted because it seems that
NDSolveValue has some difficulty in understanding y[0] == 1 /. rule *)
nsol = NDSolveValue[{y''[x] == (x^2 - 1) y[x], z[0] == 0, y'[0] == 0} /. rule,
z, {x, 0, l}]; // AbsoluteTiming
Plot[Exp@nsol[x], {x, 0, l}, PlotRange -> All]

 {0.013000, Null}


However, as mentioned above, this solution is limited, I wouldn't have thought it out if I was unaware of the analytic solution, and I really doubt if this method can be extended: I guess there's a sort of problem behind this specific example, once it's solved, more complicated problem might be solved too, for example this and this.

I'd appreciate if anyone can give an in-depth explanation for the problem or find a better solution.

I think this shows what the problem is. The ODE is inherently numerically unstable as any deviation from the exact solution goes to ±∞.

asol = DSolve[{y''[x] == (x^2 - 1) y[x], y[0] == 1, y'[0] == p}, y[x], x];
Manipulate[
Plot[y[x] /. asol /. {{p -> 0}, {p -> 10^(-error)}, {p -> -10^(-error)}} // Evaluate,
{x, 0, 7}, PlotRange -> 2],
{error, 1, 10, Appearance -> "Labeled"}]


Update - Possible workaround

Aside from increasing the working precision and increasing the "DifferenceOrder", one can try the approach I took toward the Airy equation in my answer to how to solve ODE with boundary at infinity, which is similar both in form and in having a single solution (for a given initial value for y[0]) that does not diverge to infinity. For that approach, we have to set up the problem as a boundary-value problem on $[0, \infty)$ instead of an initial-value problem. One can see from the equation that $y = 0$ is a solution, and that for some $x > 1$, if $y$, $y'$ are both positive or both negative, the solution diverges. Thus if a solution crosses the $x$-axis, it will diverge. And if $y$ is positive (negative) and reaches a local minimum (maximum resp.), then $y$ will diverge. Further $y$ cannot approach a finite nonzero value because $|y''|$ will be bounded away from $0$. One might surmise that there is a solution that starts at $y = 1$ and approaches $y = 0$. (It's harder to know that this will have the initial condition $y' = 0$, but I suppose one could check that afterwards. If it doesn't satisfy it, then I suppose the IVP actually has a diverging solution and the straightforward NDSolve solution should be used.)

First we transform the differential equation with the substitutions $x = \tan t$ and $u(t) = y(\tan t)$ to transform $[0,\infty)$ to $[0,\pi/2]$. (You might want to evaluate the NestList separately if you cannot see what it does.)

eqn = (-(x[t]^2 - 1) y[x[t]]
+ D[D[y[x[t]], t]/D[x[t], t], t]/D[x[t], t]) x'[t]^3 == 0 /.
x -> Tan;
bvp = {eqn, u[0] == 1, u[Pi/2] == 0} /.
Solve[NestList[D[#, t] &, u[t] == y[Tan[t]], 2],
{y''[Tan[t]], y'[Tan[t]], y[Tan[t]]}];

femsol = y -> (u[ArcTan[#]] &) /.
First@NDSolve[bvp, u, {t, 0, Pi/2},
Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}]

Plot[y[x] /. femsol // Evaluate, {x, 0, 7}]


Comparison with the exact solution:

Plot[Exp[-x^2/2] - (y[x] /. femsol) // Evaluate, {x, 0, 7}, PlotRange -> All]


• The "DifferenceOrder" seems not to be helpful here, at least when it's together with "ImplicitRungeKutta" or "ExplicitRungeKutta" – xzczd Apr 27 '15 at 13:19
• @xzczd I found it modestly helpful, in that increasing it made the divergence occur a little bit later (but not much later). "DifferenceOrder" -> 9 with "ExplicitRungeKutta" diverged around x == 8, as I recall. Since it was so modest, I didn't include an example. – Michael E2 Apr 27 '15 at 13:43
• You may want to take part in the discussion here :D, currently the transformation can be done by Assuming[0 < t < Pi/2, dChange[{y''[x] == (x^2 - 1) y[x], y[0] == 1, y[Infinity] == 0}, ArcTan@x == t, x, t, y[x]]]. – xzczd Apr 28 '15 at 5:29
• @xzczd Thanks. I'm sure what to add to the discussion, but thanks for linking to that useful answer! – Michael E2 Apr 28 '15 at 11:10