Interpolarea cu funcții spline#

În urma metodelor anterioare, putem concluziona că forma interpolării polinomiale devine tot mai absurdă (creșteri și descreșteri bruște) odată cu creșterea numărului de puncte în care se admite cunoscută evaluarea unei funcții oarecare.

Metodele de interpolare practice se folosesc de funcții spline, adică de funcții ce sunt alcătuite din alte funcții mai simple, mai digerabile, adică polinoame mai mici puse cap-la-cap. Acestea formează aproximări mult mai apropiate de funcția inițială, fiind folosite în practică pentru mulțimi oricât de mari de puncte.

Funcții spline#

Pornim de la o mulțime de \(n+1\) abscise distincte ordonate crescător, \((x_n)_{n\in\mathbb{N}}\) și \(x_k\lt x_{k+1}\), \(\forall k\in\overline{0,n-1}\). Totodată, fie \(f:[a,b]\subset\mathbb{R}\rightarrow\mathbb{R}\) astfel încât: \( [a,b]=[x_0,x_1)\,\cup\,[x_1,x_2)\,\cup\,\dots\,\cup\,[x_{n-1},x_n)\,\cup\{x_n\} \).

Dacă îi atribuim funcției \(f\) forma:

unde \(p_0,\dots,p_{n-1}\in\mathbb{R}[x]\) sunt polinoame, atunci \(f\) se consideră funcție spline, adică o funcție formată prin țeserea altor funcții între ele.

Proprietăți#

Fie o funcție spline \(f:[a,b]\subset\mathbb{R}\rightarrow\mathbb{R}\) alcătuită din polinoamele \(p_0,\dots,p_{n-1}\in\mathbb{R}[x]\). Considerăm gradul unei funcții spline definit drept maximul gradelor polinoamelor din care este alcătuită funcția \(f\), adică:

Spunem despre aceeași funcție \(f\) că este de clasă \(\displaystyle C^k\), \(k\in\mathbb{N}\), dacă toate derivatele sale de ordin \(0, 1\dots k\) există și sunt continue.

Această continuitate va trebui cumva comunicată între polinoamele ce alcătuiesc de fapt funcția \(f\), vedem imediat cum.

Funcții spline de clasă \(C^1\), grad \(1\)#

Începem cu cea mai simplă categorie de funcții spline, cele de clasă \(C^1\) și grad \(1\) și încercăm să descoperim de ce s-ar putea să ne mulțumească această formă.

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, pentru a ajunge la o interpolare spline de clasă \(C^1\) și grad \(1\), vom construi segmente geometrice care să unească toate punctele de forma \((x_k,y_k)\) cu \((x_{k+1},y_{k+1})\), \(\forall k\in\overline{0,n-1}\). Prin geometrice, ne referim la faptul că, dacă am fi să le desenăm, acestea ar lua forma unor segmente.

Așadar, în situația fiecărui subinterval \([x_k,x_{k+1})\), se va defini un polinom de forma \(p_k(x)=a_kx + b_k\), unde \(p_k:[x_k,x_{k+1})\rightarrow\mathbb{R}\), \(k\in\overline{0,n-1}\), care respectă următoarele criterii:

• Condițiile de interpolare (\(n+1\) la număr). Ca să putem considera funcția o posibilitate de interpolare a punctelor, este necesar ca:

  \[ \boxed{p_k(x_k)=y_k}\,,\,\,\forall k\in\overline{0,n-1},\text{ iar } \boxed{p_{n-1}(x_n)=y_n} \]

• Condițiile de racordare (\(n-1\) la număr). Deoarece funcția noastră va fi de clasă \(C^1\), înseamnă că aceasta trebuie să fie continuă peste tot, deci:

  \[ \lim_{x\nearrow x_{k+1}}p_k(x)=p_{k+1}(x_{k+1})\,,\,\forall k\in\overline{0,n-2} \]

  Relaxând rigoarea matematică, descoperim că:

  \[ \boxed{p_k(x_{k+1})=p_{k+1}(x_{k+1})}\,,\,\,\forall k\in\overline{0,n-2} \]

Toate aceste condiții însumează \((n+1)+(n-1)=2n\) criterii prin care se poate defini fiecare polinom de gradul întâi \(p_k\in\mathbb{R}[x]\), \(\forall k\in\overline{0,n-1}\).

Prin înlocuiri directe, se obțin următoarele formule care pot fi deduse ușor:

Ca exercițiu, puteți încerca să deduceți singuri aceste formule.

Vizualizare#

Vom arăta această interpolare în paralel cu interpolarea polinomială aferentă aceleiași funcții. Să considerăm faimoasa funcție Runge:

Aceasta este utilizată în textele de specialitate tocmai pentru a arăta natura oscilantă a interpolării polinomiale.

Vom interpola această funcție folosind \(5\), respectiv \(11\) puncte echidistante din intervalul \([-1,1]\) prin două metode diferite - una polinomială (de exemplu, Lagrange) și o interpolare spline cu funcții de clasă \(C^1\) și grad \(1\):

