RK4 seems broken

[mseri]

RK4 seems broken

[tachukao]

@mseri I just found the bug in the RK4 stepper function, should have been more careful when I wrote the function 🤦.
I also noticed something weird with RK4 a while back, but I never looked into it.
Basically, the problem was with the update:
let dy = M.((k1 + k2 + k3 + k4) /$ 6.)
It should be
let dy = M.((k1 + (2. $* k2) + (2. $* k3) + k4) /$ 6.)
I've made the necessary changes and pushed. But yea I agree, we should add a test suite to test all these stepper functions 👍

[ryanrhymes]

Agree, we can use matlab (or other tools) to generate a set of solutions beforehand, then use owl-ode to check against it.

[tbrk]

One example that we often use in Zélus is:

let f t y yd =
  yd.{0} <- (-. y.{1}); (* cosine *)
  yd.{1} <- y.{0} (* sine *)

let y_init = [| 1.0; 0.0 |]

[ryanrhymes]

can we add a test suite similar to owl so others can recreate the issue?

[mseri]

It's on the pipeline. I am still trying to figure out the best way of testing this. Maybe by taking a set of known ode+solutions and checking that the error is somewhat O(h^(method-order)).

[mseri]

We could even pick some equations with closed form solutions if we want to be more reliable, at that point is really about comparing numerical and analytical solutions.

[mseri]

Indeed, that's a quite good example.

PS. Zélus seems awesome, I did not know about it, thanks for mentioning it!

[mseri]

Thanks! I’ll try to prepare some tests following tbrk suggestion soon

[mseri]

I have a proposal for testing: #6
If you agree with the approach, I am going to merge and release a new version (including also @tachukao fixes)