7. Interpolarea Lagrange#

O altă modalitate de a interpola o funcție construind un polinom este interpolarea Lagrange. Vom încerca să o explicăm, alături de proveniența sa.

Vom propune mai întâi conceptul de multiplicator Lagrange, vom trece apoi prin interpolarea propriu-zisă, o vom vizualiza și îi vom determina eroarea numerică. Într-un final, vom exemplifica un cod Matlab echivalent acestei metode.

Multiplicatorul Lagrange#

Pornim tot de la o mulțime \(M=\Big\{(x_k,y_k)\in\mathbb{R}^2\,|\,k\in\overline{0,n}\Big\}\), unde \(n\in\mathbb{N}^*\) și \(x_k\lt x_{k+1}\), \(\forall k\lt n\), mulțime ce conține \(n+1\) puncte distincte ordonate crescător față de \(x_k\).

Prin definiție, multiplicatorul Lagrange \(L_k:[x_0,x_n]\rightarrow\mathbb{R}\), este acea funcție ce respectă formula:

\[\begin{split} L_k(x)&=\left(\frac{x-x_0}{x_k-x_0}\right)\dots\left(\frac{x-x_{k-1}}{x_k-x_{k-1}}\right)\left(\frac{x-x_{k+1}}{x_k-x_{k+1}}\right)\dots\left(\frac{x-x_n}{x_k-x_n}\right)\\ &=\prod_{\substack{0\leq j\leq n\\ i\neq j}}\left(\frac{x-x_j}{x_i-x_j}\right) \end{split}\]

Motivul pentru care această structură prezintă interes este următorul rezultat:

\[ \boxed{ L_k(x)=\begin{cases} 1,&\text{dacă } x=x_k\\ 0,&\text{dacă } x\neq x_k \end{cases} }\,,\,\,\forall x\in\{x_0,\dots,x_n\} \]

Cu alte cuvinte, un multiplicator Lagrange nu poate lua decât două valori, \(0\) și \(1\), bine determinate.

Interpolarea Lagrange#

Fie \(M=\Big\{(x_k,y_k)\in\mathbb{R}^2\,|\,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.

Datorită proprietăților multiplicatorului Lagrange discutate anterior, putem determina următoarea egalitate cu ușurință:

\[ y_k L_k(x)=\begin{cases} y_k,&\text{dacă } x=x_k\\ 0,&\text{dacă } x\neq x_k \end{cases} \,,\,\forall x\in\{x_0,\dots,x_n\} \]

Din aceasta, răsare polinomul \(L:[x_0,x_n]\rightarrow\mathbb{R}\), ca o însumare de \(n+1\) astfel de relații:

\[ \boxed{L(x)=\sum_{k=0}^ny_k L_k(x)} \]

Conform ecuației de mai sus, polinomul \(L\) va trece prin toate punctele mulțimii \(M\), deci va reprezenta o interpolare a mulțimii \(M\). Această formulă este o îmbunătățire numerică a interpolării Vandermonde.

Vizualizare#

Vom demonstra că această interpolare funcționează și vom compara rezultatul obținut cu interpolarea Vandermonde pentru funcția:

\[ f(x)=\frac{\cot(x)}{1+64\left(x-1\right)^{2}} \]

Aceasta este aceeași funcție descrisă în exemplul precedent. Asemănător, ne va interesa exclusiv intervalul \([0.1,1.6]\) și vom interpola cu \(4\), \(7\) și \(9\) puncte echidistante funcția:

Interpolările cu \(4\), \(7\) și \(9\) puncte

Nu este o greșeală! Rezultatul poate fi comparat cu interpolarea Vandermonde și se va observa că arată absolut identic. Oare așa trebuia să se întâmple sau am greșit ceva pe parcurs?

Teoremă - polinom unic

Așadar, interpolările polinomiale produc, în teorie, rezultate identice - diferența este dată exclusiv de eroarea numerică. În cazul interpolării Lagrange, aceasta este mai mică față de interpolarea Vandermonde.

Eroarea interpolării#

Să o calculăm riguros.

Fie \(f:[a,b]\subset\mathbb{R}\rightarrow\mathbb{R}\) o funcție de clasă \(C^{n+1}\), unde \(n\in\mathbb{N}^*\). 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\}\), \(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\).

Atunci, dacă \(L(x)\) reprezintă polinomul Lagrange ce interpolează toate punctele mulțimii \(M\), unde \(L:[x_0,x_n]\to\mathbb{R}\), are loc relația:

