Kwadratury Romberga

Sigma-ciała i zbiory borelowskie. Miary, miary zewnętrze i miara Lebesgue'a. Funkcje mierzalne. Całka Lebesgue'a. Inne zagadnienia analizy rzeczywistej.
janusz47
Użytkownik
Użytkownik
Posty: 7910
Rejestracja: 18 mar 2009, o 16:24
Płeć: Mężczyzna
Podziękował: 30 razy
Pomógł: 1670 razy

Kwadratury Romberga

Post autor: janusz47 »

Kwadratury Romberga

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

Przykład

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
 
ODPOWIEDZ