Metoda Simpson#

Elementele geometrice utilizate până la acest moment pentru a aproxima arii au fost cele mai simple cu putință - există însă metode mai bune de a găsi rezultatul căutat. Spre exemplu, metoda Simpson propune utilizarea ariilor aflate sub parabolele unor ecuații de grad 2.

Vom detalia aceasă descoperire în paragrafele ce urmează.

Metoda Simpson 1/3#

Există mai multe metode Simpson - începem cu cea mai simplă.

Fie \(f:[a,b]\subset\mathbb{R}\rightarrow\mathbb{R}\) o funcție integrabilă și continuă pe întreg intervalul de definiție. Totodată, fie \(m\in\mathbb{R}\) ales astfel încât \(m=\frac{a+b}{2}\).

Metoda Simpson utilizează aria aflată sub graficul unei ecuații de grad 2 pentru a aproxima integrala căutată. Așadar, pentru a aproxima \(\int_a^bf(x)\,dx\), se începe prin descoperirea unei parabole care trece exact prin punctele \(\big(a,f(a)\big)\), \(\big(m, f(m)\big)\) și \(\big(b,f(b)\big)\).

Pentru aceasta, se poate folosi, de exemplu, interpolarea polinomială Lagrange - se obține polinomul \(L\in\mathbb{R}[x]\):

\[ L(x)=\frac{(x-m)(x-b)}{(a-m)(a-b)}f(a)+\frac{(x-a)(x-b)}{(m-a)(m-b)}f(m)+\frac{(x-a)(x-m)}{(b-a)(b-m)}f(b) \]

Așadar, \(f(x)\approx L(x)\), \(\forall x\in\mathbb[a,b]\), deci \(\int_a^bf(x)\,dx\approx \int_a^bL(x)\,dx\), dar, spre deosebire de \(\int_a^bf(x)\,dx\), \(\int_a^bL(x)\,dx\) poate fi calculat ușor, fiind integrarea unui polinom.

Într-un final, se obține această aproximare pentru integrala inițială:

\[ \boxed{\int_a^bf(x)\,dx\approx\frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]} \]
Demonstrație

Să încercăm să calculăm \(\int_a^bf(x)\,dx\):

\[\begin{split} \int_a^bf(x)\,dx&=\int_a^b\frac{(x-m)(x-b)}{(a-m)(a-b)}f(a)+\frac{(x-a)(x-b)}{(m-a)(m-b)}f(m)+\frac{(x-a)(x-m)}{(b-a)(b-m)}f(b)\,dx \end{split}\]

Luăm termenii pe rând. Pentru început:

\[\begin{split} \int_a^b\frac{(x-m)(x-b)}{(a-m)(a-b)}f(a)\,dx&=f(a)\int_a^b\frac{(x-m)(x-b)}{(a-m)(a-b)}\,dx\\ &=f(a)\int_a^b\frac{\left(x-\frac{a+b}{2}\right)(x-b)}{\left(a-\frac{a+b}{2}\right)(a-b)}\,dx \end{split}\]

Vom face schimbarea de variabilă \(x=\frac{b-a}{2}u+\frac{a+b}{2}\), adică \(u=\frac{2x-a-b}{b-a}\), ceea ce înseamnă că \(dx=\frac{b-a}{2}du\). Întorcându-ne la integrala noastră:

\[\begin{split} f(a)\int_a^b\frac{\left(x-\frac{a+b}{2}\right)(x-b)}{\left(a-\frac{a+b}{2}\right)(a-b)}\,dx&=f(a)\int_{-1}^1\frac{\left(\frac{b-a}{2}u+\frac{a+b}{2}-\frac{a+b}{2}\right)(\frac{b-a}{2}u+\frac{a+b}{2}-b)}{\left(a-\frac{a+b}{2}\right)(a-b)}\cdot\frac{b-a}{2}\,du\\ &=f(a)\int_{-1}^1\frac{\left(\frac{a-b}{2}\right)^2u\left(1-u\right)}{\frac{a-b}{2}(a-b)}\cdot\frac{a-b}{2}\,du\\ &=\frac{a-b}{4}f(a)\int_{-1}^1u\left(1-u\right)\,du\\ &=\frac{a-b}{4}f(a)\left[\frac{x^2}{2}-\frac{x^3}{3}\right]_{-1}^1\\ &=\frac{a-b}{4}f(a)\left(\frac{1^2}{2}-\frac{1^3}{3}-\frac{(-1)^2}{2}+\frac{(-1)^3}{3}\right)\\ &=\frac{a-b}{4}\cdot\frac{-2}{3}f(a)\\ &=\frac{b-a}{6}f(a) \end{split}\]

