Extrapolação de Richardson: melhorado a estimativa de uma derivada numérica

$f(x) = e^{x/2}\cos(2\pi x)$

Existem basicamente três formas de melhorar a estimativa da derivada numérica de uma função: diminuir o valor do passo $h$, aumentar o número de pontos no cálculo da derivada ou usar a Extrapolação e Richardson. O primeira forma é (quase) trivial, basta reduzir arbitrariamente o valor do passo $h$, mas um passo muito pequeno aumenta o erro total (ver aqui). O segundo método é ilustrado a seguir (derivada centrada):
$$ f'(x) = \frac{f(x+h) - f(x-h)}{2h} + O(h^2)$$
$$ f'(x) = \frac{-f(x+2h) +8 f(x+h)-8f(x-h)+f(x-2h)}{12h} + O(h^4)$$
Já a Extrapolação de Richardson faz duas estimativas da derivada, uma ($D_1$) com passo $h$ e outra ($D_2$) com o passo reduzido $h/2$, com essas duas estimativas, a nova estimativa mais acurada é dada por:
$$D = \frac{4D_2 - D_1}{3}$$
Vejamos um exemplo. Vamos considerar a função $$f(x) = e^{x/2}\cos(2\pi x)$$ e calcular sua derivada em $x_0 = 1/4$. Os resultados, para um passo $h = 0,1$, são:
D1,                 D2,                 D,                 Derivada teórica:
 -6.6688068,  -7.0054312,    -7.1176393,  -7.1197817

Esse passo, mesmo sendo relativamente grande, gerou um resultado satisfatório com um erro relativo de cerca de 0,03\%. O código Scilab é:

---

//////////////
function f =  fb(x)
    f = exp(x/2).*cos(2*%pi*x);
endfunction

//////////////
function f =  fbd(x)
    f = x.*exp(x/2).*cos(2*%pi*x)/2 - exp(x/2)*sin(2*%pi*x)*2*%pi;
endfunction

x = 1/4;
h = 0.1;
d1 = (fb(x+h)-fb(x-h))/(2*h);
h2 = h/2;
d2 = (fb(x+h2)-fb(x-h2))/(2*h2);
d = (4*d2-d1)/3
disp([d1, d2, d, fbd(1/4)]);