Metodo di estrapolazione di Romberg

Se la funzione $f$ è molto buona si ha che il metodo dei trapezi composito con ampiezza degli intervalli $h \:I_{1\:n}(f) \equiv T_h(f)$ ammette un sviluppo del tipo $\tau _0 + \tau _1 h^2 + \ldots \tau _r h^{2r} \ldots$ con i vari $\{ \tau _i \}$ indipendenti da $h$ e se $h=0$ si ha che $T_0(f) = I(f)$. Poniamo $x=h^2$ e $g(x)= T_h(f)$. Non possiamo calcolare direttamente $g(0)=I(f)$ e non conosciamo la funzione $g(x)$, conosciamo però il valore che questa assume in $k+1$ punti distinti ed è quindi possibile ricavare il valore del polinomio $p_k(x)$ interpolante $g(x)$ nel punto $x=0$ ed assumere questo valore come approssimazione di $g(0)$ e quindi dell'integrale. Il punto $x=0$ non appartiene all'intervallo in cui si interpola la $g(x)$ e per questo si parla di estrapolazione. Scegliendo i vari $h_i$ decrescenti secondo la formula $h_i = \frac{b-a}{2^{i-1}}$, $x_i = h_i ^2$ si ottiene il metodo di Romberg, e le varie $x_i$ sono vicine al punto 0 più delle $h_i$. Dovendo calcolare il valore del polinomio interpolante in un punto si applica lo schema di Neville opportunamente modificato tenendo conto che deve essere calcolato in $x=0$ e la formula diventa

$$T_{i,j}= T_{i, j-1} + \frac{T_{i,j-1} - T_{i-1,j-1}}{4^j -1}$$

e $T_{i,0}=I_{1\; h_i}(f)$ cioè il valore restituito dal metodo dei trapezi composito con numero di punti via via crescente.

int=romberg(f,a,b,k) - Calcola $\int _a ^ b f(x) dx$ secondo il metodo di Romberg con k punti di interpolazione

%ROMBERG
% function int=romberg(f,a,b,k)
% Calcola l'integrale di f su [a,b] con il metodo di
% estrapolazione di Romberg.
% k e' il numero di valutazioni del metodo dei trapezi che
% devono essere fatte per inizializzare il metodo.
function int=romberg(f,a,b,k)
t=zeros(k,1);
quattro=t;
n=1;
for i=1:k
    t(i)=trapezicomp(f,a,b,n);
    n=2*n;
end

% Calcola le potenze di quattro in un vettore.
q=4;
for i=1:k
    quattro(i)=q;
    q=q*4;
end

for i=1:k
    for j=k:-1:i+1
	t(j)=t(j)+(t(j)-t(j-1))/(quattro(j)-1);
    end
end
int=t(k);

>> romberg('sin',0,pi,5)

ans =

    1.9747

>> romberg('sin',0,pi,12)

ans =

    2.0000

>> romberg('Runge',0,pi,12)

ans =

   1.26248605700994

>> atan(pi)                

ans =

   1.26262725567891