Let $T$ denote the target value and let $A_h$ denote our approximation. If there exists an asymptotic error expansion of the form
$$T - A_h = \alpha h^p + \beta h^q + O(h^r), \quad h \rightarrow 0_+,$$
where $\alpha$ and $\beta$ are constants independent of $h$ and the exponents $$0<p<q<r$$ are not necessarily integers, then the following is true:
- We have $$\alpha h^p = \frac{A_h - A_{2h}}{2^p - 1} + O(h^q), \quad h \rightarrow 0_+.$$
- If we define $$F_h = \frac{A_{2h} - A_{4h}}{A_h - A_{2h}},$$ then $$F_h \rightarrow 2^p, \quad h \rightarrow 0_+$$
- Moreover, there is a constant $c$ such that $$\frac{F_h - 2^p}{h^n} \rightarrow c, \quad h \rightarrow 0_+,$$ where $n=q-p$.
Item 1 contains Richardson's error estimate. Item 2 and 3 concerns Richardson's fraction $F_h$. By monitoring the computed values of Richardson's fraction you can determine the value of $p$ and identify the asymptotic range, i.e., the set of values of $h$ where rounding errors are insignificant compared with the discretization error. Item 3 show that you can determine the difference between the order of the primary error term $\alpha h^p$ and the secondary error term $\beta h^q$ by observing the valus of $F_h - 2^p$.
To be fair, these observations are hard to turn into software which automatically detects the values of $p$ and $q$. It is essentially impossible to distinguish between say $p=2$ and $p=2.001$. It is important to realize we don't need the exact value of $p$ to get a good approximation of Richardson error estimate. We need a good approximation of $2^p-1$. The asymptotic range is characterized by $$\log(|F_h - 2^p|) \approx \log|c| + n \log(h).$$
As an example of what is possible for a human we consider the problem of computing the integral $$T = \int_0^1 e^x dx = e^1 - 1.$$ We choose the composite trapezoidal rule as our approximation $A_h$. It is known that $p=2$ and $q=4$ for this rule, but we shall verify this experimentally. We have the following data:
k | Approximation A_h | Fraction F_h | Error estimate E_h | Error | Comparison
0 | 1.859140914230e+00 | 0.00000000 | 0.000000000000e+00 | -1.408590857705e-01 | 0.0000e+00
1 | 1.753931092465e+00 | 0.00000000 | -3.506994058823e-02 | -3.564926400578e-02 | -1.7891e+00
2 | 1.727221904558e+00 | 3.93908726 | -8.903062635770e-03 | -8.940076098471e-03 | -2.3830e+00
3 | 1.720518592164e+00 | 3.98447608 | -2.234437464405e-03 | -2.236763705256e-03 | -2.9830e+00
4 | 1.718841128580e+00 | 3.99610010 | -5.591545281024e-04 | -5.593001209490e-04 | -3.5845e+00
5 | 1.718421660316e+00 | 3.99902383 | -1.398227545558e-04 | -1.398318572816e-04 | -4.1864e+00
6 | 1.718316786850e+00 | 3.99975588 | -3.495782207789e-05 | -3.495839104795e-05 | -4.7885e+00
7 | 1.718290568083e+00 | 3.99993897 | -8.739588871635e-06 | -8.739624433041e-06 | -5.3905e+00
8 | 1.718284013367e+00 | 3.99998474 | -2.184905553001e-06 | -2.184907774039e-06 | -5.9929e+00
9 | 1.718282374686e+00 | 3.99999619 | -5.462269084452e-07 | -5.462270487033e-07 | -6.5904e+00
10 | 1.718281965016e+00 | 3.99999904 | -1.365567600479e-07 | -1.365567685596e-07 | -7.2053e+00
11 | 1.718281862598e+00 | 3.99999982 | -3.413919154778e-08 | -3.413919391626e-08 | -7.1588e+00
12 | 1.718281836994e+00 | 3.99999982 | -8.534798275524e-09 | -8.534799089688e-09 | -7.0205e+00
13 | 1.718281830593e+00 | 3.99999566 | -2.133701881846e-09 | -2.133693444151e-09 | -5.4029e+00
14 | 1.718281828992e+00 | 3.99999681 | -5.334258960469e-10 | -5.334157560100e-10 | -4.7210e+00
15 | 1.718281828592e+00 | 4.00004551 | -1.333549567069e-10 | -1.333508858892e-10 | -4.5153e+00
16 | 1.718281828492e+00 | 4.00022202 | -3.333688880502e-11 | -3.334021947410e-11 | -4.0004e+00
17 | 1.718281828467e+00 | 4.00046186 | -8.333260007968e-12 | -8.340439450194e-12 | -3.0651e+00
18 | 1.718281828461e+00 | 3.98094194 | -2.093288505497e-12 | -2.060573933704e-12 | -1.7992e+00
19 | 1.718281828460e+00 | 4.13782004 | -5.058916248875e-13 | -5.428990590417e-13 | -1.1664e+00
The step sizes used are $h_k = 2^{-k}$. $A_h$ is the approximation, $F_h$ is Richardson's fraction (which is only defined for $k \ge 3$) $E_h$ is Richardson's error estimate (which is only defined for $k\ge 2$). The two last two columns contain the actual error and compares Richardson's error estimate to the actual error. Specifically, the actual error is treated as a target value, Richardson's error estimate as the approximation and the last column contains $\log_{10}$ of the associated relative error.
This figure shows displays the convergence of Richardson's fraction:
From the table it is reasonably clear that if $p$ is an integer, then it has to be $p=2$ because the computed values of $F_h$ (initially) appear to converge toward $4 = 2^2$. The picture reveals that the asymptotic range is $k=2,3, \dots, 10$. The slope of the straight line is $-2$, which shows that $q=p+2=4$. At $k=11$ we have reached the point where rounding errors are comparable with the discretization error and the behavior of the computed values of $F_h$ start to deviate from the behavior of the exact values of $F_h$. At this stage, there is no reason to continue. You will note that $k=10$ is also the point where the accuracy of Richardson error estimate is maximal, the relative error is about $10^{-7.2}$ which is exceedingly good.