\[ \boxed{f(x)=L(x)+\frac{f^{(n+1)}\big(\xi(x)\big)}{(n+1)!}\prod_{k=0}^n(x-x_k)} \]

Se poate astfel deduce că eroarea interpolării Lagrange este dată de formula:

\[ \varepsilon(x)=f(x)-L(x)=\frac{f^{(n+1)}\big(\xi(x)\big)}{(n+1)!}\prod_{k=0}^n(x-x_k) \]

Menționăm că prin \(\xi(x)\) se înțelege o valoare (în general necunoscută) ce apare ca o consecință a teoremei lui Lagrange, însă se cunoaște faptul că \(\xi(x)\in[x_0,x_n]\).

Demonstrație

Vom demonstra relația dintre \(f\) și \(L\) începând cu valorile \(x=x_k\), \(\forall k\in\overline{0,n}\).

Se poate observa că, făcând înlocuirea direct în formulă, se obține \(f(x)=L(x)\), propoziție obligatoriu adevărată, căci polinomul \(L\) reprezintă interpolarea mulțimii \(M\), adică \(f(x_k)=L(x_k)\).

Dacă \(x\neq x_k\), \(\forall k\in\overline{0,n}\), definim funcția \(\gamma:[x_0,x_n]\rightarrow\mathbb{R}\), \(\gamma\in C^{n+1}\), astfel:

\[ \gamma(t)=\Big[f(t)-P(t)\Big]-\Big[f(x)-P(x)\Big]\cdot\frac{\prod_{k=0}^n(t-x_k)}{\prod_{k=0}^n(x-x_k)} \]

Observăm că \(\gamma(x_k)=0\), \(\forall k\in\overline{0,n}\). Totodată, observăm că \(\gamma(x)=0\), \(\forall x\in[x_0,x_n]\). Așadar, funcția \(\gamma\) se anulează în toate punctele mulțimii \(X=\{x,x_0,\dots,x_n\}\).

Folosind teorema generalizată a lui Rolle, se demonstrează că există o valoare \(\xi(x)\in(x_0,x_n)\) pentru care \(\gamma^{(n+1)}(\xi)=0\). Continuăm pe această idee:

\[\begin{split} \gamma^{(n+1)}(\xi)&=\Big[f^{(n+1)}(\xi)-P^{(n+1)}(\xi)\Big]-\Big[f(x)-P(x)\Big]\cdot\frac{d^{n+1}}{dt^{n+1}}\left[\frac{\prod_{k=0}^n(t-x_k)}{\prod_{k=0}^n(x-x_k)}\right]_{t=\xi}\\ 0&=f^{(n+1)}(\xi)-\Big[f(x)-P(x)\Big]\cdot\frac{(n+1)!}{\prod_{k=0}^n(x-x_k)} \end{split}\]

Din această ecuație, se poate scoate \(f(x)\) și se va obține relația ce trebuia dovedită.

Lagrange vs. Vandermonde#

Cele două interpolări produc rezultate similare din punct de vedere numeric, atâta vreme cât se interpolează un număr mic de puncte.

Cu toate acestea, deoarece în cadrul interpolării Vandermonde se rezolvă un sistem de ecuații liniare, rezultatul va fi mai instabil numeric decât în situația interpolării Lagrange.

Așadar, deși pentru puține puncte nu există nicio diferență, pentru un număr mare de puncte, utilizarea interpolării Vandermonde devine periculoasă.

Cod ilustrativ#

Codul de mai jos reprezintă o implementare a evaluării interpolării Lagrange în mai multe puncte (funcția lag_multiple) utilizând o funcție de evaluare într-un singur punct (lag_single). Vom considera:

x - abscisele mulțimii \(M\);
y - ordonatele mulțimii \(M\);
xps - punctele (abscisele) în care se vrea evaluarea multiplă;
xp - punctul (abscisa) în care se vrea unica evaluare.

Funcțiile returnează yps, rezultatele evaluărilor multiple, respectiv yp, unica evaluare a polinomului Lagrange.

function yps = lag_multiple(x, y, xps)
    for i = [1 : length(xps)]
        yps(i) = lag_single(x, y, xps(i));
    end
end

function yp = lag_single(x, y, xp)
    n = length(x);
    yp = 0;
    for i = [1 : n]
        l = 1;
        for j = [1 : n]
            if i != j
                l *= (xp - x(j)) / (x(i) - x(j));
            end
        end
        yp += y(i) * l;
    end
end

Forma matriceală#

Datorită eficienței lucrului cu matrice în Matlab, vom dezvolta și o variantă matriceală a interpolării, echivalentă cu metoda prezentată anterior.