Folosind același raționament și pentru celelalte integrale, am ajunge la formula enunțată anterior:

\[ \boxed{\int_a^bf(x)\,dx\approx\frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]} \]

Această formulă poate fi îmbunătățită prin compunere, asemănător cu ceea ce am discutat în subcapitolele anterioare. Astfel, pentru mai multe intervale consecutive, se însumează toate aceste aproximări.

Concret, fie \(f:[x_0,x_n]\subset\mathbb{R}\rightarrow\mathbb{R}\) o funcție integrabilă și continuă pe întreg intervalul de definiție. Totodată, fie mulțimea \(M=\Big\{(x_k,y_k)\in\mathbb{R}^2\,|\,y_k=f(x_k),\,k\in\overline{0,n}\Big\}\), unde \(n\in\mathbb{N}^*\) este par și \(x_k\lt x_{k+1}\), \(\forall k\lt n\), o mulțime de \(n+1\) puncte distincte și echidistante ordonate crescător, reprezentând evaluări ale funcției \(f\) (concluzionăm deci că metoda Simpson 1/3 este o metodă de tip Newton-Cotes).

Dacă notăm cu \(h=\frac{b-a}{n}\) distanța între oricare două puncte \(x_k\) și \(x_{k+1}\) consecutive, \(\forall k\in\overline{0,n-1}\), obținem următoarele două formule:

\[ \boxed{\begin{split} \int_{x_0}^{x_n} f(x) \, dx &\approx \frac{h}{3} \sum_{k=1}^{n/2}\big[f(x_{2k-2}) + 4f(x_{2k-1}) + f(x_{2k})\big] \\ &= \frac{h}{3} \bigg[f(x_0) + 4\sum_{k=1}^{n/2} f(x_{2k-1}) + 2\sum_{k=1}^{n/2-1} f(x_{2k}) + f(x_n)\bigg] \end{split}} \]

Vizualizare#

În spiritul competiției, să considerăm tot funcția \(f:[0,6]\rightarrow\mathbb{R}\), \(f(x)=\sin(x+2)+\sqrt{x}\). Vom ilustra metoda Simpson 1/3 compusă folosind \(3\), \(5\), respectiv \(7\) puncte echidistante pe întreg intervalul \([0,6]\):

O parabolă (trei puncte)

Două parabole (cinci puncte)

Trei parabole (șapte puncte)

Observăm că rezultatele obținute prin această aproximare iau următoarea formă:

\[ \int_0^6\sin(x+2)+\sqrt{x}\,dx\approx\begin{cases} \text{o parabolă:}&7.441\\ \text{două parabole:}&9.368\\ \text{trei parabole:}&9.444\\ \end{cases} \]

Amintim rezultatul corect, \(9.527\), respectiv aproximările obținute prin metoda trapezelor compuse:

\[ \int_0^6\sin(x+2)+\sqrt{x}\,dx\approx\begin{cases} \text{un trapez:}&13.044\\ \text{trei trapeze:}&9.104\\ \text{șase trapeze:}&9.359\\ \end{cases} \]

Se observă deci capacitatea acestui algoritm de a calcula rezultatul mai rapid, în ciuda utilizării unui număr similar de puncte în aproximarea sa.

Estimarea erorii parțiale#

Fie \([a,b]\) un interval luat în considerare pentru a se aplica, pe o funcție \(f\), metoda Simpson 1/3. Atunci, eroarea metodei este dată de formula lui Stirling:

