How to verify the order of DOPRI Runge-Kutta method

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.

• You mean Dormand and Price's paper? I have textbooks that describe the implementation of embedded RK methods and provide the coefficient table needed for their method. – Mr. G May 14 '13 at 3:39
• I mean "P.J. Prince and J.R.Dormand: High order embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics 7(1), 1981". Also, can you access this command to validate some of the Butcher in your code? mathworks.com/help/symbolic/mupad_ref/numeric-butcher.html – Amzoti May 14 '13 at 3:41
• Hmm, according to MatLab, the Butcher table for DOPRI87 is actually a lot more complicated than the one in my textbook and only agrees with it up to 16 digits. It isn't mentioned anywhere in the textbook that the coefficients are actually just decent approximations, but I guess it makes sense that people would just use those values since the actual ones are huge. I'm not sure if that explains the issues with the polynomials though. – Mr. G May 14 '13 at 4:25
• Yeah, I should probably take a look at that paper now, thanks. – Mr. G May 14 '13 at 4:35
• Have you seen the discussion in Hairer/Nørsett/Wanner? In particular, they explain there why they've elected to use a 5th order and a 3rd order embedded RK method along with the Dormand-Prince coefficients for error estimation purposes, instead of the seventh-order method. – J. M. is a poor mathematician May 14 '13 at 6:22

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).
If I remember correctly, in the special case where $f(t,y)$ is a polynomial in $t$ of degree $d$ and independent of $y$, every Runge-Kutta method of order $d+1$ gives the exact solution (which is a polynomial of degree $d+1$). So the Euler method solves $y' = 1$ exactly and the classical fourth-order RK-method solves $y' = t^3$ exactly (all up to round-off error).