Cuadraturi gaussiene#
Vom ruga cititorii să se familiarizeze mai întâi cu subcapitolul dedicat polinoamelor ortogonale, întrucât acestea se află la baza cuadraturilor gaussiene.
Totodată, vom face diverse trimiteri către interpolarea Lagrange, aceasta făcând parte din motorul utilizat de cuadraturile gaussiene.
Odată parcurse și înțelese, reveniți la acest subcapitol.
Explicația matematică#
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}^*\) și \(x_k\lt x_{k+1}\), \(\forall k\lt n\), o mulțime de \(n+1\) puncte distincte ordonate crescător, reprezentând evaluări ale funcției \(f\) ce vor fi alese optim ulterior.
Ne propunem să calculăm integrala:
\[ I_{\omega}(f)=\int_{-1}^1f(x)\omega(x)\,dx \]Pentru aceasta, să considerăm o aproximare de forma:
\[ \tilde{I}_{n,\omega}(f)=\sum_{k=0}^n\alpha_k f(x_k) \]unde \(\alpha_0,\dots\alpha_n\in\mathbb{R}\) sunt constante.
Folosind doar aceste două formule, devine evident că încercăm să minimizăm eroarea introdusă de aproximare, adică vrem să aducem următoarea funcție cât mai aproape de \(0\):
\[ \varepsilon_{n,\omega}(f)=I_{\omega}(f)-\tilde{I}_{n,\omega}(f) \]Dacă un polinom \(p\in\mathbb{P}_r\) (de grad \(r\in\mathbb{N}\)) respectă egalitatea \(\varepsilon_{n,\omega}(p)=0\), spunem despre aproximarea \(\tilde{I}_{n,\omega}\) că are gradul de valabilitate \(r\) în raport cu ponderea \(\omega\).
Utilizarea interpolării Lagrange#
Să considerăm interpolarea Lagrange a celor punctelor din \(M\), întrucât nu cunoaștem funcția \(f\) decât din punct de vedere numeric.
Notăm interpolarea cu \(L_f(x)=\sum_{k=0}^ny_k l_k(x)\) (păstrăm semnificațiile descrise în subcapitolul dedicat interpolării Lagrange); amintim că:
\[ l_k(x)=\prod_{\substack{0\leq j\leq n\\ i\neq j}}\left(\frac{x-x_j}{x_i-x_j}\right) \]Dacă vrem să calculăm integrala \(I_{\omega}(L_f)\) printr-o aproximare cu gradul de valabilitate (cel puțin) \(n\), putem găsi o formulă relativ ușor pornind de la integrala inițială:
\[\begin{split} I_{\omega}(L_f)&=\int_{-1}^1L_f(x)\omega(x)\,dx\\ &= \int_{-1}^1\left(\sum_{k=0}^ny_k l_k(x)\right)\omega(x)\,dx \end{split}\]Separând apoi integralele și însumându-le pe rând, obținem:
\[\begin{split} I_{\omega}(L_f)&=\sum_{k=0}^n\left(\int_{-1}^1y_k l_k(x)\omega(x)\,dx\right)\\ &=\sum_{k=0}^ny_k\left(\int_{-1}^1l_k(x)\omega(x)\,dx\right) \end{split}\]Aproximarea la care încercăm să ajungem are forma descrisă anterior:
\[ \tilde{I}_{n,\omega}(L_f)=\sum_{k=0}^n\alpha_k L_f(x_k)=\alpha_ky_k \]Așadar, este evident (prin identificare) că ar trebui să alegem \(\alpha_k\) astfel încât:
\[ \boxed{\alpha_k=\int_{-1}^1l_k(x)\omega(x)\,dx}\,,\,\,\forall k\in\overline{0,n} \]Setarea gradului de valabilitate#
Se poate demonstra matematic că, pentru o valoare prestabilită \(m\in\mathbb{N}^*\), aproximarea gaussiană a unei funcții \(M=\Big\{(x_k,y_k)\in\mathbb{R}^2\,|\,y_k=f(x_k),\,k\in\overline{0,n}\Big\}\), \(n\in\mathbb{N}^*\), are gradul de valabilitate egal cu \(n+m\) dacă și numai dacă se utilizează o tehnică de interpolare și se respectă egalitatea:
\[ \int_{-1}^1\left(\prod_{k=0}^n(x-x_k)\right)p(x)\omega(x)\,dx=0,\,\forall p\in\mathbb{P}_{m-1} \]Dacă notăm \(\Psi_{n+1}(x)=\prod_{k=0}^n(x-x_k)\), ecuația anterioară se poate scrie sub forma:
\[ \boxed{\int_{-1}^1\Psi_{n+1}(x)p(x)\omega(x)\,dx=0}\,,\,\,\forall p\in\mathbb{P}_{m-1} \]Cu alte cuvinte, polinoamele \(\Psi_{n+1}(x)\) și \(p(x)\) trebuie să fie ortogonale în raport cu o funcție pondere \(\omega\), indiferent de alegerea polinomului \(p\in\mathbb{P}_{m-1}\) de grad \(m-1\).
Voi schița doar demonstrația în sens invers, adică doar condiția suficientă.
Fie \(f\in\mathbb{P}_{n+m}\) un polinom (de grad \(n+m\)) oarecare. În acest context, prin teorema împărțirii cu rest a polinoamelor, funcția \(f\) poate fi scrisă ca o sumă compusă din produsul între polinomul \(\Psi_{n+1}=\prod_{k=0}^n(x-x_k)\) și un cât \(\pi_{m-1}\in\mathbb{P}_{m-1}\), respectiv un polinom rest, \(q_n\in\mathbb{P}_{n}\). Matematic:
\[ f(x)=\Psi_{n+1}(x)\pi_{m-1}(x) + q_n(x) \]Deoarece am demonstrat deja în paragrafele introductive faptul că interpolarea Lagrange are drept efect un grad de valabilitate al aproximării integralei egal cu cel puțin \(n\) (pentru \(n+1\) puncte), putem calcula integrala lui \(q_n\) astfel:
\[ \sum_{k=0}^n\alpha_kq_n(x_k)=\int_{-1}^1q_n(x)\omega(x)\,dx \]Știind că \(f(x)=\Psi_{n+1}(x)\pi_{m-1}(x)+q_n(x)\), obținem \(q_n(x)=f(x)-\Psi_{n+1}(x)\pi_{m-1}(x)\). Înlocuim acum în formula anterioară:
\[ \sum_{k=0}^n\alpha_kq_n(x_k)=\int_{-1}^1f(x)\omega(x)\,dx-\int_{-1}^1\Psi_{n+1}(x)\pi_{m-1}(x)\omega(x)\,dx \]Cu toate acestea, am început presupunând că \(\int_{-1}^1\Psi_{n+1}(x)p(x)\omega(x)\,dx=0\) (nu uitați că demonstrăm condiția suficientă!), deci integrala din partea dreaptă este nulă:
\[ \sum_{k=0}^n\alpha_kq_n(x_k)=\int_{-1}^1f(x)\omega(x)\,dx \]Totodată, integrala din partea dreaptă poate primi același tratament și poate fi aproximată folosind:
\[ \sum_{k=0}^{n+m}\beta_kf(x_k)=\int_{-1}^1f(x)\omega(x)\,dx \]unde \(\beta_0,\dots,\beta_{n+m}\in\mathbb{R}\) sunt constante. Cu toate acestea, egalasem inițial integrala cu un polinom de grad \(n\), așadar și acesta trebuie să fie un polinom de grad \(n\), deci \(\beta_{n+1}=\dots=\beta_{n+m}=0\), așadar:
\[ \sum_{k=0}^n\alpha_kq_n(x_k)=\sum_{k=0}^{n}\beta_kf(x_k) \]Cum funcția \(f\) este aleasă la întâmplare, iar constantele pot fi egalate, putem concluziona că eroarea \(\varepsilon_{n,\omega}\) este nulă, adică \(\varepsilon_{n,\omega}=0\), indiferent de polinomul \(f\in\mathbb{P}_{n+m}\).
Gradul maxim de valabilitate#
Drept corolar al teoremei anterioare, se garantează faptul că gradul maxim de valabilitate al unei aproximări de \(n+1\) puncte (a funcției \(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 \(2n+1\).
Vom presupune că se poate alege un \(m\) mai mare decât \(n+1\), adică \(m\geq n+2\), și vom demonstra că acest lucru conduce către o contradicție.
Setăm \(m=n+2\) și alegem polinomul \(p\in\mathbb{P}_{m-1}\Leftrightarrow p\in\mathbb{P}_{n+1}\) egal chiar cu \(\Psi_{n+1}\) - amintiți-vă că teorema trebuia să se respecte indiferent de alegerea lui \(p\), așadar, dacă se găsește un polinom care strică relația, este evident că am întâmpinat o contradicție.
Înlocuind acum în formula integralei ce trebuie egalată cu \(0\), se obține:
\[ \int_{-1}^1\Psi_{n+1}(x)p(x)\omega(x)\,dx=0\Leftrightarrow \int_{-1}^1\Psi^2_{n+1}(x)\omega(x)\,dx=0 \]Această ecuație implică însă că \(\Psi^2_{n+1}(x)\omega(x)=0\), deci că \(\Psi_{n+1}(x)=0\), \(\forall x\in[x_0,x_n]\). Această egalitate este însă absurdă - este suficient să alegem \(x\in(x_0,x_1)\) și să vedem că \(\Psi_{n+1}(x)=(x-x_0)(x-x_1)\dots\neq0\). Am ajuns la contradicția căutată.
Stabilirea punctelor de evaluare#
Pentru a obține gradul maxim de valabilitate al unei aproximări de integrală pentru o funcție \(M=\Big\{(x_k,y_k)\in\mathbb{R}^2\,|\,y_k=f(x_k),\,k\in\overline{0,n}\Big\}\) (puncte alese optim), unde \(n\in\mathbb{N}^*\), am observat că trebuie să forțăm:
\[ \int_{-1}^1\Psi_{n+1}(x)p(x)\omega(x)\,dx=0,\,\forall p\in\mathbb{P}_{n} \]Cu alte cuvinte, polinomul \(\Psi_{n+1}(x)\) de grad \(n+1\) trebuie să fie ortogonal față de orice alt polinom de grad mai mic decât acesta (cel mult \(n\)).
Astfel, concluzionăm că polinomul monic \(\Psi_{n+1}(x)\) poate fi scris drept un multiplu al unui polinom ortogonal de grad \(n+1\), să zicem \(P_{n+1}\). Cunoaștem însă rădăcinile lui \(\Psi_{n+1}(x)=\prod_{k=0}^n(x-x_k)\), așadar le cunoaștem și pe ale lui \(P_{n+1}\):
\[ \boxed{P_{n+1}(x_k)=0}\,,\,\,\forall k\in\overline{0,n} \]Astfel, folosind această ultimă formulă, se pot calcula acele valori \(x_0,\dots,x_n\) pentru care integrarea folosind cuadraturi gaussiene produce rezultate optime.
Conversia domeniului de integrare#
Toate explicațiile anterioare presupun calcularea unei integrale pe domeniul \([-1,1]\). În realitate însă, acest domeniu ar putea să difere, motiv pentru care se poate realiza foarte ușor o schimbare de variabilă.
De exemplu, să zicem că vrem să calculăm integrala:
\[ I=\int_a^bf(x)\,dx \]unde \(f:[a,b]\subset\mathbb{R}\rightarrow\mathbb{R}\) este o funcție integrabilă și continuă. În acest context, se poate realiza schimbarea de variabilă:
\[ u=\frac{2x-a-b}{b-a}\Leftrightarrow\boxed{x=\frac{(b-a)u+a+b}{2}\Rightarrow dx=\frac{b-a}{2}du} \]Astfel, întorcându-ne la integrala noastră, obținem:
\[ \boxed{I=\int_a^bf(x)\,dx=\frac{b-a}{2}\int_{-1}^1f\left(\frac{(b-a)u+a+b}{2}\right)\,du} \]Sinteza algoritmului#
Sintetizând toate informațiile anterioare, alegem o familie de polinoame ortogonale \(P_0,P_1,\dots\) în raport cu o funcție pondere \(\omega\), iar apoi definim următorul algoritm generic de rezolvare a integralelor definite:
-
Conversia domeniului. În primul rând, se transformă capetele de integrare astfel încât să se scrie integrala în funcție de domeniul \([-1,1]\). Cu alte cuvinte, se face trecerea:
\[ I=\int_a^bf(x)\,dx=\frac{b-a}{2}\int_{-1}^1f\left(\frac{(b-a)x+a+b}{2}\right)\,dx \] -
Rădăcinile polinomului ortogonal. Cunoscând polinoamele ortogonale cu care lucrăm, se calculează rădăcinile polinomului \(P_{n+1}\) - acestea vor reprezenta abscisele \(x_0,\dots,x_n\).
Din punct de vedere strict algoritmic, aceste valori sunt adesea precalculate și reținute ca atare (eventual, alternativ, memoizate).
-
Coeficienții \(\alpha_k\). Se calculează apoi coeficienții \(\alpha_0,\dots,\alpha_n\) utilizând formula:
\[ \alpha_k=\int_{-1}^1l_k(x)\omega(x)\,dx,\,\forall k\in\overline{0,n} \]Scriind definiția lui \(l_k\), se ajunge la formula:
\[ \alpha_k=\int_{-1}^1\Big[\prod_{\substack{0\leq j\leq n\\ i\neq j}}\frac{x-x_j}{x_i-x_j}\Big]\omega(x)\,dx,\,\forall k\in\overline{0,n} \]Aceste valori sunt de asemenea adesea precalculate și utilizate ca atare.
-
Calcularea integralei. Utilizând toate aceste valori, se poate calcula integrala:
\[ I=\frac{b-a}{2}\int_{-1}^1f\left(\frac{(b-a)x+a+b}{2}\right)\,dx\approx \frac{b-a}{2}\sum_{k=0}^n\alpha_kf\left(\frac{(b-a)x_k+a+b}{2}\right) \]
Cuadratura Gauss-Legendre#
Până la acest moment, am preferat o explicație matematică cât mai generală posibil. Propunem însă utilizarea polinoamelor Legendre normalizate (am descris pe scurt formarea acestor polinoame în subcapitolul dedicat polinoamelor ortogonale):
| \(n\) | Al \((n+1)\)-lea polinom Legendre | Polinomul normalizat | |
| 0 | \(1\) | \(1\) | |
| 1 | \(x\) | \(x\) | |
| 2 | \(\dfrac{1}{2}(3x^2-1)\) | \(x^2-\dfrac{2}{3}\) | |
| 3 | \(\dfrac{1}{2}(5x^3-3x)\) | \(x^3-\dfrac{6}{5}x\) | |
| 4 | \(\dfrac{1}{8}(35x^4-30x^2+3)\) | \(x^4-\dfrac{48}{7}x^2+\dfrac{24}{35}\) |
Acestea prezintă o formulă simplificată pentru coeficienții \(\alpha_0,\dots,\alpha_n\) pe care nu o vom demonstra aici[4]:
\[ \alpha_k=\frac{2}{(1-x_k^2)\left[P_n'(x_k)\right]^2},\,\forall k\in\overline{0,n} \]Vom ilustra în tabelul de mai jos primele câteva valori ale absciselor \(x_0,\dots,x_n\), respectiv ale ponderilor individuale \(\alpha_0,\dots,\alpha_n\) (toate valorile se pot obține prin calcul brut):
| \(n\) | Valorile absciselor (\(x_k\)) | Valorile ponderilor individuale (\(\alpha_k\)) | |
| 0 | \(0\) | \(2\) | |
| 1 | \(-\dfrac{1}{\sqrt{3}}\) | \(1\) | |
| \(-\dfrac{1}{\sqrt{3}}\) | \(1\) | ||
| 2 | \(-\sqrt{\dfrac{3}{5}}\) | \(\dfrac{5}{9}\) | |
| \(0\) | \(\dfrac{8}{9}\) | ||
| \(-\sqrt{\dfrac{3}{5}}\) | \(\dfrac{5}{9}\) |
Aceste valori se continuă la infinit, însă nu suntem interesați la acest moment să dezvoltăm mai mult. Detalii suplimentare se găsesc oricum în subcapitolul dedicat, după cum v-am obișnuit.
Licență#
The book "Metode Numerice", written by Valentin-Ioan Vintilă, is licensed under CC BY-NC-SA 4.0 
