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}

enter image description here

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}

enter image description here

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}

enter image description here

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.


1 Answer 1


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];
 Plot[y[x] /. asol /. {{p -> 0}, {p -> 10^(-error)}, {p -> -10^(-error)}} // Evaluate,
  {x, 0, 7}, PlotRange -> 2],
 {error, 1, 10, Appearance -> "Labeled"}]

Mathematica graphics

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}]

Mathematica graphics

Comparison with the exact solution:

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

Mathematica graphics

  • $\begingroup$ The "DifferenceOrder" seems not to be helpful here, at least when it's together with "ImplicitRungeKutta" or "ExplicitRungeKutta" $\endgroup$
    – xzczd
    Commented Apr 27, 2015 at 13:19
  • 1
    $\begingroup$ @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. $\endgroup$
    – Michael E2
    Commented Apr 27, 2015 at 13:43
  • 2
    $\begingroup$ 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]]]. $\endgroup$
    – xzczd
    Commented Apr 28, 2015 at 5:29
  • $\begingroup$ @xzczd Thanks. I'm sure what to add to the discussion, but thanks for linking to that useful answer! $\endgroup$
    – Michael E2
    Commented Apr 28, 2015 at 11:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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