3

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?

user490393
  • 91
  • First idea that comes to my mind: If you inject your solution into the ODE and take the norm (integral of the absolute value over your range), maybe it will give you a way to compare your different solutions ? – Jeanba Feb 25 '20 at 08:06

2 Answers2

4

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.

Lutz Lehmann
  • 126,666
  • Thank you. So, when I divide the step by two, it makes sense to compare the points that are shared? So, if I calculate two solutions with steps dt and dt/2, I can use every other point from the finer solution to compare it with the one where the step is dt? – user490393 Feb 25 '20 at 11:02
  • Yes, exactly that. Depending on the language used, what array operations are available, this can be done rather easily, like in matlab or python-numpy. – Lutz Lehmann Feb 25 '20 at 11:47
  • If you use a Taylor Series Method solver, you can try the same trick with integrators of steadily increasing order. You will find that the "clean simulation time" (where the standard and half-step answers agree) increases with integrator order, until it reaches a maximum value, determined entirely by the precision of the floating point library that you are using. I have been implementing this technique for years now, but sadly it is not very well known (it dates back at least to the 1970s). See https://arxiv.org/abs/1111.7149 for a review of the method. – m4r35n357 Aug 11 '23 at 15:39
  • I can (if requested) provide a 60-line version of the Lorenz ODE problem written in GNU bc (which uses arbitrary precision via the GNU MPFR library) which implements the arbitrary-order Taylor Series Method. – m4r35n357 Aug 11 '23 at 15:44
  • 1
    @m4r35n357 : See also the TADIFF part of the FADBAD C++ module, it also is used to Taylor expand the step and provide an error bound i an easy way. // Comparing methods of different orders is, if one wants to remain close to RK4, directly realized by embedded methods, for instance the low-order ones by Fehlberg. – Lutz Lehmann Aug 11 '23 at 15:56
  • 1
    I have written a full "playground" for this in c with a YAD user interface, including "clean numerical simulation", and "chaos scanning" (bifurcation diagrams). I looked at the standard offerings on autodiff.org, but they were all too "heavyweight" for my tastes (you can add TIDES to the list). So I implemented the algorithm and dependencies in a "close to the metal" way, so for example you have to write your ODEs in terms of function calls. My code is tiny as a result, it is also as efficient as it can be, and shockingly accurate (in its MPFR incarnation). Let's just say I am proud of it! – m4r35n357 Aug 11 '23 at 17:10
  • For some reason (i.e. lack of attention) I wrote "dependencies" instead of "recurrences" in my comment above! – m4r35n357 Aug 12 '23 at 08:43
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.

Mohammad Riazi-Kermani
  • 68,728
  • 1
    Is this a standard technique? I've never heard of doing this before. – littleO Feb 25 '20 at 08:14
  • If I understand this, I should use the points that Runge-Kutta determined to fit a polynomial, and then see how these polynomials fit the ODE? How to determine if the polynomial is good? – user490393 Feb 25 '20 at 10:57
  • No, this is wrong. In this situation, the error can grow exponentially over time. A simple example is given by $x'(t) = \lambda x(t)$. If, say, $y'(t) = \lambda y(t) + \epsilon$ and $y(0) = x(0)$, then the error $e(t) = x(t) - y(t)$ satisfies $e'(t) = \lambda e(t) - \epsilon$. It follows that $e(t) = \epsilon(1-e^{\lambda t})/\lambda$. – Carl Christian Feb 25 '20 at 11:16
  • Using the "dense output" interpolation, you actually get a third degree polynomial for each RK4 step that is 4th order accurate (only the computed next point is 5th order accurate). One could use this polynomial to extrapolate the next step and examine the error to the actual RK4 step. – Lutz Lehmann Feb 25 '20 at 11:53