\[ \boxed{\varepsilon_i(a,b)\approx-\frac{f^{(4)}(\xi)}{180}(b-a)h^4=-\frac{f^{(4)}(\xi)}{90}h^5}\,,\,\,\xi\in(a,b) \]

unde \(h\in\mathbb{R}\) reprezintă media numerelor \(a\) și \(b\). Folosind teorema lui Lagrange, există cel puțin un punct \(\xi\in(a,b)\) astfel încât:

\[ f^{(4)}(\xi)=\left(f^{(3)}(\xi)\right)'=\frac{f^{(3)}(b)-f^{(3)}(a)}{b-a} \]

deci ecuația se mai poate scrie și sub forma:

\[ \boxed{\varepsilon_i(a,b)\approx\frac{f^{(3)}(a)-f^{(3)}(b)}{180}h^4}\,,\,\forall x_0\in(a,b) \]
Demonstrație

Dovada este inspirată de munca lui Scarborough[51].

Fie \(F(x)\) definită astfel încât:

\[ \int_a^xf(x)\,dx=F(x)+\mathcal{C} \]

Atunci, pentru a calcula \(\int_{k-h}^{k+h}f(x)\,dx\), unde \(k-h\leq x\leq k+h\) și \(a+h\leq k\leq b-h\), se poate scrie formula următoare:

\[ \int_{k-h}^{k+h}f(x)\,dx=F(k+h)-F(k-h) \]

Vom transforma în serie Taylor termenii \(F(k+h)\) și \(F(k-h)\), amintind că \(F'(x)=f(x)\):

\[\begin{split} F(k+h)&=F(k)+f(k)h+\frac{f'(k)}{2!}h^2+\frac{f''(k)}{3!}h^3+\dots\\ F(k-h)&=F(k)-f(k)h+\frac{f'(k)}{2!}h^2-\frac{f''(k)}{3!}h^3+\dots \end{split}\]

Din aceste două ecuații, înlocuind în cea de mai sus, se obține:

\[ \int_{k-h}^{k+h}f(x)\,dx=2\left(f(k)h+\frac{f''(k)}{3!}h^3+\frac{f^{(4)}(k)}{5!}h^5+\dots\right) \]

Reținem această ecuație. Folosind metoda lui Simpson, valoarea integralei este dată de:

\[ \int_{k-h}^{k+h}f(x)\,dx\approx\frac{h}{3}\left[f(k-h)+4f(k)+f(k+h)\right] \]

Prin același artificiu de a scrie în serie Taylor funcțiile \(f(k+h)\) și \(f(k-h)\), se obține:

\[ \int_{k-h}^{k+h}f(x)\,dx\approx\frac{h}{3}\left[6f(k)+f''(k)h^2+\frac{2f^{(4)}(k)}{4!}h^4+\dots\right] \]

Prin scăderea acestei relații din cea reținută și considerând \(a=k+h\) și \(b=k-h\), se obține:

\[ \varepsilon_i(a,b)\approx-\frac{f^{(4)}(k)}{90}h^5 \]

Din moment ce \(h=\dfrac{a+b}{2}\), aici se încheie demonstrația.

Estimarea erorii totale#

Folosind ecuația erorii parțiale pentru întreg intervalul \([a,b]=[x_0,x_n]\), se poate calcula o formulă pentru obținerea unei erori absolute maxime a metodei Simpson 1/3 compusă:

\[ \varepsilon(f,h)=-\frac{x_n-x_0}{180}h^4\cdot f^{(4)}(\xi)\Rightarrow\boxed{\varepsilon(f,h)\leq\frac{x_n-x_0}{180}h^4\cdot\max_{\xi\in[x_0,x_n]}\left\{|f^{(4)}(\xi)|\right\}} \]

Aceasta nu va fi demonstrată, ci se va accepta drept adevărată ca atare.

Coduri ilustrative#

Vom separa cazul implementării algoritmului Simpson simplu de cel compus.

Regula lui Simpson simplă#

Funcția simpson13_simplu prezintă, în Matlab, o implementare posibilă a metodei Simpson 1/3 simple, așa cum a fost explicată mai sus. Vom considera:

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

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

  • 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 = simpson13_simplu(a, b, f)
    I = (b - a) * (f(a) + 4 * f((a+b)/2) + f(b)) / 6;
end

Metoda lui Simpson compusă#

Funcția simpson13_compus va conține o implementare în Matlab pentru metoda Simpson 1/3 compusă, în urma explicațiilor de mai sus. În acest context, folosim parametrii:

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

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

  • n - numărul de parabole utilizat;

  • 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 = simpson13_compus(a, b, n, f)	
    h = (b - a) / (2 * n);

    s1 = 0;
    for i = [1 : n]
        s1 += f(a + (2*i-1) * h);
    end

    s2 = 0;
    for i = [1 : n-1]
        s2 += f(a + 2*i * h);
    end

    I = h * (f(a) + f(b) + 4*s1 + 2*s2) / 3;
end

(Bonus) Metoda Simpson 3/8#

Ideea din spatele metodei Simpson 1/3 a fost revoluționară pentru că a utilizat interpolările unor polinoame de gradul \(2\); această noțiune se poate extinde pe polinoame de grad 3 și este cunoscută drept metoda Simpson 3/8.

Pe scurt, în cazul metodei compuse: fie \(f:[x_0,x_n]\subset\mathbb{R}\rightarrow\mathbb{R}\) o funcție integrabilă și continuă pe întreg intervalul de definiție.

Totodată, fie mulțimea \(M=\Big\{(x_k,y_k)\in\mathbb{R}^2\,|\,y_k=f(x_k),\,k\in\overline{0,n}\Big\}\), unde \(n\in\mathbb{N}^*\) este divizibil cu \(3\) și \(x_k\lt x_{k+1}\), \(\forall k\lt n\), o mulțime de \(n+1\) puncte distincte și echidistante ordonate crescător, reprezentând evaluări ale funcției \(f\) - deci o altă metodă tip Newton-Cotes.

Dacă notăm cu \(h=\frac{b-a}{n}\) distanța între oricare două abscise \(x_k\) și \(x_{k+1}\) consecutive, \(\forall k\in\overline{0,n-1}\), obținem următoarea formulă pentru metoda compusă:

\[ \boxed{\begin{split} \int_{x_0}^{x_n} f(x) \, dx &\approx \frac{3h}{8}\bigg[f(x_0)+3f(x_1)+3f(x_2)+2f(x_3)+3f(x_4)+3f(x_5)+2f(x_6)+ \\ &\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\dots+3f(x_{n-2}) + 3f(x_{n-1}) + f(x_n)\bigg]\\ &= \frac{3h}{8} \bigg[f(x_0) + 3 \sum_{i \ne 3k}^{n-1} f(x_i) + 2 \sum_{j=1}^{n/3 - 1} f(x_{3j}) + f(x_n) \bigg] \end{split}} \]

Estimarea erorii#

Nu vom intra în detalii, însă eroarea provocată de această metodă la nivel local (pe \([a,b]\subset\mathbb{R}\)) este dată de formula:

\[ \varepsilon(x)=-\frac{h^4}{405}(b-a)\cdot f^{(4)}(\xi) \]

unde \(\xi\in[a,b]\) este necunoscut și \(h=\frac{b-a}{2}\). Evident, o limită superioară a erorii absolute poate fi determinată ușor:

\[ \boxed{|\varepsilon(x)|\leq\frac{h^4}{405}(b-a)\cdot\max_{\xi\in[a,b]}\left|f^{(4)}(\xi)\right|} \]

Concluzii#

Metoda Simpson 1/3 compusă este, atât în teorie, cât și în practică, mult mai rapidă, convergând de câteva ori bune mai ușor decât metoda trapezelor compuse. Singurul dezavantaj al acestei metode este creșterea numărului de puncte ce sunt necesare pentru a calcula estimarea.

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

Results will appear as you type