Jeśli w złożonej kwadraturze trapezów
\(\displaystyle{ T_{n} = \frac{h}{2}[ f(a) +f(b)] +h\sum_{i=1}^{n-1}f(x_{i}) }\)
zastąpimy \(\displaystyle{ n }\) przez potęgę \(\displaystyle{ 2 , \ \ n = 2^{k-1} , \ \ k = 1,2,... }\)
wówczas długość kroku całkowania wyniesie \(\displaystyle{ h = \frac{b-a}{n} = \frac{b-a}{2^{k-1}}. }\)
i otrzymamy kwadraturę Romberga
\(\displaystyle{ R_{k,1} = \frac{b-a}{2^{k-1}} [ f(a) + f(b)] + \frac{b-a}{2^{k-1}} \sum_{i=1}^{2^{k-1}-1} f\left( a + \frac{b-a}{2^{k-1}}i \right), \ \ k=1,2,... \ \ (1)}\)
Stąd
\(\displaystyle{ R_{1,1} = \frac{b-a}{2} [ f(a) + f(b)],}\)
\(\displaystyle{ R_{2,1} = \frac{b-a}{4} [ f(a)+ f(b)]+ \frac{b+a}{2} f\left( a + \frac{b-a}{2}\right), }\)
\(\displaystyle{ R_{3,1} = \frac{b-a}{8}[ f(a) + f(b)] + \frac{b-a}{4}\sum_{i=1}^{3}f\left( a + \frac{b-a}{4}i\right),}\)
itd. ....
Zauważmy, że:
\(\displaystyle{ R_{2,1} = \frac{R_{1,1}}{2} + \frac{b-a}{2} f\left( a + \frac{b-a}{2}\right) }\)
\(\displaystyle{ R_{3,1} = \frac{R_{2,1}}{2}+ \frac{b-a}{4}\sum_{=1}^{2} f\left [a + \frac{b-a}{2}\left(i-\frac{1}{2}\right)\right] }\)
itd. ...
Stosując indukcję wzgkędem \(\displaystyle{ k }\) otrzymujemy równanie rekurencyjne
\(\displaystyle{ R_{k,1} = \frac{R_{k-1,1}}{2} + \frac{b-a}{2^{k-1}} \sum_{i=1}^{2^{k-2}} f\left [a + \frac{b-a}{2^{k-2}}\left(i-\frac{1}{2}\right)\right] \ \ (2) }\)
Ten schemat rekurencyjny może być użyty do obliczenia wartości elementów ciągu \(\displaystyle{ R_{2,1}, R_{3,1}, ..., R_{n,1} }\) pierwszej
kolumny tablicy kwadratur Romberga.
Można pokazać (*) że błędy zaokrągleń w schemacie Romberga spełniają równanie
\(\displaystyle{ C(h^2 ) + D(h^4) + E(h^6) + \ \ ... }\) z \(\displaystyle{ h = \frac{b-a}{2^{n-1}} }\).
W celu ich zmniejszenia stosuje się ekstrapolację Richardsona.
Dla kwadratur \(\displaystyle{ R_{1,1}, R_{2,1}, }\) przekształcając równania:
\(\displaystyle{ I = R_{1,1} + C(h^2)+D(h^4) +... }\)
\(\displaystyle{ I = R_{2,1} + C\left(\frac{h^4}{4}\right) + D \left( \frac{h^4}{16}\right) +... }\)
otrzymujemy:
\(\displaystyle{ R_{2,2} = \frac{4R_{2,1}-R_{1,1}}{3} }\)
Ta ekstrepolacja eliminuje błąd \(\displaystyle{ C(h^2).}\)
Podobnie, jeśli ekstrapolujemy kwadratury \(\displaystyle{ R_{2,2}, R_{3,2}, }\)
dostajemy nowe przybliżenie całki \(\displaystyle{ I }\)
\(\displaystyle{ R_{3,3} = \frac{16R_{3,2}-R_{2,2}}{15}.}\)
Z kolei ta ekstrepolacja eliminuje błąd \(\displaystyle{ D(h^4).}\)
Kontynuując procedurę - otrzymujemy ogólny wzór ekstrapolacyjny :
\(\displaystyle{ R_{i,k} = \frac{4^{k-1} R_{i,k-1} - R_{i-1,k-1}}{4^{k-1} -1} \ \ (3) }\) dla \(\displaystyle{ i=2, ...,n.}\)
Błąd obcięcia dla \(\displaystyle{ R_{k,k} }\) wynosi \(\displaystyle{ O\left(h^{2k}\right)}\) (*)
(*) Andrzej Kiełbasiński, Janusz Chojnacki. Błędy zaokrągleń w algorytmie Romberga. Wydawnictwa Uniwersytetu Warszawskiego 1985.
Miałem przyjemność jako współautor uczestniczyć w opracowaniu tej monografii.
Program w MATLAB
Kod: Zaznacz cały
function romberg(f,a,b,n)
% Obliczenie całki funkcji f w przedziale [a,b] metodą Romberga
fprintf('\n')
disp(' Tablica Romberga ')
disp('________________________________________________________')
disp(' i h Ri1 Ri2 Ri3 ...')
h=b-a;
R(1,1)=h*(feval(f,a)+feval(f,b))/2;
fprintf('%2.0f %8.4f %12.4f\n',1,h,R(1,1));
m=1;
for i=1:n-1
h=h/2;
S=0;
for j=1:m
x=a+h*(2*j-1);
S=S+feval(f,x);
end
R(i+1,1)=R(i,1)/2+h*S;
fprintf('%2.0f %8.4f %12.4f',i+1,h,R(i+1,1));
m=2*m;
for k=1:i
R(i+1,k+1)=R(i+1,k)+(R(i+1,k)-R(i,k))/(4^k-1);
fprintf('%12.4f',R(i+1,k+1));
end
fprintf('\n');
end
Kod: Zaznacz cały
function f=f1(x)
f=7+14*x.^6;
>> romberg('f1',0,1,6)
Tablica Romberga
________________________________________________________
i h Ri1 Ri2 Ri3 ...
1 1.0000 14.0000
2 0.5000 10.6094 9.4792
3 0.2500 9.4285 9.0348 9.0052
4 0.1250 9.1088 9.0023 9.0001 9.0000
5 0.0625 9.0273 9.0001 9.0000 9.0000 9.0000
6 0.0313 9.0068 9.0000 9.0000 9.0000 9.0000 9.0000