Metoda Romberg#

Pentru a putea înțelege mai bine cum funcționează metoda Romberg, este necesar un detur prin extrapolarea Richardson și aplicabilitatea sa în universul numeric.

De aceea, insist să (re)parcurgeți mai întâi acel capitol.

Intuiția metodei#

Metoda trapezelor compuse propune o soluție pentru aproximarea unei integrale folosind o sumă de arii de trapeze ce obține un rezultat cu o eroare de ordinul \(O(h^2)\), unde \(h\in\mathbb{R}_+\).

Prin extrapolarea Richardson, se poate trece de la acest ordin al erorii către orice altă eroare de ordin \(O(h^{2n+2})\), pentru \(n\in\mathbb{N}^*\).

Metoda Romberg va folosi inițial un singur trapez, apoi două, apoi patru și tot așa, continuând să dubleze numărul de trapeze utilizat în timp ce va îmbunătăți acuratețea soluției prin îmbinarea tuturor acestor rezultatelor într-unul singur.

Metoda Romberg#

Fie \(\Psi_i^j\in\mathbb{R}\), \(\forall i\in\overline{0,n}\), \(j\in\overline{0,i}\) și \(n\in\mathbb{N}^*\), valori intermediare ce vor fi folosite pentru a aproxima integrala \(\int_a^bf(x)\,dx\), unde \(f:[a,b]\subset\mathbb{R}\rightarrow\mathbb{R}\) este o funcție aleasă astfel încât \(f\in C^2\).

Dacă prin \(h_i\in\mathbb{R}\), \(\forall i\in\overline{1,n}\), se înțelege \(h_i=\dfrac{b-a}{2^i}\), atunci metoda Romberg se caracterizează prin următoarele valori:

• Prin excepție, dacă \(i=j=0\), se obține:

  \[ \boxed{\Psi_0^0=h_1\big(f(a)+f(b)\big)} \]

  Această valoare reprezintă o aproximare a integralei folosind metoda trapezului.

• Dacă \(j=0\), iar \(i\gt 0\), se respectă formula:

  \[ \boxed{\Psi_i^0=\frac{1}{2}\Psi_{i-1}^0+h_i\sum_{k=1}^{2^{i-1}}f\big(a+(2k-1)h_i\big)}\,,\,\,\forall i\geq1 \]

  Rezultatul astfel obținut este de fapt regula trapezelor compuse utilizând \(2^i\) trapeze.

  Demonstrație

  Folosind regula trapezelor compuse, obținem:

  \[ \int_a^bf(x)\,dx\approx\frac{h_i}{2}\left(f(a)+2\sum_{k=1}^{2^i-1}f(x+kh_i)+f(b)\right) \]

  Printr-o simplă rescriere, separăm termenii sumei:

  \[\begin{split} \int_a^bf(x)\,dx&\approx\frac{h_i}{2}\left[\left(f(a)+2\sum_{k=1}^{2^{i-1}-1}f(x+2kh_i)+f(b)\right)+2\sum_{k=1}^{2^{i-1}}f(x+(2k-1)h_i)\right]\\ &= h_{i-1}\left(f(a)+2\sum_{k=1}^{2^{i-1}-i}f(x+kh_{i-1})+f(b)\right)+h_i\sum_{k=1}^{2^{i-1}}f\big(a+(2k-1)h_i\big) \end{split}\]

  Se observă așadar recurența descrisă anterior.

• Într-un final, recurența ce apare în situația \(i\gt 0\), \(j\gt 0\) este:

  \[ \boxed{\Psi_i^j=\frac{1}{4^j-1}\left(4^j\cdot\Psi_i^{j-1}-\Psi_{i-1}^{j-1}\right)}\,,\,\,\forall i\geq j\geq 1 \]

  Aceasta apare prin utilizarea extrapolării Richardson.

Astfel, fiecare valoare de forma \(\Psi_i^j\) va fi o aproximare a integralei \(\int_a^bf(x)\,dx\) ce va avea un ordin de eroare egal cu \(O(h_i^{2j+2})\), \(\forall i\in\overline{0,n}\), \(j\in\overline{0,i}\).

Exemple#

Să considerăm funcția \(f:[-2,1.5]\rightarrow\mathbb{R}\):

Matematic, știm că \(\displaystyle\int_{-2}^{1.5}f(x)\,dx=\left[\dfrac{x^5}{5}+\dfrac{x^4}{4}-x^3+6x\right]_{-2}^{1.5}=14.809375\). Folosind \(n=2\), se obțin următoarele aproximări, trunchiate la \(6\) zecimale:

Se observă cum rezultatul a convers foarte repede către valoarea corectă. Amintim că singurele valori ce au necesitat evaluarea funcției \(f\) sunt cele de pe coloana \(0\).

Să considerăm și o altă funcție integrabilă și continuă, \(g:[0,1]\rightarrow\mathbb{R}\):

Se poate calcula matematic integrala următoare, poate nu la fel de ușor, cu rezultatul:

Vom folosi \(n=4\) și vom aproxima folosind metoda Romberg - rezultatele trunchiate sunt:

Așadar, în doar câțiva pași, s-a obținut o eroare foarte mică de calcul.

Cod ilustrativ#

Funcția romberg prezintă, în Matlab, o implementare posibilă a metodei Romberg. Drept parametri, primește:

• a - capătul stâng al integralei de calculat;

• b - capătul drept al integralei de calculat;

• n - numărul de pași utilizați;

• f - funcția ce trebuie integrată.

În această situație, funcția va returna I, valoarea integralei căutate, adică \(I\approx\int_a^bf(x)\,dx\).

function I = romberg(a, b, n, f)
    n += 1; # pentru simplitate
    psi = zeros(n, n);
    h = zeros(n, 1);

    h(1) = (b - a) / 2;
    for i = [2 : n]
        h(i) = h(i-1) / 2;
    end

    psi(1,1) = h(1) * (f(a) + f(b));
    for i = [2 : n]
        S = 0;
        for k = [1 : 2^(i-2)]
            S += f(a + (2*k-1) * h(i-1));
        end
        psi(i,1) = psi(i-1,1) / 2 + h(i-1) * S;
    end

    for j = [2 : n]
        for i = [j : n]
            psi(i,j) = 4^(j-1) * psi(i,j-1) - psi(i-1,j-1);
            psi(i,j) /= 4^(j-1) - 1;
        end
    end

    I = psi(n, n);
end

Licență#

The book "Metode Numerice", written by Valentin-Ioan Vintilă, is licensed under CC BY-NC-SA 4.0