Numerical solution of IVP for linear ODE with variable coefficient blows up

Question

Cross posted in Mathematica.SE, I'll try to rephrase it in a more general way here.

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 and the solution is:

$$y(x)=e^{-\frac{x^2}{2}}$$

But when I tried to solve it with the classical Runge–Kutta method, the numerical solution blows up:

The step size used here is 0.001.

Why did this happen? If I've chosen an improper method, what's a suitable one?

Answer 1

This is a second order ODE so there are two solutions. The second solution is unstable in the sense that it goes either to plus or minus infinity as $x$ becomes large. For $x$ near zero the solution is obviously stable. Even if this solution is the analytic one you show above, in the numerical integration algorithm, due to round off error, a component of the second solution will creep in. Eventually, this component will grow and dominate the solution.

I believe that all standard integration algorithms for solution of initial value ODE will exhibit this behavior.

Answer 2

Your analytical solution assumes that $y>=0$ for $x>0$. However, for initial conditions where $y<0$, $y(x) = -e^{-x^2/2}$. Analyzing the stability of this system in detail is likely off-topic, and would be better suited for a mathematics site.

The solution asymptotically converges to $0$ exponentially. Numerically, most integration methods will have difficulty with this. Numerical noise around small values can cause a solver to take a step that is too large, causing $y(x)$ to go negative. One way to get around this is to use a solver that is able to intelligently constrain solution components to be nonnegative. Doing this correctly isn't trivial – see Shampine, et al. 2005. Matlab's ODE solvers, e.g., ode45, have this feature via the 'NonNegative' output property that can be set with odeset. I don't think Mathematics's NDSolve has such an option. This Mathematica.StackExchange post offers several ways you could try to add it.

Here's a quick example in Matlab:

f = @(x,y)[y(2);(x^2-1)*y(1)];
opts = odeset('NonNegative',1); % Keep y(1) >= 0
[x,y] = ode45(f,[0 10],[1 0],opts);

Plotting, you get something like this:

Answer 3

The complete solution is $$ y(x)=C·e^{-x^2/2}+D·e^{-x^2/2}·\int_0^x e^{s^2}ds $$ The first is an even function and bounded, the second an odd function, has a not further reducible integral and is exploding like $\sim x·e^{x^2/2}$.