Why is the numerical solution of this equation unstable? Is this equation stiff?
Solution 1:
Yes, this initial value problem is made to do such things. At the request of Robert Lewis, here is a picture: red curve is the analytic solution, blue curve is the textbook implementation of RK4 with the parameters from your post.
(The blue curve is pretty much vertical from there.) I don't think that this equation qualifies as stiff (the term I never really understood), it's just hypersensitive to initial condition. Indeed, the exact solution with initial condition $y(0)=1+\delta$ is $$y(t) = t\sin t + \cos t + \delta\, e^{t^2/2} \tag{1}$$ where the last term spells trouble. As soon as we step off the exact trajectory the tiniest bit, the $e^{t^2/2}$ term kicks in and dooms the solution, even if the subsequent computations have no error at all.
This is not really a fault of the method. Say, the initial value $y(0)$ is known up to $2^{-53}$ (in double precision). Then what computer really knows about $\delta$ in (1) is that $|\delta|\le 2^{-53}$. But $2^{-53} e^{t^2/2} > 500,000$ when $t=10$. This is for the analytic solution with machine-precision initial value. So, the problem is not (only) in the RK4 method, we don't have enough precision at the beginning to get a decent result at the end. Carrying out the computations in quadruple precision should help.
In a way, the trigonometric stuff here is mostly a decoration. The simpler equation $y'=ty$ already has the same feature. With initial value $0$ the solution is $y\equiv 0$. But with a nonzero initial value $C$ it's $y= C \exp (t^2/2)$, which is huge.
By the way, I tried stiff ODE solvers like ode15s in Matlab; that did not help at all.
Just as an aside: consider what happens if we reverse the order of time axis. That is, let $t = 10-s$, so the unknown function is $z(s) = y(10-t)$. Then solve $$z' = - (10-s) (z-(10-s)\sin(10-s))$$ with some initial condition $z(0)=C$. Then the numerical method works fine; you don't need very small steps. For example, below I solved with $C = -6.1$ using mere $40$ steps of RK4 and plotted the solution in reverse:
In a way, this is cheating: how did I know that $C=-6.1$ is the right value to use? (I took it from the endpoint value of the exact solution.)
Solution 2:
I know this question is old, but I thought I would post some extra info around this.
One reason your numerical solution to the Ordinary Differential Equation (ODE) is exploding, kind of regardless of the numerical method used, is because of the ODE itself. If you go and study some theory on numerical integration of ODEs, you will find a common model problem used to study numerical methods is:
\begin{align} y' = \lambda y \end{align}
Where $\lambda$ can be a complex number. The general solution to this ODE is obviously
\begin{align} y(t) = A e^{\lambda t} \end{align}
What you will find is that numerical methods basically always diverge from the real solution when $Re(\lambda) > 0$.
If you look back at your ODE, the $ty$ term is analogous to the $\lambda y$ term in the model ODE. One interesting thing to note is since $t$ is positive during the lifetime of this simulation, this means, analogously, your ODE has a $Re(\lambda) > 0$, meaning you should expect the solution to diverge and thus blow up basically no matter what numerical method you use.