Interpolarea Newton#

Există un algoritm asemănător algoritmului Neville, și anume interpolarea cu diferențe divizate, cunoscută și sub denumirea de interpolare Newton.

Această metodă este folosită atunci când se vrea întregului polinom unic ce interpolează o mulțime de puncte, nu doar o evaluare a sa într-un punct fix.

Pentru a o putea explica însă, va fi nevoie să explicăm ce sunt acelea diferențele divizate și de ce pot fi utile în calculele noastre.

Diferențe divizate#

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.

În acest context, se definesc recursiv diferențele divizate ale acestei funcții numerice (din punct de vedere numeric, așa cum am tot arătat pe parcursul acestui capitol, funcțiile sunt doar o mulțime de puncte) folosind formula:

\[ \boxed{[y_i,\dots,y_j]=\begin{cases} y_i,&\text{dacă } i=j\\ \dfrac{[y_{i+1},\dots,y_{j}]-[y_{i},\dots,y_{j-1}]}{x_{j}-x_i},&\text{dacă } i\lt j \end{cases}}\,,\,\,\forall i\in\overline{0,n},\,j\in\overline{i,n} \]

Dacă ar fi să scriem relația pentru o funcție \(f\) care interpolează toate aceste puncte, se mai folosește și notația \(f[x_i,\dots,x_j]\) pentru diferența divizată \([y_i,\dots,y_j]\).

Polinoamele Newton de bază#

Definim polinoamele Newton de bază prin formula:

\[ \boxed{\pi_k(x)=(x-x_0)(x-x_1)\dots(x-x_{k-1})=\prod_{k=0}^{j-1}(x-x_k)}\,\,,\,\forall k\in\overline{0,n} \]

Observăm cum există o legătură evidentă între acestea și multiplicatorii Lagrange.

Fără a intra în detalii suplimentare, avem tot ce ne trebuie pentru a putea descoperi interpolarea Newton.

Interpolarea Newton#

Încercăm acum să construim un polinom \(P_n\), interpolarea căutată, în baza unei combinații liniare de polinoame Newton de bază:

\[ P_n(x)=\alpha_0+\alpha_1\pi_1(x)+\dots+\alpha_n\pi_{n}(x)=\sum_{k=0}^n\alpha_k\pi_k \]

unde \(\alpha_0\dots\alpha_n\in\mathbb{R}\) sunt constante. Vom construi însă acest polinom treptat, recursiv.

Construcția recursivă#

Începem cu \(P_0\), un polinom de grad \(0\) care interpolează doar punctul \((x_0,y_0)\). Evident:

\[ P_0(x)=y_0 \]

Trecem la \(P_1\), un polinom de grad \(1\) ce interpolează \((x_0,y_0)\) și \((x_1,y_1)\). Pentru a genera o recurență, vom încerca să îl scriem pe \(P_1\) în funcție de \(P_0\) și de un alt polinom, \(Q_1\):

\[ P_1(x)=P_0(x)+Q_1(x) \]

Pentru că \(P_1\) trebuie să treacă prin \((x_0,y_0)\), știm că \(P(x_0)=y_0\). Totodată, deoarece așa am definit polinomul \(P_0\), cunoaștem \(P_0(x_0)=y_0\). Urmând acest raționament, deducem că \(Q_1(x_0)=0\), ceea ce ne face să credem că poate lua forma:

\[ Q_1(x)=\alpha_1(x-x_0) \]

Ne întoarcem în \(P_1\) evaluat în \(x_1\) și extragem valoarea lui \(\alpha_1\):

\[ P_1(x_1)=P_0(x_1)+\alpha_1(x_1-x_0)\Rightarrow\boxed{\alpha_1=\frac{y_1-P_0(x_1)}{x_1-x_0}} \]

Ne îndreptăm acum atenția asupra lui \(P_2\), un polinom de grad \(2\) ce interpolează \((x_0,y_0)\), \((x_1,y_1)\) și \((x_2,y_2)\). În baza celor discutate anterior, polinomul poate fi scris în funcție de \(P_1\) și de un alt polinom, \(Q_2\):

\[ P_2(x)=P_1(x)+Q_2(x) \]

Ne vom folosi de punctele \((x_0,y_0)\) și \((x_1,y_1)\) pentru a afla detalii suplimentare despre \(Q_2\):

\[ \begin{cases} P_2(x_0)=P_1(x_0)+Q_2(x_0)\Rightarrow y_0=y_0+Q_2(x_0)\\ P_2(x_1)=P_1(x_1)+Q_2(x_1)\Rightarrow y_1=y_1+Q_2(x_1) \end{cases}\Rightarrow Q_2(x_0)=Q_2(x_1)=0 \]

Concluzionăm deci că \(Q_2\) ia forma:

\[ Q_2(x)=\alpha_2(x-x_0)(x-x_1) \]

Ne întoarcem în \(P_2\) evaluat în \(x_2\) și extragem valoarea lui \(\alpha_2\):

\[ P_2(x_2)=P_1(x_2)+\alpha_2(x_2-x_0)(x_2-x_1)\Rightarrow\boxed{\alpha_2=\frac{y_2-P_1(x_2)}{(x_2-x_0)(x_2-x_1)}} \]

