Richardson extrapolation is a well known technique in scientific computing. It is used to compute error estimates when our approximations can be viewed as a function of a parameter $h$. Common applications are numerical differentiation, integration, and the solution of ordinary as well as partial differential equations. Typically, we seek to compute a target value $T$ and we have the ability to compute approximations $A = A_h$ for many or all $h>0$. Now, 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, \quad h > 0,$$ where $$0 < p < q < r, \quad \alpha \not = 0, \quad \beta \not = 0,$$ then Richardson's error estimate is given by $$E_h^{est} = \frac{A_h - A_{2h}}{2^p-1}$$ and $$T - A_h = E_h^{est} + O(h^q), \quad h \rightarrow 0, \quad h >0.$$ This is true in exact arithmetic only. In practice, Richardson's error estimate suffers from subtractive cancellation when $h$ is sufficiently small and $A_h \approx A_{2h}$.
Is there any published literature on the subject of estimating the accuracy of Richardson's error estimate in presence of rounding errors?
Restrictions:
- The typical paper will compare Richardson's error estimate to the actual error for problems where the exact target value is known. I am not interested in such papers. It is the ability to accurately estimate the accuracy of Richardson's error estimate in the general case which is of interest to me.
- In the general case, we do not have a formula for the coefficients, i.e., $\alpha$ and $\beta$, which define the asymptotic error expansion.
- I am not interested in using the method $M_h$ given by $$M_h = A_h + E_h^{est}.$$ It is known that this method satisfies an asymptotic error expansion of the form $$T-M_h = \alpha' h^q + O(h^r), \quad h \rightarrow 0, \quad h>0.$$ Therefore, the error associate with $M_h$ can be estimated as $$ T - M_h \approx \frac{M_h - M_{2h}}{2^q-1},$$ amd this would, in principle, allow us to gauge the accuracy of the original error estimate. However, this derived method is even more vulnerable to subtractive cancellation than the parent $A_h$.
The very practical problem of computing the values of $p$ and $q$ is not an issue. It is straightforward to verify that Richardson's fraction $F_h$ given by $$F_h = \frac{A_{2h} - A_{4h}}{A_h - A_{2h}}$$ satisfies $$F_h - 2^p = \Theta(h^{q-p}), \quad h \rightarrow 0, \quad h > 0.$$ Therefore, the values of $p$ and $q$ can (frequently, but not universally) be obtained by observing the computed values of Richardson's fraction for several different values of $h$.