Să definim \(n+1\) funcții, \(v_k(x)\), unde \(v_k:\{x_0,\dots,x_n\}\rightarrow\mathbb{R}^{n+1}\), și \(n+1\) vectori, \(u_k\in\mathbb{R}^{n+1}\), iar \(k\in\overline{0,n}\):

\[ v_k(x)=\begin{bmatrix} x-x_0\\ x-x_1\\ \vdots\\ x-x_{k-1}\\ 1\\ x-x_{k+1}\\ \vdots\\ x-x_n \end{bmatrix}\text{ și } u_k=\begin{bmatrix} x_k-x_0\\ x_k-x_1\\ \vdots\\ x_k-x_{k-1}\\ 1\\ x_k-x_{k+1}\\ \vdots\\ x_k-x_n \end{bmatrix}\,,\,\forall k\in\overline{0,n} \]

Observăm că, înmulțind toate elementele din \(v_k(x)\) și împărțind rezultatul la produsul elementelor din \(u_k\), obținem chiar \(L_k(x)\).

Așadar, folosind aceste definiții, obținem:

\[ \boxed{L_k(x)=\frac{\prod_{i\in\overline{0,n}}\Big(v_k(x)\Big)_i}{\prod_{i\in\overline{0,n}}\Big(u_k\Big)_i}} \]

Să construim acum două matrice, \(V(x)\), \(U\in\mathbb{R}^{(n+1)\times(n+1)}\), folosind vectorii definiți anterior. Matricea \(V\) este de fapt o funcție ce îl are pe \(x\) drept parametru. Așadar, acestea arată astfel:

\[ \begin{cases} V(x)&=\begin{bmatrix} v_0(x) & \dots & v_n(x) \end{bmatrix}\\ U&=\begin{bmatrix} u_0 & \dots & u_n \end{bmatrix} \end{cases} \]

Alternativ, scrise desfășurat:

\[ \begin{cases} V(x)=\begin{bmatrix} 1 & x-x_0 & \dots & x-x_0\\ x-x_1 & 1 & \dots & x-x_1\\ \vdots & \vdots & \ddots & \vdots\\ x-x_n & x-x_n & \dots & 1 \end{bmatrix}\\\\ U=\begin{bmatrix} 1 & x_1-x_0 & \dots & x_n-x_0\\ x_0-x_1 & 1 & \dots & x_n-x_1\\ \vdots & \vdots & \ddots & \vdots\\ x_0-x_n & x_1-x_n & \dots & 1 \end{bmatrix} \end{cases} \]

Folosind matricele anterioare, se poate demonstra relația:

\[ L_k(x)=\frac{\prod_{i\in\overline{0,n}}\Big(v_k(x)\Big)_i}{\prod_{i\in\overline{0,n}}\Big(u_k\Big)_i}=\frac{\prod_{i\in\overline{0,n}}\Big(V(x)_k\Big)_i}{\prod_{i\in\overline{0,n}}\Big(U_k\Big)_i} \]

Vrem acum o modalitate de a scrie matricea:

\[ A(x)=\begin{bmatrix} L_0(x) & L_1(x) & \dots & L_n(x) \end{bmatrix} \]

Strict din punct de vedere al unui pseudocod Octave, se obține formula:

\[ \boxed{A(x)= \text{prod}(V) \text{ ./ } \text{prod}(U)} \]

În acest moment, evaluarea \(L(x)\) se face prin produsul următor:

\[ L(x)=A\cdot\begin{bmatrix}y_0&y_1&\dots&y_n\end{bmatrix}\Rightarrow\text{ prin notație: } L(x)=A\bm{y} \]

Cod matriceal#

Implementarea sub formă matriceală a interpolării Lagrange poate fi făcută prin funcția lag_mat, ce utilizează:

x - abscisele mulțimii \(M\);
y - ordonatele mulțimii \(M\);
xp - punctul (abscisa) în care se vrea unica evaluare.

Funcția returnează yp, evaluarea polinomului Lagrange în punctul de abscisă xp.

function yp = lag_mat(x, y, xp)	
    n = length(x);
    U = diag(x) * ones(n) - ones(n) * diag(x) + eye(n);
    V = diag(xp - x) * ~eye(n) + eye(n);
    yp = prod(V) ./ prod(U) * y;
end

Concluzii#

Această metodă este robustă în comparație cu soluția Vandermonde, însă rămâne instabilă numeric, întrucât oscilează în capete. Nu se va folosi în practică niciodată această metodă pentru mai mult de \(50\) puncte.

Licență#

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