Generalizând acest raționament pentru polinomul \(P_k\) (riguros, s-ar folosi inducția), observăm cum acesta trebuie să interpoleze punctele \((x_0,y_0),\dots,(x_k,y_k)\), așadar \(P_k(x_i)=P_{k-1}(x_i)+Q_k(x_i)\) se va reduce la \(Q_k(x_i)=0\), \(\forall i\in\overline{0,k-1}\). Din acest motiv, \(Q_k\) va lua forma:

\[ Q_k(x)=\alpha_k(x-x_0)\dots(x-x_{k-1}) \]

Când se evaluează \(P_k\) în \(x_k\), se poate extrage valoarea lui \(\alpha_k\):

\[ P_k(x_k)=P_{k-1}(x_k)+\alpha_k(x_k-x_0)\dots(x_k-x_{k-1})\Rightarrow\boxed{\alpha_k=\frac{y_k-P_{k-1}(x_k)}{\prod_{j=0}^{k-1}(x_k-x_j)}} \]

Prin aceste calcule, am descoperit că polinomul \(P_n\) poate fi scris ca o combinație liniară de polinoame Newton de bază, întrucât fiecare polinom \(Q_k\) poate fi de fapt rescris:

\[ \boxed{Q_k=\alpha_k\pi_k(x)}\,,\,\,\forall k\in\overline{1,n} \]

Utilizarea diferențelor divizate#

Fără a demonstra[36], recurențele din spatele diferențelor divizate se pot colapsa în valorile de forma \(\alpha_k\), \(k\in\overline{1,n}\). Acest lucru ne permite simplificarea notațiilor:

\[\begin{split} P(x)&=[y_0]+[y_0,y_1]\cdot\pi_1(x)+[y_0,y_1,y_2]\cdot\pi_2(x)+\dots+[y_0,\dots,y_n]\cdot\pi_n(x)\\ &=\sum_{k=0}^n[y_0,\dots,y_k]\cdot\pi_k(x) \end{split}\]

Rescrisă dezvoltat, ecuația anterioară ia forma:

\[ \boxed{P(x)=[y_0]+\sum_{k=1}^n\left([y_0,\dots,y_k]\cdot\prod_{i=0}^{k-1}(x-x_i)\right)} \]

Acest polinom este cunoscut sub denumirea de polinom Newton.

Exemplu practic#

Să se calculeze aproximarea funcției \(f(x)=\sqrt{x}\) în \(x=71\) folosindu-se interpolarea cu diferențe divizate, unde valorile utilizate sunt \(x_1=16\), \(x_2=64\) și \(x_3=100\).

Soluție

Să începem cu funcțiile simple: \( \begin{cases} [y_0]=\sqrt{16}&=4\\ [y_1]=\sqrt{64}&=8\\ [y_2]=\sqrt{100}&=10\\ \end{cases} \). Facem acum un pas mai departe și calculăm: \( \begin{cases} [y_0,y_1]=\frac{8-4}{64-16}&=\frac{1}{12}\\ [y_1,y_2]=\frac{10-8}{100-64}&=\frac{1}{18}\\ \end{cases} \). Ajungem și la cea mai complicată funcție intermediară:

\[ [y_0,y_1,y_2]=\frac{\frac{1}{18}-\frac{1}{12}}{100-16}=-\frac{1}{3024} \]

Concluzionăm așadar că polinom de interpolare este:

\[ P(x)=4+\frac{1}{12}\cdot(x-16)-\frac{1}{3024}\cdot(x-16)(x-64) \]

Prin această interpolare, aproximăm \(P(71)\approx{8.456}\). În realitate, ne-am apropiat foarte mult de valoarea reală, căci \(\sqrt{71}\approx8.426\).

Suplimentar, vom ilustra polinomul și eroarea funcției:

Funcția \(\sqrt{x}\) și interpolarea \(P(x)\)

Eroarea interpolării

Putem deci concluziona că am ajuns la un răspuns foarte apropiat de cel real.

Cod ilustrativ#

Implementarea interpolării Newton poate fi făcută sub forma funcției interpolare_newton ce primește parametrii:

M - o matrice de două coloane, unde fiecare linie reprezintă o pereche de puncte;
xp - punctul în care se vrea evaluarea interpolării;

Funcția returnează yp, evaluarea în punctul xp.

function yp = interpolare_newton(M, xp)
    n = size(M)(1);
    dif_div = zeros(n,n);
    for i = [1 : n]
        dif_div(i, i) = M(i, 2);
    end
    for dj = [1 : n-1]
        for i = [1 : n-dj]
            j = i + dj;
            dif_div(i,j) = (dif_div(i+1,j) - dif_div(i,j-1));
            dif_div(i,j) /= (M(j,1) - M(i,1));
        end
    end
    yp = dif_div(1, 1);
    p = 1;
    for k = [2 : n]
        p *= xp - M(k-1,1);
        yp += p * dif_div(1, k);
    end
end

Evident, în exemplul anterior, se consideră cunoscută o funcție dif_div() ce calculează diferențe divizate precum a fost explicat în paragrafele anterioare. Algoritmul poate fi îmbunătățit din punct de vedere al complexității sale temporale, însă această remarcă nu face obiectul subcapitolului nostru.

Licență#

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