Robust Numerical ODE Solver?

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: • BTW, there's an ExplicitRungeKutta method for NDSolve, it seems to be just what you're trying to implement. Try it to see if it works there. It's likely that you have some bug in your implementation. – Ruslan Jun 6 '14 at 15:00
• I've been using NDSolve for a while and it's great, but I always like seeing what goes on in the "blackbox" so wanted to create my own numerical solver. – InquisitiveInquirer Jun 6 '14 at 15:02
• Your ODE seem not too bad. Can you give some details on the implementation of the solver? What kind step sizes are you using, what of error do you see, in the size of the step length or significantly larger? How do the solutions differ when halving the step size from one to the next? – LutzL Jun 6 '14 at 15:08
• I was using various step sizes between 0.005 and 0.0000005 the errors either didn't match NDSolve's output or "looked wrong". I'll attach some example output in my main post for you to look at. – InquisitiveInquirer Jun 6 '14 at 15:28
• It could be that both your graphs are correct. At the scale you use, the phase portrait looks rather complex. The least you have to do to deal with this is a dynamical step size control, either using embedded methods or two steps of half the step size to gauge the local error. – LutzL Jun 6 '14 at 19:33

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.

• Brilliant, thank you Ruslan, getting the output for $y(t)$ fixed it! That brings me to two questions: How robust is the RK4 solver, and do you know where I can find a good source for creating an implicit RK4 solver, as I hear the implicit solver will be better than an explicit one if you're trying to integrate an n-body simulation. – InquisitiveInquirer Jun 6 '14 at 22:35
• Not sure about "how robust" it is (the question is very broad), but I'd recommend a book from "Numerical Recipes" series. It's expensive, but covers a large set of topics about numerical methods, including implicit solvers and Runge-Kutta methods of various orders. – Ruslan Jun 7 '14 at 8:58
• Sorry, maybe I should rephrase the question. It seems like RK4 is the most universally known type of numerical ODE solver, but it might not necessarily be the best. Are there solvers that are better than RK4 (judging by how quickly NDSolve finds an accurate solution compared to my RK4 solver it seems this is definitely the case)? One solver already mentioned below is the adaptive RK-Fehlberg solver but there could be something even better than this, and I'm guessing there is. – InquisitiveInquirer Jun 13 '14 at 15:42
• There's plenty of numerical algorithms, which could perform better, including RK of higher orders. Many of algorithms are covered in Numeric Recipes. – Ruslan Jun 13 '14 at 17:41

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.

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