I've written code in Fortran based on the RK8(7)-13 method by Dormand and Prince to solve the system $\mathbf{y}'=\mathbf{f}(t,\mathbf{y})$. The method is Runge-Kutta 8$^\text{th}$ order with an embedded 7$^\text{th}$ order method; I want to verify that the error on the output agrees with this. For a method of order $p$ and timestep $h$, we expect the local truncation error to be $\mathcal{O}(h^{p+1})$. Therefore, if I advance the program by one time step and calculate $\text{error} = |y_{\text{RK}}-y_{\text{exact}}|$ (say I know $y_{\text{exact}}$) and plot $\log(\text{error})$ vs. $\log(h)$ for several values of $h$, the slope of the line should be 9.
Here's the problem: it isn't. I've tried simple polynomials like $f(t,y) = 1$ and $f(t,y) = 2t$, but the slope of the line always seems to be between 1 and 2; this is also the case for $f(t,y) = y$. However, for $y'' + y = 0$ (which translates into the coupled system $\mathbf{f}(t,y_{1}, y_{2}) = (y_{2}, -y_{1})$), I found that the slope is almost exact 9 for larger values of $h$ – for small values the slope goes back to ~1. Furthermore, the magnitude of the error is $~10^{-20}$ for this small $h$. I was wondering if anyone had any ideas as to what the issue could be.
My ideas:
- I've been misinformed as to what order means, and the method I've described is not representative of the order of the RK scheme. This is entirely possible, but the fact that one of the slopes was 9 seems too coincidental.
- Something to do with precision. This could explain why the small $h$ values have unreliable errors since I'm "only" carrying 30 digits of precision, but this wouldn't explain the disagreement for the polynomial $f$.
- Something specific to RK methods that I'm unaware of. I'm hoping it's this.
- Programming error. I'm hoping it isn't this.
- Mistake entering in the Butcher table. I'm really hoping it isn't this.
Side note: I've also tried plotting the global error as a function of $h$, but the results were even stranger, so I think the problem may be more fundamental.
Best Answer
You're just hitting machine precision. Verify the slope of the log for error values larger than say $10^{-12}$ as for smaller errors, you'll see rounding errors due to the double precision arithmetic of your computer, and to the finite number of digits of the coefficients of the butcher table you have (quite often in textbooks, you don't even have the first 16 digits).