How to know if the solution given by Runge Kutta is good?

Question

We're trying to solve a nonlinear ODE problem using the Runge-Kutta method (4th order). How to see how accurate the solution is?

Since we don't know the exact solution, there is nothing to compare the numerical solution to. We thought about decreasing the step size and comapring solutions for varying step sizes, but I don't think that makes sense since different step sizes give vectors of different lengths as solutions.

Basically, we solved a system of nonlinear ODEs using Runge-Kutta, and don't know how good the solution is. Any suggestions?

Answer 1

With the given data use "polyfit" and "polyder" to fit a polynomial into the solution and check to see if your data satisfy the differential equation well enough.

If the error is within your tolerance then the numerical solution is acceptable.

Answer 2

You can compute a second numerical solution with half the step size or double the step size. Then you will have solution points with differing accuracy at the same time points and can compute their difference which will fully consist of the method error.

To be more precise, apply Richardson extrapolation. That is, assume the error for the numerical solution with step-size $h$ is of the form $$ E(t,h)=y_h(t)-y(t)=C(t)h^p+O(h^{p+1}). $$ Then comparing to the solution with twice the step-size gives for the difference $$ y_{2h}(t)-y_h(t)=E(t,2h)-E(t,h)=C(t)(2^p-1)h^p+O(h^{p+1}). $$ This allows to compute the leading error term and thus a refined error estimate as $$ E_R(t,h)=\frac{y_{2h}(t)-y_h(t)}{(2^p-1)}+O(h^{p+1}) $$ as estimate of the current error. Similar applies for using the solution with half the step-size.

• An example for a third order method with error estimation on the chaotic Lorenz system