Cuadraturi adaptive#

În general, algoritmii explicați anterior utilizează metode de tip Newton-Cotes, adică se folosesc evaluări echidistante ale funcției ce trebuie integrată.

Există însă funcții ale căror integrale nu pot fi calculate în această manieră - spre exemplu, funcția \(\cos(x)\) va trece prin axa \(Ox\) la fiecare \(\pi\) pași, iar dacă se utilizează puncte aflate la o distanță \(\pi\) unul de celălalt, atunci aproximarea converge către \(0\), ceea ce nu este mereu adevărat.

Din acest motiv, există o îmbunătățire a algoritmilor precedenți, și anume utilizarea cuadraturilor adaptive.

Descoperirea algoritmului#

Fie \(f:[a,b]\subset\mathbb{R}\rightarrow\mathbb{R}\) o funcție integrabilă și continuă pe întreg intervalul de definiție.

Se pornește apoi de la calculul integralei definite de intervalul \([a,b]\), cu pasul \(h=\dfrac{b-a}{2}\), utilizând orice metodă studiată (vom nota generic \(\text{approx}\) o funcție ce primește ca parametru funcția ce trebuie integrată, capetele de integrare și pasul) - notăm această aproximare cu \(I\):

Putem înțelege prin funcția \(\text{approx}\), de exemplu, metoda Simpson 1/3 - se poate utiliza orice metodă studiată (inclusiv o îmbunătățire a lui Simpson 1/3 prin extrapolare Richardson, asemănător metodei Romberg), însă, pentru simplitate, o vom folosi pe aceasta. Amintim că formula de aproximare a metodei Simpson 1/3 este:

Se împarte apoi intervalul \([a,b]\) în două subintervale de mărime egală, anume \(\left[a,\dfrac{a+b}{2}\right]\) și \(\left[\dfrac{a+b}{2},b\right]\), pe care se calculează din nou integralele definite ale funcției \(f\), de data aceasta însă cu pasul \(h'=\dfrac{h}{2}=\dfrac{b-a}{4}\) - vom nota cele două aproximări cu \(I_1\) și \(I_2\):

Utilizând metoda Simpson 1/3 pentru aproximare, se obține:

Întrucât eroarea introdusă de metoda Simpson este dată de formula \(-\dfrac{v-u}{180}h^4\cdot f^{(4)}(\xi)\) pentru un interval oarecare \([u,v]\subset\mathbb{R}\) și un pas \(h\), concluzionăm că sunt valabile următoarele două egalități:

• Pentru intervalul complet, adică \([a,b]\), integrala va fi:

  \[\begin{split} \int_a^bf(x)\,dx&=\text{approx}\left(f,a,b,h\right)-\frac{b-a}{180}h^4\cdot f^{(4)}(\xi)\\ &=I-\frac{h^5}{90}\cdot f^{(4)}(\xi) \end{split}\]

• Pentru subintervale, adică \(\left[a,\dfrac{a+b}{2}\right]\) și \(\left[\dfrac{a+b}{2},b\right]\), se obține:

  \[\begin{split} \int_a^bf(x)\,dx&=\text{approx}\left(f,a,\frac{a+b}{2},\frac{h}{2}\right)+\text{approx}\left(f,\frac{a+b}{2},b,\frac{h}{2}\right)\\&\,\,\,\,\,\,-\frac{b-a}{180}\left(\frac{h}{2}\right)^4\cdot f^{(4)}(\xi)\\ &=I_1+I_2-\frac{1}{16}\cdot\frac{h^5}{90}\cdot f^{(4)}(\xi) \end{split}\]

Menționăm că valorile \(\varepsilon\) nu sunt neapărat egale, însă, pentru această demonstrație, le vom considera aproximativ egale.

Astfel, cele două formule prezentate sunt și ele aproximativ egale:

Prin simpla jonglare cu termenii ecuației anterioare, obținem următoarea formulă:

Așadar, modulul erorii obținute prin integrarea folosind două subintervale poate fi calculat întorcându-ne la ecuația ce descrie aproximarea pentru subintervale:

Metoda cuadraturilor adaptive#

Utilizând formula anterioară pentru modulul erorii, putem găsi un algoritm ce calculează aproximări din ce în ce mai bune pentru integrala căutată - astfel, obținem următoarea metodă:

1. Se pornește de la un interval \([a,b]\) ce se înjumătățește, obținându-se intervalele \(\left[a,\dfrac{a+b}{2}\right]\) și \(\left[\dfrac{a+b}{2},b\right]\).

   Se calculează apoi aproximările \(I\), \(I_1\) și \(I_2\), conform formulelor:

   \[\begin{split} I&=\frac{h}{3}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]\\ I_1&=\frac{h}{6}\left[f(a)+4f\left(a+\frac{h}{2}\right)+f(a+h)\right]\\ I_2&=\frac{h}{6}\left[f(a+h)+4f\left(a+\frac{3h}{2}\right)+f(b)\right] \end{split}\]

2. Se verifică dacă eroarea \(\varepsilon=\dfrac{1}{15}(I-I_1-I_2)\) este suficient de mică (de exemplu, mai mică decât un \(\varepsilon_0\) prestabilit, precum \(\varepsilon_0=10^{-6}\))

   În practică, se utilizează adesea numai jumătate din valoarea prestabilită, deci \(\dfrac{\varepsilon_0}{2}\) în loc de \(\varepsilon_0\), pentru a avea o marjă de eroare suficientă.

   Dacă nu am ajuns la un \(\varepsilon\) suficient de mic, se repetă algoritmul recursiv pentru fiecare subinterval din cele două.

Cod ilustrativ#

Vom considera deja cunoscută funcția simpson13_simplu ce calculează o integrală definită pe un domeniu dat folosind metoda Simpson 1/3 simplă (a se aminti codul aferent acesteia).

Pentru funcția cuad_adap ce ilustrează metoda cuadraturilor adaptive, vom utiliza:

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

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

• f - funcția ce trebuie integrată;

• eps - eroarea maximă acceptată, de exemplu \(10^{-6}\) - în practică, se va utiliza chiar jumătate din aceasta, așa cum am explicat în paragrafele anterioare.

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

function I = cuad_adap(a, b, f, eps)
    m = (a + b) / 2;
    I = simpson13_simplu(a, b, f);
    I1 = simpson13_simplu(a, m, f);
    I2 = simpson13_simplu(m, b, f);
    if (I - I1 - I2) / 15 > eps / 2
        I = cuad_adap(a, m, f) + cuad_adap(m, b, f);
    end
end

Licență#

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