How close observed order of accuracy should be to theoretical order of accuracy?

Question

I try to write tests to the implementations of numerical methods. The best way to do this is to study the observed order of accuracy and check that it matches with the theoretical order of accuracy.

However, I noticed that the observed order of accuracy lies in a band $$ [p_{\mathrm{theor}} - \Delta, p_{\mathrm{theor}} + \Delta], $$ where $p_{\mathrm{theor}}$ is the theoretical order of accuracy and $\Delta$ is some absolute error that must be chosen for any particular method. For some of the methods I implement I can take $\Delta=0.05$, for other $\Delta=0.1$ and so forth.

So, my question is: how should one choose $\Delta$ to be sure that there are no mistakes in the implementation?

Answer 1

From practical experience, I find plotting the results of discretization refinement studies best to figure out whether or not the observed order of accuracy is "close enough". However, writing an automated test of a method's order of accuracy requires some sort of regression to calculate the observed order of accuracy from data, and then defining an appropriate $\Delta$ -- a comparison tolerance -- such that if the observed order of accuracy is within $\Delta$ of the theoretical order of accuracy, the test "passes".

I don't know of any theory that suggests an appropriate $\Delta$ or an appropriate way to bracket the asymptotic order of convergence region a priori. I have always written these kinds of tests by reverse-engineering them from plots, so I typically choose $\Delta = .2$ or $\Delta = .3$, because I see the order of accuracy deviate by at most that amount.

In certain cases (e.g., discontinuous initial conditions, or stiff problems), theory can predict how the order of accuracy of a method will degrade; if this additional theory isn't taken into account, then you might need to choose a larger $\Delta$ to compensate. For instance, in Section 8.7 of "Finite-Volume Methods for Hyperbolic Problems", Randy LeVeque discusses how convergence theory for a first-order upwind method assumes smooth solutions. If instead, a weak solution is discontinuous, the order of accuracy in space (in the $L^{1}$ norm) is one-half instead of one. So in a case like that, if you assume that the true order should be one, you might need to make $\Delta$ as large as 0.6, but in some cases, the degradation of the error is known, and you can adjust the theoretical order of accuracy accordingly.