Funcția \(f(x)\) și interpolarea sa polinomială

Funcția \(f(x)\) și interpolarea sa spline de clasă \(C^1\) și grad \(1\) (segmente)

Observăm deci cum această metodă simplă de interpolare spline obține rezultate cu mult superioare metodelor polinomiale discutate până la acest moment.

Cod ilustrativ#

Interpolarea cu funcții spline de clasă \(C^1\) și grad \(1\) poate fi implementată rapid în Octave sub forma funcției interpolare_spline_C1_g1 ce utilizează:

• x - abscisele punctelor din mulțimea \(M\);

• y - ordonatele punctelor din mulțimea \(M\).

Funcția returnează o pereche de doi vectori, [a,b], utilizați pentru a reprezenta polinoame de forma:

function [a, b] = interpolare_spline_C1_g1(x, y)
    n = length(x) - 1;
    a = b = zeros(1, n);
    for k = [1 : n]
        a(k) = (y(k+1) - y(k)) / (x(k+1) - x(k));
        b(k) = (x(k+1) * y(k) - x(k) * y(k+1)) / (x(k+1) - x(k));
    end
end

Funcții spline de clasă \(C^2\), grad \(3\)#

Poate însă nu suntem mulțumiți cu funcțiile ce se limitează la clasa \(C^1\) - poate pentru că nu curg, de exemplu. Să încercăm să rezolvăm acest impediment mergând un pas mai departe

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.

Pentru a se obține o interpolare spline de clasă \(C^2\) și grad \(3\), există mai multe forme pe care le pot lua polinoamele intermediare, însă o vom prefera pe următoarea:

unde \(p_k:[x_k,x_{k+1})\rightarrow\mathbb{R}\), \(p_k\in\mathbb{R}[x]\), și \(a_k,b_k,c_k,d_k\in\mathbb{R}\), \(\forall k\in\overline{0,n-1}\).

În situația fiecărui subinterval \([x_k,x_{k+1})\), se vor respecta următoarele restricții:

• Condițiile de interpolare (\(n+1\) la număr). Deoarece funcția reprezintă interpolarea unor puncte, trebuie ca:

  \[ \boxed{p_k(x_k)=y_k}\,,\,\,\forall k\in\overline{0,n-1},\text{ și } \boxed{p_{n-1}(x_n)=y_n} \]

• Condițiile de racordare (\(3n-3\) la număr). Deoarece funcția noastră va fi de clasă \(C^2\), înseamnă că aceasta trebuie să fie continuă, derivata sa trebuie să fie continuă și derivata sa de ordin 2 trebuie să fie continuă:

  • Continuitatea (\(n-1\) condiții):

    \[ \lim_{x\nearrow x_{k+1}}p_k(x)=p_{k+1}(x_{k+1})\,,\,\,\forall k\in\overline{0,n-2} \]

  • Continuitatea derivatei de ordin \(\bm{1}\) (\(n-1\) condiții):

    \[ \lim_{x\nearrow x_{k+1}}\frac{dp_k}{dx}(x)=\frac{dp_{k+1}}{dx}(x_{k+1})\,,\,\,\forall k\in\overline{0,n-2} \]

  • Continuitatea derivatei de ordin \(\bm{2}\) (\(n-1\) condiții):

    \[ \lim_{x\nearrow x_{k+1}}\frac{d^2p_k}{dx^2}(x)=\frac{d^2p_{k+1}}{dx^2}(x_{k+1})\,,\,\,\forall k\in\overline{0,n-2} \]

Toate acestea pot fi scrise, mai relaxat, astfel:

În total, până la acest moment, am introdus \(4n-2\) ecuații, însă avem \(4n\) necunoscute. Pentru ultimele două ecuații, acceptăm trei tipuri diferite de interpolare:

• Interpolarea naturală. Pentru aceasta, vom considera că derivata de ordin \(2\) în punctele de capăt, \(x_0\) și \(x_n\), este egală cu \(0\):

  \[ \boxed{ \frac{d^2p_0}{dx^2}(x_0)=0 }\,,\,\, \boxed{ \frac{d^2p_{n-1}}{dx^2}(x_n)=0 } \]

• Interpolarea tensionată. Pentru aceasta, vom avea nevoie de informații suplimentare, însă vom considera că derivata de ordin \(1\) în punctele de capăt, \(x_0\) și \(x_n\), este cunoscută, adică:

  \[ \boxed{ \frac{dp_0}{dx}(x_0)=F'(x_0) }\,,\,\, \boxed{ \frac{dp_{n-1}}{dx}(x_n)=F'(x_n) } \]

  unde funcția \(F\) este cunoscută și reprezintă (adesea) funcția de la care s-a pornit când a fost construită mulțimea \(M\), adică \(F(x_k)=y_k\), \(\forall k\in\overline{0,n}\);

