A generalized Richardson extrapolation routine that can perform multiple extrapolation steps, can accommodate arbitrary leading error terms, can handle any uniform refinement scheme (e.g. h/2, h/3, h/5, etc), and can report a full extrapolation table including comparing against a known good solution. Able to handle vector inputs. Compatible both with MATLAB and Octave. Includes Octave-style unit tests.
A contrived usage example:
Suppose your sequence of approximations was [1, 2, 3, 4] where each successive value refined the grid by a factor of 3. Suppose the three leading error terms in the approximation were O(h^2), O(h^4), and O(h^5) and that you knew the exact value was 4.14. Then you'd invoke
octave:78> [result, normtable] = richardson([1,2,3,4],[2,4,5], 3, 4.14)
result = 4.14163223140496
3.13999999999999968 NaN NaN NaN
2.13999999999999968 2.01499999999999968 NaN NaN
1.13999999999999968 1.01499999999999968 1.00249999999999950 NaN
0.13999999999999968 0.01499999999999968 0.00249999999999950 0.00163223140495905
which shows you the exact result and a full extrapolation table showing the l2 error at each substep.
For more normal examples with h/2 refinement, more standard leading error terms (e.g. starting at h^2 and increasing by h or h^2), and unknown exact solutions the function has much default magic to simplify input parameters.
Fixed a typo or two in the text. No code change.