Robust Numerical ODE Solver?

Question

I made a little explicit Runge-Kutta 4th order solver a few days ago, but when testing it against various 1st and 2nd order ODEs chosen at random (for example $d^{2}y/dt^{2} = -y \sin(y)$, $d^{2}y/dt^{2} = -yt$ or $d^{2}y/dt^{2} = -y + t^{2}$) it seemed like most were stiff ODEs (unless the algorithm I'm using is incorrect), by comparing my output to that of Mathematica's NDSolve, and hence rendered my RK4 solver to be fairly useless. As such, I've decided to try and find a numerical solver that I can create, that is robust and can solve stiff and non-stiff ODEs. Does such an algorithm exist, or is it a case of the more robust a solver the more abstruse its algorithm becomes. Even better, is there such a thing as a universal solver that is able to solve any ODE you throw at it?

EDIT: Here's an example of my RK4 solver output for $d^{2}y/dt^{2} = -y \sin(y)$ using a step size of $h=0.005$:

And here's what I get from NDSolve:

Answer 1

I've been able to reproduce both your pictures with NDSolve. The second, smoothly-looking one is the solution of $y''(t)=-y(t)\sin(y(t))$ with initial conditions $y(0)=0$ and $y'(0)=50$.

I get the first one if I plot the derivative of the solution:

So, looks like you're taking wrong output from your correct solver. As it's a Runge-Kutta method, you're most likely splitting the equation into system of two equations, one for $y'(t)$ and another for $y(t)$. You're taking the former as the solution to the original equation, while you have to take the latter.

Answer 2

I tried to make an second order ODE solver in ocaml:

let secondode n=
    let rec iode i a b c=
        if i=n then a
        else iode (i+1) (a+.0.001*.b) (b+.0.001*.c) (-.a*.(sin a))
    in iode 0 0. 1.75 0.
;;

secondode 200000;;

and the result is $\textrm{float} = -75208.630048495936$, which significantly differ from the real value $\approx -1.75$. So I assume my naive iterative algorithm is not stable. Here I deliberately choose $0.001$ as it is not too small and not too large. If I choose accuracy to be $10^{-5}$, then the result is $\textrm{float} = -2.050306346414887$, which still differs a lot. And for accuracy at $10^{-6}$, my result is about $-1.80$.

Therefore I think the Runge-Kutta method is largely stable (in comparison with my naive method). But I do not really know how Mathematica does that....

Answer 3

Yes, such algorithms exist.

Implicit multistep methods have been developed deliberately to solve stiff ODE problems. They might require an explicit method to initialize, but this is not an issue because you can scale your step size of your explicit solver in the stiff region rather unobtrusively.

In the non-stiff scales of a problem, the implicit multistep methods are computationally somewhat more expensive but ultimately will yield as an accurate a solution as something like an RK4 method.

Alternatively, you could just use an adaptive solver such as Runge-Kutta-Fehlberg or a Dormand-Prince 4(5) pair.