How do you compute relative error when the exact solution is unknown?

Question

I'd have a rather complex system of non linear ODEs and with a lot of help I've written an algorithm that solves them. I'd now like to compute the relative error, but I do not have a known solution to compare it with.

One thought was to compare my algorithm to MATLABs ode45 function (I wrote an RK4 method in python). However, I don't know what step size ode45 is using. And, my algorithm is in python, so I'd probably have to print to file and compare data in excel or something.

I have also read about "Approximate Relative Error", where you compare one step in the solution to the last. I don't intuitively understand how this depicts error. It just seems like it measures the smoothness of the estimate rather than any error.

I have also thought of decreasing step size and comparing the solution at a single point with two different step sizes.

I like my first idea best, because it's nice to have a "known" solution. Not sure if that is the most revealing approach, though.

Answer 1

If your derivative function is called odesys and you have a time array tspan that contains the step times of the RK4 method, then you can use the scipy package and use, in the main points, something like

import scipy.integrate as spi

res = spi.solve_ivp(odesys, tspan[[0,-1]], y0, t_eval=tspan, method="RK45", atol=1e-9, rtol=1e-12)
print res.message

diff = y-res.y

---

I hope there is an odesys(t,y) system function, vector y to vector \dot y. If not, the code for any integration method becomes complex and error prone.

To the first question, ode45 does not have one fixed step size, it has a step size controller that adapts the step size to the local behavior of the system. All values in-between are interpolated using a method-adapted polynomial interpolation. You can see these internal steps in matlab ode45 when giving the integration interval only by the end-points and python solve_ivp by not specifying the t_eval parameter. Then the time and value list returned by the method contain those internal points.