• Interpolarea periodică. Nu în ultimul rând, avem posibilitatea construirii unei funcții periodice - această situație introduce \(3\) ecuații în loc de \(2\):

  \[ \boxed{ p_0(x_0)=p_{n-1}(x_n) }\,,\,\, \boxed{ \frac{dp_0}{dx}(x_0)=\frac{dp_{n-1}}{dx}(x_n) }\,,\,\, \boxed{ \frac{d^2p_0}{dx^2}(x_0)=\frac{d^2p_{n-1}}{dx^2}(x_n) } \]

Observați cum, pentru moment, am evitat găsirea unor formule analitice pentru \(a_k\), \(b_k\), \(c_k\) și \(d_k\), motivul fiind matematica anevoioasă.

Pentru cei ce își doresc să înțeleagă însă proveniența formei matriciale ce va fi explicată mai jos, vă invit să lecturați următoarea explicație.

Forma analitică#

Cunoaștem următoarea formă a unuia dintre polinoamele compononente ale interpolării (\(k\in\overline{0,n-1}\)):

Așadar, dacă vrem să evaluăm polinomul \(p_k\) în punctul \(x_k\), vom obține chiar \(\boxed{p_k(x_k)=a_k}\). Continuăm, derivând relația polinomială în funcție de \(x\) (\(f'(x)=\frac{dp_k}{dx}(x)\)):

Din această relație, observăm cum, prin evaluarea în punctul \(x_k\), se obține \(f'(x_k)=b_k\). Facem o ultimă derivare a polinomului \(f'(x)\):

Din ultima relație extragem \(f''(x_k)=2c_k\). Așadar, punând toate aceste informații cap la cap, alături de relațiile de continuitate (racordare) forțate de interpolarea de clasă \(C^2\), concluzionăm că:

În continuare, vom defini \(h_k=x_{k+1}-x_k\), \(\forall k\in\overline{0,n-1}\). Aceste valori pot fi utilizate pentru a simplifica notațiile din cadrul metodei.

Așadar, înlocuind ecuațiile în sistemul anterior conform polinoamelor de forma \(p_k(x)\), \(p_k'(x)\), respectiv \(p_k''(x)\), se obține:

Din ultima egalitate, extragem o relație pentru \(d_k\):

Vom folosi această egalitate astfel încât să înlocuim instanțele lui \(d_k\) în restul sistemului:

Vom prelucra acest sistem, astfel încât să rămânem în prima ecuație doar cu elementele \(a_k\) și \(a_{k+1}\) în partea stângă, respectiv în a doua ecuație doar cu elementele \(b_k\) și \(b_{k+1}\):

Îl vom scoate acum pe \(b_k\) din prima ecuație:

Folosind această relație, este evidentă valoarea lui \(b_{k+1}\):

Cu aceste noi informații, revenim la sistemul de ecuații anterior și înlocuim \(b_k\) și \(b_{k+1}\) în ultima ecuație:

Această ecuație, deși foarte violentă, poate fi îmblânzită prin scrierea sa sub următoarea formă digerabilă:

Alternativ, derulăm cu o unitate mai în spate indicii, adică \(k\rightarrow(k-1)\), ceea ce implică:

Această ultimă ecuație implică existența unui sistem de \((n-1)\) ecuații liniare în care se încearcă aflarea valorilor \(c_0,\dots,c_n\). Cu toate acestea, lipsesc două ecuații - acestea se vor obține în funcție de tipul interpolării; ne vom limita la cea naturală și la cea tensionată:

• Pentru interpolarea naturală, se obține, prin înlocuire, \(c_0=c_n=0\);

• Pentru interpolarea tensionată, cunoaștem \(b_0=f'(x_0)\) și \(b_n=f'(x_n)\), de unde se pot deduce și valorile pentru \(c_0\) și \(c_n\).

Forma matriceală#

În baza tuturor informațiilor anterioare, concluzionăm că valorile \(a_k\), \(b_k\) și \(d_k\) pot fi calculate ușor odată ce sunt cunoscute valorile de forma \(c_k\).

Totodată, am extras câte un sistem tridiagonal pentru fiecare tip de interpolare ce soluționează valorile \(c_k\).

Amintim mai întâi formulele pentru \(a_k\), \(b_k\) și \(d_k\):

Vrem acum să găsim un sistem de forma \(H\bm{c}=\bm{h}\), unde \(H\in\mathbb{R}^{(n+1)\times(n+1)}\) este o matrice tridiagonală, iar \(\bm{c},\bm{h}\in\mathbb{R}^{n+1}\) sunt vectori coloană astfel încât \(\bm{c}\) să conțină toate elementele \(c_0,\dots,c_n\).

Interpolarea naturală#

Am demonstrat anterior următoarea relație:

Folosind-o, alături de faptul că \(c_0=c_n=0\), obținem următorul sistem tridiagonal:

Interpolarea tensionată#

Analog, se obține următorul sistem tridiagonal pentru interpolarea tensionată:

Licență#

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