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

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?

• There's obviously some numerical instability at play here. But without knowing what Mathematica is actually doing when you call NDSolve, it's impossible to give an explanation, much less a solution. Is there some diagnostic information you can obtain (e.g., which integrator is used, how are the time steps chosen, are any error estimators used)? Otherwise I would say the question is off-topic here, since it requires in-depth Mathematica knowledge to answer. Apr 27 '15 at 11:22
• @ChristianClason I admit my question wasn't clear enough. NDSolve isn't really the issue here. I've clarified a little, have a look. Apr 27 '15 at 12:12
• This is precisely my point: Your question only contains details about the software, but none about the mathematics. So you can only expect answers related to the former (for which Mathematica.SE is the right place), not the latter (for which this is the right place). Apr 27 '15 at 12:49
• @ChristianClason How about now? Apr 27 '15 at 13:13
• Much better! I'm assuming you mean RK4, right? Apr 27 '15 at 14:12

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.

• You mean general solution for the ODE itself is a linear combination of two linearly independent solutions by saying "This is a second order ODE so there are two solutions", right? :) Apr 28 '15 at 6:23
• Yes. You have one solution that decays. The other explodes and will ultimately dominate in a general numerical solution algorithm. Apr 28 '15 at 12:49

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: • @DavidKetcheson: Maybe if this were a first order system, but it's second order. The plot in the OP's question only shows the (mis)behavior of $y(x)$, but doesn't show $\dot{y}(x)$. Apr 28 '15 at 4:48
• ... additionally, this system is nonautonomous so if the initial conditions were different from those in the question it might be possible to cause $y(x)$ to go to $-\infty$ instead. A different integration scheme might also cause the solution to diverge negatively due to different numerical errors in each state variable. Apr 28 '15 at 5:05
• Sorry, but I also think that your analysis is at least unclear. Bill Greene's answer is in agreement with Michael E2's. What causing trouble here is likely to be the value of $y'(x)$ rather than $y(x)$. Anyway, forcing a change on the value of the function is indeed a solution. (This can be done in Mma by l = 10; nsol = NDSolve[{y''[x] == (x^2 - 1) y[x], y == 1, y' == 0, WhenEvent[y[x] == 0, y'[x] -> 0]}, y, {x, 0, l}]; Plot[y[x] /. nsol, {x, 0, l}, PlotRange -> 1]) +1 Apr 28 '15 at 6:30

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}$.