Fix RungeKutta23 interpolator

[hgrecco]

RungaKutta23 interpolator does not match the SciPy result. Test hast been temporarily disabled.

numbakit-ode/nbkode/testsuite/test_against_scipy.py

Line 86 in 67c186a

# TODO: RK23 interpolation is not working correctly, the results do no match SciPy

[Yash-10]

Hello @hgrecco I have been looking into this for a comparison with scipy's implementation.

SciPy's implementation uses:

C = np.array([0, 1/2, 3/4])
A = np.array([
    [0, 0, 0],
    [1/2, 0, 0],
    [0, 3/4, 0]
])
B = np.array([2/9, 1/3, 4/9])
E = np.array([5/72, -1/12, -1/9, 1/8])
P = np.array([[1, -4 / 3, 5 / 9],
              [0, 1, -2/3],
              [0, 4/3, -8/9],
              [0, -1, 1]])

whereas `numbakit-ode's implementation uses

C = np.array([0, 1 / 2, 3 / 4], dtype=float)
A = np.array([[0, 0, 0], [1 / 2, 0, 0], [0, 3 / 4, 0]], dtype=float)
B = np.array([2 / 9, 1 / 3, 4 / 9], dtype=float)
B2 = np.array([7 / 24, 1 / 4, 1 / 3, 1 / 8], dtype=float)
P = np.array([[1, -4 / 3, 5 / 9], [0, 1, -2 / 3], [0, 4 / 3, -8 / 9], [0, -1, 1]])

Other implementations such as RK45 and others seem to use the same values but the use of B2 instead of E (with different values) made me wonder whether it is a cause of this issue.

Moreover, I tried to check it locally whether removing B2 and replacing it with E of scipy:

E = np.array([5/72, -1/12, -1/9, 1/8])

solve the error or not. Unfortunately, it still gives an error after running pytest but the difference between nb_y and scipy_sol.y.T seems to have changed.

Could you help me in the above? I am eager to learn more about this implementation.

Thanks!

[maurosilber]

The AdaptiveRungeKutta class (from which RungeKutta23 derives) uses E, which is calculated from B and B2 coefficients here. There seems to be just a difference in sign with respect to Scipy's. But it shouldn't matter, as it is later squared here.

[Yash-10]

Thank you @maurosilber for your reply and guidance. It seems slightly tricky to figure out the cause of mismatch with the scipy's implementation but I would take a look at it.

Thanks!