Elemente de Teoria Aproximarii Functiilor

CUPRINS

Introducere

Lucrarea de fațǎ se înscrie în domeniul aproximǎrii și interpolǎrii funcțiilor, o problemǎ importantǎ a disciplinei Calcul Numeric ce reprezintǎ o ramurǎ a matematicii aplicative ce cautǎ sǎ elaboreze principii, metode și procedee pentru aflarea efectivǎ, cu o precizie determinatǎ a soluțiilor numerice pentru diverse probleme formulate și studiate în matematicǎ. Aceasta pe lângǎ algoritmii de calcul pe care îi stabilește în vederea aflǎrii soluțiilor numerice cǎutate, urmǎrește sǎ elaboreze și o teorie relativǎ la studiul preciziei cu care se pot obține aceste soluții aproximative.

Se poate intampla ca diverse fenomene din jurul nostru care se cerceteazǎ cu ajutorul metodelor matematice sǎ fie caracterizate prin funcții de una sau mai multe variabile reale, însa pentru care nu existǎ posibilitatea determinǎrii lor complete ci doar a aflǎrii valorilor pe anumite valori particulare ale argumentelor. Cunoscând aceste valori se dorește adesea sǎ se aproximeze fie valorile funcțiilor respective pe alte noi puncte din domeniul lor de definiție, fie anumite caracteristici numerice cum ar fi: valorile pe puncte date ale derivatelor de anumite ordine sau integralelor pe anumite domenii.

Pentru aceste scopuri în analiza matematicǎ se elaboreazǎ teoria interpolǎrii, care face apel la calculul cu diferențe divizate și finite.

O problemǎ importantǎ a teoriei aproximǎrii constǎ în aflarea pentru o funcție, dintr-un spațiu , a unei funcții de aproximare simplǎ, definitǎ pe o submulțime

a lui Spațiul este de obicei un spațiu liniar normat, cum ar fi spațiile și sau un alt spațiu de funcții. Distanța de la la se poate mǎsura prin norma în spațiul Pentru funcțiile definite pe un interval compact al axei reale de obicei se alege , spațiul polinoamelor algebrice de grad cel mult O altǎ clasǎ importantǎ studiatǎ și in lucrarea de fațǎ este reprezentatǎ de funcțiile splie polinomiale de gradul pe care o vom nota , unde reprezintǎ o diviziune a intervalului de bazǎ Acestea reprezintǎ o clasǎ de funcții segmentar-polinomiale , cu un anumit grad de netezime, în sensul cǎ graficul lor este form,at din arce de polinoame care se racordeazǎ pe anumite puncte de joncțiune, în așa fel încât acestea sǎ dmitǎ derivate continue pânǎ la un anumit ordin . Un rol important în teoria aproximǎrii îl au teoremele lui K,Weierstrass, care ne asigurǎ cǎ o funcție continuǎ pe un interval compact poate fi aproximatǎ în mod uniform, oricât de bine, prin polinoame algebrice .

Capitolul 1 Elemente de teoria aproximǎrii funcțiilor

1.1 Notații. Spații de funcții

Se poate intampla ca diverse fenomene din jurul nostru care se cerceteazǎ cu ajutorul metodelor matematice sǎ fie caracterizate prin funcții de una sau mai multe variabile reale, însa pentru care nu existǎ posibilitatea determinrii lor complete ci doar a aflǎrii valorilor pe anumite valori particulare ale argumentelor. Cunoscând aceste valori se dorește adesea sǎ se aproximeze fie valorile funcțiilor respective pe alte noi puncte din domeniul lor de definiție, fie anumite caracteristici numerice cum ar fi: valorile pe puncte date ale derivatelor de anumite ordine sau integralelor pe anumite domenii.

Intrucat cea mai mare parte a materialului se referǎ la aproximarea funcțiilor, a unor funcționale liniare definite pe spații de funcții sau la aproximarea soluțiilor unor ecuații operatoriale prin metode generate de procedee de aproximare a funcțiilor, accentul se va pune pe prezentarea cu precǎdere, a spațiilor de funcții.

In general, într-o problemǎ de aproximare intervin trei elemente esențiale:

o multime B ale cǎrei elemente urmeazǎ a fi approximate

o mulțime A a aproximațiilor (cu elementele cǎreia se aproximeazǎ elementele lui B)

un criteriu prin care se selecteazǎ aproximația din A a unui element dat din B

De cele mai multe ori criteiul de selecție are la bazǎ eroarea de aproximare. Nu sunt excluse nici situațiile în care aproximația din A a unui element a lui B se alege astfel încât ea sǎ fie ușor de calculat sau cumularea erorilor de calcul sǎ fie cât mai redusǎ.

Un loc aparte îl ocupǎ problemele în care pentru un element dat se cere determinarea elementului lui A care aproximeazǎ cel mai bine pe f sau elementul din A care aproximeazǎ pe f cu o eroare dinainte datǎ. Dacǎ mulțimea A nu conține nici un element care sǎ aproximeze pe f cu eroarea prescrisǎ, ea poate fi extinsǎ astfel încat condiția pusǎ sa fie satisfacutǎ. De exemplu dacǎ A este mulțimea funcțiilor polinomiale de grad cel mult , , extinderea ei se poate face prin mǎrirea gradului , prin trecerea la mulțimea funcțiilor spline, a funcțiilor raționale, etc.

Apare astfel necesitatea definirii erorii de aproximare a lui prin elementul din A, a comparǎrii mai multor aproximații ale lui f cu scopul de a decide care dintre ele este mai convenabilǎ sau chiar la determinarea celei mai bune aproximații din A a lui f.

In unele cazuri definirea erorii poate depinde si de natura multimii B, de structurile cu care este înzestratǎ. De exemplu dacǎ este un spațiu metric, și , eroarea de aproximare a lui prin se masoarǎ prin distanța . Elementul pentru care oricare ar fi este cea mai bunǎ aproximație a lui f în A.

Amintim cǎ este un spațiu metric dacǎ pe mulțimea B este definiție o aplicație (metricǎ) care satisfac condițiile:

1. si dacǎ și numai dacǎ

2. )

3. pentru orice

Fie B un spațiu liniar peste corpul K (vezi definiția la spațiu liniar)

Daca K=R atunci B este spațiu liniar real, iar dacǎ K=C, B este spațiu liniar complex.

O aplicație a lui B în K se numește funcționalǎ.

O funcționalǎ se numește funcționalǎ aditivǎ dacǎ oricarea ar fi

O funcționalǎ se numește funcționalǎ omogenǎ dacǎ

oricare ar fi si

O funcționalǎ aditivǎ și omogenǎ se numește funcționala liniarǎ.

O funcționala realǎ nenegativǎ definitǎ pe spațiul liniar B cu proprietǎțile:

1. oricare ar fi

2. , oricare ar fi și

se numește seminormǎ pe B.

O seminormǎ p definitǎ pe spațiul liniar B se numește norma dacǎ din

rezultǎ .

Exemple de spații de funcții

1. – spatiul funcțiilor continue pe intervalul având norma

2. -spațiul funcțiilor care admit derivate continue pâna la ordinul inclusiv. In raport cu norma

spațiul este spațiu Banach.

3. spațiul funcțiilor marginite pe [a,b] înzestrat cu norma

este un spațiu Banach.

4. – spațiul funcțiilor continue și mǎrginite

5. spațiul funcțiilor polinomiale de grad cel mult . Avem

1.2 Aproximarea uniformǎ a funcțiilor continue

Aproximarea funcțiilor continue cu ajutorul polinoamelor reprezintǎ un capitol foarte important în cadrul teoriei aproximǎrii funcțiilor. La baza acesteia stǎ o o teoremǎ celebrǎ, enunțatǎ în anul 1885 de Karl Weirstrass (1815-1897), pe care acesta a descoperit-o și publicat-o când avea vârsta înaintatǎ de 70 de ani. Teorema exprimǎ posibilitatea aproximarii polinomiale uniforme a unei funcții , continuǎ pe un interval mǎrginit și închis [a,b] al axei reale.

Sa prezentam în continuare câteva enunțuri echivalente ale acestei teoreme.

Teorema 1.2.1 Dacǎ atunci, oricare este un numar pozitiv , existǎ un polinom algebric , care depinde atât de cât și de , astfel încat sǎ avem

, (1.2.1)

oricare ar fi

Acestei teoreme i se mai spune, de obicei, teorema întâia de aproximare a lui Weierstrass, pentru a o deosebi de a doua teoremǎ a lui Weierstrass , care se referǎ la aproximarea funcțiilor continue periodice, de perioadǎ , prin polinoame trigonometrice.

Aceastǎ teoremǎ, care exprimǎ faptul cǎ mulțimea polinoamelor algebrice este densǎ în spațiul de funcții , se mai poate enunța și în forma echivalentǎ urmǎtoare

Teorema 1.2.2 Orice functie este limita unui șir de polinoame algebrice uniform convergent la funcția pe

O altǎ variantǎ echivalentǎ se regasește în enunțul urmǎtor

Teorema 1.2.3 Pentru orice funcție se poate gǎsi o serie de polinoame care pe intervalul este absolutǎ și uniform continuǎ cǎtre

In virtutea Teoremei 1.2.2 funcției îi corespunde un șir de polinoame algebrice uniform convergent la Formând șirul unde , avem

astfel cǎ seria de termen general este uniform convergentǎ cǎtre pe .

Reciproc, dacǎ se dezvoltǎ în seria de polinoame , uniform convergentǎ cǎtre pe, aceasta înseamna cǎ, oricare ar fi numǎrul pozitiv , existǎ un numar astfel încat sǎ avem pentru , oricare este , unde . Observând cǎ este un polinom, este doveditǎ și echivalentǎ Teoremelor 1.2.2 și 1.2.3.

Interpretarea geometricǎ a teoremei întâia a lui Weierstrass

Având în vedere cǎ inegalitatea (1.2.1) se mai poate exprima prin dubla inegalitate

,

valabilǎ pe intervalul , constatǎm cǎ putem da urmǎtoarea interpretare geometricǎ.

Oricât de “îngustǎ” ar fi banda mǎrginita de curbele de ecuații , , se poate gǎsi un polinom astfel încat graficul lui sǎ fie situate în intregime în acea bandǎ, pe intervalul

1.3 Diferențe divizate

Definiție 1.3.1 Sǎ considerǎm . Diferențele divizate de ordinul p se definesc astfel:

Pentru

Pentru

Pentru

Sǎ prezentǎm câteva proprietǎți importante ale acestei mǎrimi:

i)

ii)este determinantul Vandermonde ce corespunde nodurilorși este determinantul obținut din cel precedent prin inlocuirea coloanelor cu valorile.

Dacǎ atunci

Pentru , mǎrimea

se numește diferența divizatǎ a funcției , relativ la punctele

Definiție 1.3.2 Fie Mǎrimea

se numește diferența divizatǎ de ordinul a funcției pe punctul .

Teorema 1.3.3 Dacǎ atunci

Observație Fie Aplicația

se numește operatorul de diferența divizatǎ.

In cazul diferențelor divizate și-a dovedit utilitatea tabloul cu diferențe divizate:

unde .

Prin urmare elementele tabloului de mai sus se determinǎ dupǎ formula recursivǎ

.

Implementǎri în Matlab

function [A,df]=difdivizate(X,Y)

A=Y;

n=length(X);

for j=2:n

for k=n:-1:j

A(k)=(A(k)-A(k-1))/(X(k)-X(k-j+1))

end

end

x0 =X(1);

df =A(2);

prod =1;

n1=length (A)-1;

for k=2:n1

prod =prod*(x0-X(k));

df =df+prod*A(k+1);

end

Tabelul diferențelor divizate poate fi calculat și folosind urmǎtorul script

x=[0 1 2 4 5];

y=[0 16 48 88 0];

n = length(x)

d= zeros(n,n-1);

for k = 1 : n-1

d (k, 1) = (y(k+1) – y(k))/(x(k+1) – x(k))

end

for j = 2 : n-1

for k = 1 : n – j

d (k, j) = (d(k+1, j – 1) – d(k, j – 1))/(x(k+j) – x(k));

end

end

d

1.5 Diferențe finite

Fie .

Definiție 1.4.1. Pentru ,

se numeste diferența finitǎ de ordinul întâi, cu pasul a funcției pe punctul a.

Observație Aplicația numește operator de diferențǎ finitǎ cu pasul .

Teorema 1.4.2. este un operator liniar.

Demonstrație Dacǎ

=A

=A (

Definiție 1.4.3 Fie 0 si Mǎrimea definitǎ prin ( cu convențiile si , se numește diferența finitǎ cu pasul h de ordinul k al funcției , pe punctul .

Astfel , (f(a

De exemplu, dacǎ avem (

(

Iar prin metoda inducției matematice complete se demonstreaza ca

( (2.5.1.)

Teorema 1.4.4 Dacǎ 0atunci

(

Teorema 1.4.5.Dacǎ 0atunci

(

Demonstratie .Se folosește definiția diferenței finite de ordin superior.

Observație Proprietatea de liniaritate se pǎstreazǎ și pentru operatorul .

In continuare vom prezenta unele formule sau procedee de exprimare a diferențelor finite de ordin superior cu ajutorul valorilor funcției, precum și a valorii unei funcții pe un punct cu ajutorul diferențelor ei finite , formule utile în calculul cu diferențe.

Teorema 1.4.6. Dacǎ , atunci ( (2.5.2.)

unde coeficientii nu depind de .Ca urmare, funcția poate fi particularizatǎ confuncția poate fi particularizatǎ convenabil. Luând și având în vedere relația (2.5.1.) formula (2.5.4.) devine (e

Simplificând cu și aplicând primului membru formula binomului, se obține

Rezultǎ cǎ

Observație Pe baza proprietǎții de simetrie a combinǎrilor, formula (2.5.2.) se mai poate scrie

( (2.5.4.)

Teorema 1.4.7. Dacǎ și , atunci

F (a (2.5.5.)

Demonstrație. Eliminând din relațiile (

valorile , se obține

(2.5.6.)

unde B nu depind de . Luând , din nou, ,exprimand pe sub forma

e si având in vedere (2.5.1.), formula (2.5.6.), dupǎ simplificarea cu , devine Rezultǎ deci cǎ

Teorema 1.4.8. Dacǎ , atunci ( (2.5.7.)

care este adevaratǎ. Intr-adevar,

Adunând și scǎzând termenul , se obține

,

adicǎ formula (2.5.8).

Presupunând formula (2.5.7) adevaratǎ pentru m-1,adicǎ

Folosind liniaritatea operatorului ,resultǎ cǎ

Fǎcând în al doilea termen schimbarea avem

sau, daca grupǎm termenii asemenea,

Cum

teorema este complet demonstratǎ.

Dacǎ atunci ,

relatie care stabilește legatura dintre diferențele divizate și diferențele finite, în cazul nodurilor echidistante.

De mare utilitate pentru calcul este tabloul cu diferențe finite.

Folosind notațiile ,avem:

unde

(2.5.9)

Tabloul cu diferențe oferǎ o metoda simplǎ, algoritmicǎ, de calcul a tuturor diferențelor finite ale unei funcții, relativ la o multime de puncte echidistante date.

Dup cum rezultǎ și din relațiile (2.5.9) elementele fiecǎrei coloane se obțin din elementele coloanei precedente prin simple scǎderi, prin urmare, tabloul se completeazǎ pe coloane.

Relațiile (2.5.9) permit și o programare simplǎ pentru generarea tabloului cu diferențele finite la calculator.

In multe aplicații concrete intervin numai elementele de pe latura superioarǎ sau inferioarǎ a acestui tablou triunghiular. In astfel de cazuri este necesarǎ reținerea numai a unei pǎrți din elementele tabloului, ceea ce permite, având în vedere cǎ odatǎ elementele coloanei k fiind calculate, elementele coloanei k-1 nu mai intervin în calculele urmǎtoare, o esențialǎ economie de memorie.

Mai concret, dacǎ trebuie reținute numai diferențele de pe latura superioara,

atunci pe lânga zona de memorie care înregistreazǎ datele inițiale mai sunt necesare m locații care în diferențele deja calculate,de exemplu pentru precum și restul elementelor coloanei i, adicǎ , j = 1,…,m-i.

Capitolul 2 Aproximarea prin interpolare polinomialǎ

2.1 Preliminarii

Pentru o funcție problema aproximării ei printr-un polinom se pune fie când este dificil de evaluat f, fie când nu se cunoaște expresia analitică a lui f ci doar valorile ei în anumite puncte , , obținute în general ca urmare a unor măsurări și prezentate într-un tabel. Mulțimea de puncte , , cu proprietatea , o vom nota prin și o vom numi diviziune a intervalului . Punctele , vor fi numite nodurile diviziunii. Problema interpolǎrii funcției în nodurile , , constǎ în determinarea unei funcții dintr-o clasǎ de funcții cunoscute, cu proprietatea Pusǎ sub aceastǎ formǎ generalǎ problema poate sǎ nu aibǎ soluție sau sǎ aibǎ o infinitate de soluții.

Cea mai utilizatǎ clasǎ de funcții de interpolare este clasa polinoamelor, datoritǎ ușurinței cu care se integreazǎ și se deriveazǎ.

2.2 Interpolare Lagrange

Teorema 2.2.1. Fie f : [a,b]R si , (n+1) noduri din intervalul [a,b]. Atunci existǎ un polinom unic , de gradul , care interpoleazǎ funcția în nodurile .Acest polinom se numește polinomul de interpolare al lui Lagrange.

Demonstrație

Cǎutam un polinom sub forma urmǎtoare: ,

unde sunt polinoame de gradul ce urmeazǎ sǎ fie determinate. Deoarece dorim ca vom pune condițiile:

Deoarece pentru rezultǎ cǎ admite rǎdǎcinile

Așadar,

Cum , rezultǎ

In concluzie avem

(2.2.1)

unde

(2.2.2)

Evident polinomul (2.2.1) are gradul n si are proprietatea Fie un alt polinom de gradul n cu proprietatea și sǎ fie . Deoarece grad R rezultǎ cǎ R este polinom identic nul, deci ca

Exemplu. Fie nodurile . Atunci

Efectuând calculele obținem

In continuare vom nota eroarea în fiecare punct cu

(2.2.3)

Evident , . Introducem deasemenea notația :

(2.2.4)

Teorema 2.2.2 Dacǎ atunci pentru orice existǎ astfel incât

(2.2.5)

Demonstrație

Considerǎm funcția auxiliarǎ

.

Observǎm cǎ g se anuleazǎ în (n+2) puncte distincte Din teorema lui Rolle rezultǎ cǎ existǎ

Cum

rezultǎ

Corolar. Daca existǎ astfel încât atunci :

Observație

Pentru a minimiza eroarea la interpolarea Lagrange pentru o funcție

definitǎ pe intervalul , nodurile de interpolare se aleg astfel

.

Interpolarea cu aceste noduri se numește interpolare Cebâșev.

Intr-adevăr, daca este polinomul de interpolare pentru si este polinomul de interpolare pentru g, atunci este polinomul de interpolare pentru f+g, si deci

In continuare vom presupune ca nodurile sunt echidistante, deci ca , unde

(2.2.6)

Considerâm de asemenea schimbarea de variabilǎ

(2.2.7)

Înlocuind (2.2.6) si (2.2.7) in (2.2.2) obținem:

Folosind notația:

(2.2.8)

obținem

Obținem astfel expresia polinomului lui Lagrange pentru noduri echidistante

(2.2.9)

Eroarea devine:

(2.2.10)

In continuare considerǎm un șir de diviziuni ale intervalului [a,b] cu

Notǎm cu polinomul lui Lagrange care interpolează funcția f în nodurile . Daca n este mai mare, coincide cu f într-un număr mare de noduri, deci ne așteptam ca eroarea sǎ fie micǎ, eventual ca

Ajungem astfel la următoarea întrebare:

In ce condiții șirul de polinoame converge punctual (eventual uniform) la funcția f pe intervalul {a,b}?

In anul 1912 S. N. Bernstein a arătat ca pentru funcția , daca alegem nodurile echidistante , atunci daca

S-ar putea crede ca acest lucru se datorează faptului ca funcția-modul nu este derivabila in origine. Următorul exemplu dat de C. Runge in 1901 aratǎ cǎ existǎ funcții indefinit derivabile pentru care nu converge la

Fie

Evident Fie nodurile echidistante

.

Se poate arata ca dacǎ si dacǎ , unde este o rădăcina a ecuației:

In anul 1914 S. N. Bernstein a arătat cǎ pentru orice sistem de noduri , din intervalul , dat dinainte, existǎ o funcție continuǎ , astfel încât șirul polinoamelor lui Lagrange care interpolează funcția f în aceste noduri nu converge uniform la f pe .

Exista totuși si situații când convergenta are loc. Se poate demonstra următoarea teorema:

Teorema 2.2.3. Dacǎ și se dezvoltǎ în serie Taylor pe R, atunci pentru orice sistem de noduri distincte și echidistante din [a,b], șirul polinoamelor care interpolează funcția f in aceste noduri converge uniform la pe [a,b].

Se pune întrebarea daca interpolarea cu polinoame Lagrange este utilǎ în practica din moment ce, cum am văzut, în general șirul polinoamelor de interpolare nu converge la .

Răspunsul este cǎ interpolarea cu polinoame Lagrange este utilǎ. Se constatǎ în practicǎ faptul cǎ pentru un punct , eroarea scade pânǎ la un punct pe măsura ce n crește, și deci, pentru n relativ mic aproximează acceptabil valoarea. Pentru valori mari ale lui n, interpolarea Lagrange nu este recomandatǎ.

Din cele prezentate pânǎ acum, rezultǎ cǎ șirul polinoamelor de interpolare asociate unei funcții continue nu converge uniform, în mod necesar, la aceasta funcție. Se pune întrebarea dacǎ o funcție continuǎ poate fi aproximatǎ uniform, cu polinoame. Răspunsul a fost dat de K. Weierstrass in anul 1885.(vezi sectiunea 1.2)

Teorema 2.2.4. Fie continuǎ. Atunci pentru orice , existǎ un polinom , astfel încât

.

Evident, dacǎ luǎm rezultǎ cǎ existǎ un șir de polinoame , care converge uniform pe sunt, în raport cu funcțiile continue pe , in aceeași relație ca numerele raționale Q fațǎ de numerele raționale R.

Teorema lui Weierstrass este extrem de importantǎ în analiza matematicǎ, în general, și în analiza numericǎ, în special. Dintre numeroasele demonstrații date acestei teoreme, cea mai cunoscutǎ este cea datǎ de S. N. Bernstein, în anul 1912. Bernstein a arătat cum se poate construi șirul de polinoame care aproximează funcția f si anume:

Acest șir de polinoame, care se numesc polinoame Bernstein, au proprietatea , pe [0,1] la [a,b] se face cu ușurința printr-o schimbare de variabila. Evident, polinoamele Bernstein nu sunt polinoame de interpolare. Din păcate, convergenta șirului câtre este destul de înceată, și din aceastǎ cauza, în practicǎ, polinoamele Bernstein nu se folosesc la aproximarea directǎ a funcțiilor. Teorema lui Weierstrass este importantǎ prin implicațiile sale teoretice, dar și practice, de exemplu, la integrarea numericǎ.

Exemplu. Fie funcția si nodurile 0.4; 0.5; 0.7; 0.8. Evaluǎm eroarea în punctul x=0.6.

Rezultǎ

,

acest numǎr fiind doar un majorant al erorii.

Dacǎ folosim urmǎtoarele valori de noduri

și calculǎm polinomul lui Lagrange obținem :

Pe de alt parte ln(0.6)=-0.510826. Rezultǎ cǎ , ceea ce confirmǎ afirmația de mai sus.

Observatie. Dacǎ f=Q este un polinom de grad cel mult n, atunci E (f; x)=0, oricare [a,b].

Afirmația rezultǎ din Teorema 2.2.2 deoarece, în acest caz

Implementǎri în Matlab

m=5;

x=0.2:1.8/m:2;

y=1./x.^2;

xmin = min(x);

xmax = max(x);

XX=xmin:0.1:xmax;

YY=1./XX.^2;

plot (XX,YY,'r')

for k=1:length(XX)

[z l]=valLmf(XX(k),m,x,y);

ZZ(k)=z

end

hold on

plot (x,y,'om')

hold on

plot (XX,ZZ,'b')

title ('functia si polinomul Lagrange')

legend ('rosu- Functia ','magenta- valorile functiei in X dat','albastru – polinomul Lagrange')

function [v,l] = valLmf(x,m,X,Y)

% valLmf calculeaza valoarea poliomului de interpolare Lagrange intr-un punct x;

% m- gradul polinomului Lagrange

% X- vectorul nodurilor de interpolare

% Y- vectorul valorilor fuctiei pe noduri

l=zeros(0,m);

for k=0:m

p1=1;

p2=1;

for j=0:m

if (j~=k)

p1=p1*(x-X(j+1));

p2=p2*(X(k+1)-X(j+1));

end

l(k+1)=p1/p2;

end

end

v=sum(l.*Y)

Fig 2.1

Presupunând cǎ , vom fi tentați sǎ credem cǎ numarul al nodurilor crește indefinite, astfel ca distanța între douǎ noduri consecutive oarecare sǎ tindǎ la zero, atunci polinomul va tinde cǎtre funcția Situația nu este însǎ mereu așa. Astfel este binecunoscut contraexemplul lui C. Runge dat inca din 1901 pentru funcția , oricare este din intervalul . Astfel acesta a dovedit cǎ polinomul de interpolare a lui Lagrange corespunzǎtor , care folosește noduri ce impart în pǎrți egale intervalul , converge cǎtre numai in intervalul , unde în timp ce în intervalele diverge, deși este o funcție analiticǎ.

2.3 Interpolare Newton

2.3.1 Polinomul Newton de interpolare de speța I (ascendent)

Fie o funcție, fie , n + 1 puncte distincte din intervalul [a, b], și pentru orice . Punctele se presupun echidistante adică

fiind pasul rețelei. Pentru determinarea polinomului de grad mai mic sau egal cu ce satisface relațiile:

Se pleacă de la un polinom de forma:

. Coeficienții se determinǎ plecând de la relațiile utilizând noțiunile de diferențe finite și putere generalizatǎ. Câteva din proprietǎțile diferențelor finite se regǎsesc in sectiunea 1.4. Produsul se numește putere generalizatǎ a lui Pentru determinarea coeficientului al polinomului de interpolare se utilizează diferența finită la stanga de ordinul a lui în . Se obține

Aproximarea funcției prin polinomul Newton ascendent este avantajoasa pentru valorile din vecinatatea lui , eroarea fiind micǎ pentru un in modul, mic.

2.3.2 Polinomul Newton de interpolare de speța II (descendent)

Fie o funcție, fie , n + 1 puncte ehidistante din intervalul [a, b], și pentru orice .Ca și in cazul polinomului Newton ascendent pentru determinarea polinomului Newton descendent se pleacǎ de la un polinom de forma:

Pentru exprimarea coeficienților se utilizează diferențele finite la stânga lui .

Expresia

se numește diferență finită (la stânga) de ordinul întâi. Diferențele finite la stânga de ordin superior se definesc cu ajutorul relației:

Se observă că

In aplicații, pentru calculul diferențelor finite la stânga se utilizează tabelul diagonal următor:

. …………..

……………….

…………………………

………………………………

. .

. .

.

Pentru determinarea coeficientului al polinomului de interpolare se utilizează diferența finită la stanga de ordinul a lui în . Se obține:

Expresia polinomului de interpolare de mai sus poartă denumirea de polinom Newton de speța II (sau descendent). Expresia polinomului Newton descendent se poate pune sub o formă mai comodă pentru aplicații, dacă se face schimbarea de variabilă

Se obține

Aproximarea funcției prin polinomul Newton descendent este avantajoasă pentru valorile din vecinătatea lui , eroarea fiind mică pentru în modul mic.

Metoda de interpolare Newton are urmǎtoarele avantaje :

realizeazǎ un compromise maxim între efortul de construcție și cel de evaluare

algoritmul este stabil din punct de vedere numeric, având erori numerice acceptabile ale rezultatelor;

permite mǎrirea gradului polinomului de interpolare prin inserarea unui nod nou în rețeaua de interpolare, cu reutilizarea coeficienților de la gradul anterior, care nu se modificǎ, deci și un effort minim de calcul;

are posibiltatea de a controla eroarea de interpolare;

coeficienții polinomului de interpolare Newton reprezintǎ diferențele divizate de interpolare, ceea ce faciliteazǎ calculul numeric al polinomului;

Algoritm pentru determinarea polinomului de interpolare Newton (ascendent)

Date de intrare

– valorile nodurilor de interpolare;

– pasul rețelei de noduri;

-valorile funcției pe noduri;

– 26ewton26l în care se calculeazǎ polinomul;

Date de ieșire

– valoarea polinomului de interpolare Newton ascendent

Algoritm pentru determinarea polinomului Newton descendent

Date de intrare

– valorile nodurilor de interpolare;

– pasul rețelei de noduri;

-valorile funcției pe noduri;

– 26ewton26l în care se calculeazǎ polinomul;

Date de ieșire

– valoarea polinomului de interpolare Newton descendent

Implementǎri în Matlab

function fp = newton_interpolation(x,y,p)

% x- vectorul nodurilor de interpolare;

% y- vectorul valorilor functiei pe noduri

% p- 26ewton26l in care se calculeaza polinomul lui Newton cu diferente divizate

n = length(x);

a(1) = y(1);

for k = 1 : n — 1

d(k, 1) = (y(k+1) — y(k))/(x(k+1) — x(k));

end

for j = 2 : n — 1

for k = 1 : n — j

d(k, j) = (d(k+1, j — 1) — d(k, j — 1))/(x(k+j) — x(k)); % calculul diferentelor divizate

end

end

d %afisarea tabelului diferentelor divizate

for j = 2 : n

a(j) = d(1, j-1);

end

Df(1) = 1;

c(1) = a(1);

for j = 2 : n

Df(j)=(p — x(j-1)) .* Df(j-1);

c(j) = a(j) .* Df(j);

end

fp=sum(c);

Problema

Sǎ se evalueze polinomul Newton cu diferențe divizate pentru datele de mai jos

Vom introduce în fereastra de comandǎ a mediului Matlab parametrii de mai sus pentru care vom apela apoi functia “newton_interpolation.m”

>> x = [1 2 3 5 7];

>> y = [-1 1 5 13 52];

>> p=3;

>> newton_interpolation(x,y,p)

d =

2.0000 1.0000 -0.2500 0.1708

4.0000 0 0.7750 0

4.0000 3.8750 0 0

19.5000 0 0 0

ans =

5

Să reprezentăm graficul polinomului de interpolare Newton, în cazul funcției , definita pe intervalul [0.2 ,2]. Pentru aceasta va trebui să calculăm valoare polinomului Newton pe o mulțime de puncte din acest interval. Vom folosi un cod de program Matlab care va apela funcția newton_interpolation(), definită mai sus :

x=0.2 :.8 :2 ;% nodurile

y=1./x.^2 ; % valorile functiei pe noduri

n=length(x)-1;

xmin=min(x);

xmax=max(x);

XX=xmin:0.1:xmax;% punctele in care se avalueaza polinomul Newton de interpolare

YY=1./XX.^2;% valorile functiei in aceste puncte

plot(XX,YY,’r’);

for k=1:length(XX)

[nl]=28ewton_interpolation(x,y,XX(k))% valorile polinomului Newton in XX

NP(k)=nl

end

plot(XX,NP,’b’)% graficul polinomului Newton de interpolare

In urma rularii codului de program de mai sus salvat cu numele plot_newton.m se obțin valorile polinomului si graficul acestuia

NP =

Columns 1 through 7

25.0000 20.7253 16.8148 13.2685 10.0864 7.2685 4.8148

Columns 8 through 14

2.7253 1.0000 -0.3611 -1.3580 -1.9907 -2.2593 -2.1636

Columns 15 through 17

-1.7037 -0.8796 0.3086

Fig. 2.2

Adăugând codului de mai sus instrucțiunile

hold on

plot(XX,YY,’—‘);

legend(‘polinomul Newton’,’– functia 1/x^2’);

putem obține în același 29ewton de axe și reprezentarea funcției pe aceleași puncte, putând face astfel o comparație a aproximării funcției cu acest polinom. Obținem:

Fig. 2.3

Așa cum aminteam și la interpolarea Lagrange, procesul de interpolare nu întotdeauna produce șiruri de polinoame ce converg la functia interpolate când gradul polinomului tinde la infinit. Vom exemplifica în Matlab binecunoscutul exemplu al lui Runge (vezi pag. 9 )

m=input(‘dati numarul interpolantilor=’);

for k=1:m

n=input(‘dati gradul pol de interpolare=’);

hold on

x=linspace (-5, 5,n+1);

y=1./(1+x.*x);

z=linspace (-5.5,5.5);

t=1./(1+z.^2);

h1_line=plot(z,t,’-.’);

set(h1_line,’LineWidth’,1.25)

t=30ewton_interpolation(x,y,z);

h2_line=plot(z,t,’r’);

set(h2_line,’LineWidth’,1.3,’Color’,[0 0 1])

axis ([-5.5 5.5 -.5 1])

title(sprintf(‘exemplu de divergenta (n=%2.0f)’,n))

xlabel(‘x’);

ylabel(‘y’);

legend (‘y=1/(1+x^2)’,’interpolant’)

hold off

end

Fig. 2.4

Observație Evident polinoamele Lagrange și Newton diferǎ numai prin formǎ, pentru aceeași rețea de puncte, restul fiind același în ambele forme. Ca volum de calcul în schimb se preferǎ polinomul lui Newton care necesitǎ un numǎr de operații aritmetice de fațǎ de pentru polinomul Lagrange. Necesarul de memorie este același pentru ambii algoritmi. Așa cum s-a vǎzut în cazul polinomului de interpolare Newton este nevoie de o matrice suplimentarǎ în cares e memoreaza tabelul diferențelor divizate, însǎ din table se folosesc doar N coeficienți , existând posibilitatea refolosirii celulelor de memorie. Aceasta este partea cea mai costisitoare a algoritmului Newton.

2.4 Interpolare Hermite

Fie , , pentru si

astfel încât . Problema de interpolare Hermite consta în determinarea polinomului P de grad minim care verificǎ

Fie

Expresia polinomului de interpolare Hermite este:

unde sunt polinoamele fundamentale de interpolare Hermite date de

Vom da un procedeu de obținere a polinomului de interpolare Hermite cu noduri duble. Expresia sa este

unde

iar

Aplicarea expresiilor de mai sus este dificilǎ datoritǎ apariției derivatelor. Vom utiliza o metodǎ derivata din formula de interpolare a lui Newton.

Presupunem cǎ se dau nodurile si valorile ,, . Definim secvența prin

Construim acum tabela diferențelor divizate utilizând nodurile , . Deoarece pentru orice , este o diferența divizatǎ cu un nod dublu și este egalǎ cu , deci în locul diferențelor divizate de ordinul 1 vom utiliza valorile derivatelor . Restul diferențelor divizate se obțin în maniera obișnuitǎ, așa cum se aratǎ în tabelul de mai jos. Se obține astfel un algoritm care poate fi folosit și pentru alte interpolari osculatoare. Se pare cǎ metoda este datoratǎ lui Powell.

Algoritmul de interpolare Hermite. Obține coeficienții polinomului Hermite , cu noduri duble , , corespunzǎtor funcției .

Intrare. ,

Ieșire. Numerele unde

P1. Pentru executǎ pașii 2 și 3.

P2.

P3. Dacǎ atunci

9

P4. Pentru executǎ

Pentru j executǎ

P5. Extrage . Stop.

Exemplu numeric. Se considerǎ urmǎtoarele date de intrare

Cu ajutorul diferențelor divizate din tabel (unde datele de intrare sunt subliniate) se obține:

(1.5)=0.6200860 + (1.5-1.3) (-0.5220232) + (1.5-1.3)2(-0.897427) + (1.5-1.3)2(1.5-1.6)(0.0663657) + (1.5-1.3)2(1.5-1.6)2(0.0026663) + (1.5-1.3)2(1.5-1.6)2(1.5-1.9)(-0.0027738)=

=0.5118277

Implementǎri în Matlab

function [val]=interpolare_Hermite(x,x1,f,f1)

%functie care calculeaza valoarea unei functii in punctul x

%daca se cunosc:

% x1-nodurile

% f-valorile functiei in nodurile x1

% f1-valorile derivatei intai in nodurile x1

%rezultatul este pus in val

m=length(x1);

for i=1:m

z(2*i-1)=x1(i);

z(2*i)=x1(i);

Q(2*i-1,1)=f(i);

Q(2*i,1)=f(i);

Q(2*i,2)=f1(i);

if (i~=1)

Q(2*i-1,2)=(Q(2*i-1,1)-Q(2*i-2,1))/(z(2*i-1)-z(2*i-2));

end

end

for i=3:2*m,

for j=3:i,

Q(i,j)=(Q(i,j-1)-Q(i-1,j-1))/(z(i)-z(i-j+1));

end

end

s=Q(1,1);

p1=1;

for i=2:2*m

if (mod(i,2)==1)

p1=p1*(x-x1((i-1)/2));

else

p1=p1*(x-x1(i/2));

end

s=s+p1*Q(i,i);

end

val=s;

function grafic_Hermite(x,v,x1,f,f1)

%functie care reprezinta pe acelasi grafic f si polinomul de interpolare Hermite

%parametrii

% x – intervalul pe care se face graficul

% v – valorea reala a functiei pe intervalul dat de x

% x1 – nodurile in care se cunoaste valoarea functiei

% f – valorea functiei in nodurile x1

% f1 – valorea derivatei de ordinul I in nodurille date de x1

for i=1:length(x)

val(i)=interpolare_Hermite(x(i),x1,f,f1);

end

clf

hold on

plot(x,v,'r-','LineWidth',5); %rosu -> graficul real al functiei

plot(x,val,'g-','LineWidth',2);%verde -> graficul cu interpolarea Hermite

plot(x1,f,'bo','LineWidth',6); %albastru -> nodurile in care se cunoaste valoarea functiei

hold off

Sa dam in continuare un exemplu de reprezentare al unei functiei si a polinomului Hermite asociat

x1=1:0.5:10;

f=sin(x1);

f1=cos(x1);

x=1.2:0.3:10.1;

v=sin(x);

clf

grafic_Hermite(x,v,x1,f,f1);

clear

Fig. 2.5

Problemǎ

O mașinǎ ce se deplaseazǎ pe un drum drept este cronometratǎ în mai multe puncte. Datele obținute din observații sunt date în tabelul de mai jos. Utilizați interpolarea Hermite pentru a prevedea poziția mașinii și viteza sa la momentul t=10.

Soluție Matlab

Pentru rezolvarea problemei se va scrie urmatorul fisier Matlab, “Problema.m”

clc

x=10;

x1=[0 3 5 8 13];

f=[0 225 383 623 993];

f1=[75 77 80 74 72];

rez = interpolare_Hermite(x,x1,f,f1);

fprintf ('Distanta=%3.15f\n',rez);

fprintf ('Viteza=%3.13f\n',rez/x);

La apelarea în fereastra de comandǎ a fișierului cu numele “Problema.m” se va obține rezultatul

Distanta=742.502839098770890

Viteza=74.[anonimizat]

Dacǎ se dorește interpolarea datelor cu ajutorul funcțiilor cubice de tip Hermite vom folosi în Matlab, funcția pchip.

Sintaxa acestei funcții este urmǎtoarea:

yi=pchip(x,y,xi),

unde

– x și y reprezintǎ vectorii ce conțin abscisele și ordonatele datelor de interpolare

– xi- reprezintǎ vectorul în care se dorește interpolarea, cu pas mai fin

-yi- este vectorul returnat de funcție asociat datelor xi

Capitolul 3 Interpolare spline

Considerând un numar de N+1 puncte prin care se impune sǎ treacǎ graficul unui polinom al cǎrui coeficienți trebuie determinați sau a cǎrui valoare trebuie calculatǎ a rezultat cǎ gradul acestuia este N în cazul interpolarii Lagrange respective 2N+1 în cazul interpolarii Hermite. Dupǎ cum s-a vǎzut numǎrul de operații necesar pentru calcularea valorii unui astfel de polinom este destul de mare, crescând cu creșterea numǎrului de noduri. În plus, deseori polinomul de interpolare între noduri, are abateri fațǎ de funcția care se interpoleazǎ.

Acestea sunt motivele pentru care în unele situații se preferǎ sǎ se lucreze cu funcții diferite pe diferite intervale, cazul cel mai simplu fiind cel al interpolǎrii liniare, când graficul funcției este înlocuit cu o linie poligonalǎ ce trece prin cele N+1 puncte.

Istoria funcțiilor spline începe în anul 1946 prin publicarea lucrǎrilor aparținând lui I.J. Schoenberg, în care se introduce și denumirea de funcție spline, preluatǎ de la liniarul elastic (termenul din limba englezǎ) utilizat pentru trasarea unor curbe ce trebuie sǎ treacǎ prin mai multe puncte de pe desen.

3.1 Interpolare cu funcții spline de ordinul întâi

Interpolare cu funcții spline de ordinul întâi

Fie și o diviziune a intervalului

Se numește funcție spline de ordinul întâi, corespunzǎtoare diviziunii o funcție care satisface condițiile:

i) este polinom de gradul întâi pentru

ii) .

Se noteazǎ cu mulțimea funcțiilor spline polinomiale de ordinul întâi, corespunzǎtoare diviziunii

Pentru orice existǎ și sunt unice astfel încât:

unde pentru .

Mai mult, avem

2) Funcțiile constituie o bazǎ în .

Funcțiile definite prin:

constituie o bază în și oricare ar fi avem

(3.1.1)

Problema de interpolare cu spline-uri de ordin întâi constă în:

se dau care pot fi valorile unei funcții continue

se caută

Dacă

În cazul general, problema de interpolare are soluție unică dacă și numai dacă e îndeplinită una dintre condițiile:

;

;

3.2 Interpolare cu funcții spline cubice

Fie o diviziune a intervalului si .

Se numește funcție spline cubică, corespunzătoare diviziunii o funcție care satisface condițiile:

Se noteaza cu mulțimea funcțiilor spline cubice, corespunzătoare diviziunii .

Funcția spline cubică asociată sistemului de puncte este o funcție cubică corespunzătoare diviziunii cu proprietatea

.

În condițiile formulate mai sus, , funcția s se reprezintă în mod unic sub forma:

, (3.2.1)

unde

(3.2.2)

iar sunt soluțiile sistemului:

(3.2.3)

Pentru determinarea unică a spline-ului cubic, sistemul (3.2.3) trebuie completat cu condiții la capete. De exemplu putem presupune că sunt cunoscute valorile pentru

(3.2.4)

sau că sunt cunoscute valorile pentru :

(3.2.5)

Relatiile (16) sunt echivalente cu:

(3.2.6)

Remarcă. Din (3.2.3) și (3.2.5) sau din (3.2.3) și (3.2.5) se pot determina , apoi din (3.2.2) se determină pentru și deci, conform cu (3.2.1), funcția este cunoscută pe [a,b]. Pentru de două ori derivabilă determinăm funcția spline cubică de interpolare asociată diviziunii și funcției, adică, o funcție spline cubică cu proprietătile și

Exemplu: Să se determine valorile funcției f în punctul folosind interpolarea spline cubică, știind că:

și că valorile lui f’ în punctele și sunt: .

R. Aplicăm teorema 1 și vom găsi funcția spline cubică s ce interpolează funcția f , dacă determinăm coeficienții .

Coeficienții se determină prin rezolvarea sistemului liniar AM=b, unde

și Obținem că .

Punctul . Scriem funcția de interpolare s pe acest interval:

Calculăm valoarea funcției s în punctul dat:

Deci, valoarea aproximativă a lui f în punctul este 0.60875.

In continuare vom defini funcțiile B-spline cubice și vom arăta că orice funcție spline cubică care interpoleză funcția in nodurile se reprezintă ca o combinație liniară unică de funcții B-spline cubice. Fie o diviziune a intervalului cu noduri echidistante (unde ). Asociem acestei diviziuni, diviziunea care are in plus șase noduri auxiliare, de asemenea echidistante. :.

Definim orice , funcția B-spline cubică astfel:

Graficul funcției arată astfel:

Din definiția funcțiilor rezultă că

Se verifică ușor că si dacă . Se observă de asemenea că , deci este o funcție spline cubică.

Propoziția 3.2.1. Funcțiile sunt liniar independente pe R.

Demonstrație. Fie combinația liniară

Daca dăm lui x succesiv valorile obținem sistemul:

Deoarece matricea sistemului este strict diagonal dominantă, deci nesingulară, rezultă că sistemul admite numai soluția banală.

În continuare notăm cu spațiul liniar generat de funcțiile , .

Teorema 3.2.2. Există o funcție unică B care interpolează funcția in nodurile

Demonstrație. Fie

Punând condiția ca B să interpoleze funcția in nodurile rezultă

Se obține sistemul:

(3.2.7)

Sistemul (3.2.7 ) are (n+3) ecuații liniare și (n+3) necunoscute

Eliminând necunoscuta din primele două ecuații și necunoscuta din ultimele două ecuații, obținem sistemul echivalent:

(3.2.8)

Matricea coeficienților sistemului este strict diagonal dominantă,deci sistemul (3.2.8) are soluție unică. Așadar, sistemul (3.2.7) are soluție unică.

Înlocuind această soluție in (3.2.6) obținem funcția B căutată, care evident este unică.

Teorema 3.2.3 Fie o funcție spline cubică care interpolează funcția în nodurile . Atunci există unic determinate, astfel incât

Demonstrație:

Din teorema 3.2.2 rezultă că există astfel incât B interpolează funcția in nodurile . Funcția este de clasă pe și este polinomială de gradul trei pe porțiuni. Rezultă că restricția lui B la intervalul [a,b] este o funcție spline cubică, care interpolează in nodurile . Cum asemenea funcție este unică (conform Teoremei 3.2.2 ) rezultă

Unicitatea coeficienților este asigurată de Teorema 3.2.2.

Pachetul de programe MATLAB conține funcția spline care permite interpolarea unei funcții în punctele , , finit, printr-o funcție spline cubică , dacă se cunosc valorile ale funcției in nodurile . Secvența de apelare este: .

Exemplu. Să se determine valorile funcției în punctele folosind interpolarea spline cubică, știind că:

utilizând MATLAB.

În mediul MATLAB se scriu comenzile:

% interpolarea cu functii spline cubice folosind pachetul de programe Matlab

function [x,y,xi,yi]=cub

% Nodurile

x=[0,pi/6,pi/4,pi/3,pi/2];

% Valorile functiei in noduri

y=[0,1/2,1/2^(1/2),3^(1/2)/2,1];

% Valorile in care se interpoleaza functia

xi=[pi/12,pi/8,pi/5];

% Apelarea functiei Matlab spline care face interpolarea

yi=spline(x,y,xi);

Funcția considerată este și putem compara valorile de interpolare cu cele "exacte":

zi=sin(xi);

Pentru reprezentarea grafică a funcției interpolate se poate folosi următoarea secvența MATLAB.

% Reprezentarea grafica a functiei interpolate

plot(x,y,xi,yi,'*',xi,zi,'o');

axis([0,pi/2,0,1.2]); % se stabilesc intervalele de reprezentare pe

axe

title('Interpolarea cu spline cubice');

xlabel('Unghiul');

ylabel('Valorile functiei');

grid

Fig 3.1

Curba obținutǎ prin interpolare spline cubicǎ este o curbǎ netedǎ, definitǎ de un set de polinoame de gradul trei. Curba pe porțiuni, între fiecare pereche de puncte este un polinom de gradul trei calculat astfel încât sǎ conducǎ la tranziții netede de la un polinom de gradul trei la altul. De pildǎ șase puncte sunt conectate prin intermediul a cinci curbe diferite de grad trei, ce constituie o funcție netedǎ între toate cele șase puncte.

În Matlab curba de interpolare spline cubicǎ se calculeazǎ cu funcția spline care se apeleazǎ cu sintaxa:

yi=spline(x,y,xi)

unde

x si y sunt vectori care conțin abscisele si ordonatele datelor inițiale

xi este un vector care conține noile abscise, de regulǎ cu pas mai fin

yi este vectorul returnat, asociat lui

Comparând funcția spline cu funcția pchip, folositǎ la interpolarea cu polinoame Hermite, se desprind urmǎtoarele:

polinoamele spline sunt construite în același mod ca și polinoamele Hermite de interpolare;

spline alege pantele diferite, cu alte cuvinte fǎcând ca continuǎ, dar determinǎ un rezultat incert, pentru continuǎ;

spline determinǎ multe rezultate corecte, dacǎ datele conțin o funcție incertǎ;

pchip nu depǎșește cea mai micǎ oscilație, dacǎ data este incertǎ;

pchip este mai puțin costisitoare, ca timp de calcul, de reglat;

cele douǎ funcții sunt la fel de costisitoare de evaluat;

>> x =-3:3;

>> y =[-1 -1 -1 0 1 1 1];

>> t =-3:.01:3;

>> p=pchip(x,y,t);

>> s=spline(x,y,t);

>> plot (x,y,'o',t,p,'-',t,s,'-.');

>> legend ('date','p-pchip','s-spline')

Fig. 3.2

Matlab utilizează funcția spline pentru a găsi curba spline cubică asociată unei funcții f pentru care se cunosc valorile pe secvența de puncte amintită mai sus.

Să exemplificăm cu cazul funcției

» x = 0:2*pi;

» y = 4*cos(x)+1;

» xx = 0:.001:2*pi;

» yy = spline(x,y,xx);

» plot(x,y,'o',xx,yy)

» legend ('nodurile','spline cubic') ;

Fig. 3.3

Vom ilustra mai jos pentru acelasi domeniu de valori, interpolarea liniarǎ, pentru care vom utiliza funcția interp1, și interpolarea cubicǎ.

x=0:5;

» y=[0,20,60,68,72,100];

» xi=0:0.1:5;

» ylin=interp1(x,y,xi);

» ycub=spline(x,y,xi);

» plot(xi,ylin,':',xi,ycub,x,y,'o')

» legend('liniar','cubic');

» title('interpolare liniara si cubica');

Fig. 3.4

Sǎ urmǎrim în continuare cum calculǎm coeficienții unui spline cubic natural când se dau nodurile și valorile pe noduri . Vom aminti urmǎtoarea formulǎ de calcul pentru spline-ul cubic natural:

unde .

Expresia de mai sus dǎ expresia spline-lui cubic natural pe intervalul , în termenii primei și a doua derivate astfel :. Pentru calculul coeficienților , vom rezolva ecuația:

Pentru spline-ul cubic natural se pot impune și condițiile la granițǎ:

Având astfel coeficienții calculați se poate calcula spline-ul cubic natural .

Capitolul 4 Curbe de aproximare și interpolare

Metodele de aproximare cu polinoame determinate prin utilizarea celor mai mici pǎtrate sau aproximarea funcțiilor periodice folosind seriile Fourier pot fi considerate “clasice”, ele nelipsind din nici o lucrare care prezintǎ diferite metode numerice.

Odata cu folosirea calculatorului în proiectarea formelor a apǎrut necesitatea de a gasi noi metode atât din punct de vedere al aspectului matematic cât și al implementari lor în cadrul unor soft-uri ușor de folosit de neinițiați.

Astfel în ultimii 20 de ani a apǎrut o nouǎ disciplinǎ care îmbinǎ elemente de geometrie analiticǎ, diferențialǎ, desen tehnic, analizǎ matematicǎ și evident programarea – este vorba despre “grafica pe calculator” a cǎrei denumire în engleza este binecunoscuta “Computer Aided Geometric Design”. Curbele și suprafețele folosite în disciplina CAGD au inceput sǎ fie studiate în anul 1959 de matematicieni care lucrau pentru firme constructoare de automobile sau de avioane: P. de Casteljau (Citroen), P. Bezier (Renault), S. Coons (Ford), W. Gordon (General Motors), J. Ferguson (Boeing). Lucrǎrile lor de pionierat au fost urmate de numeroase alte studii si aplicații. Aceste realizǎri în grafica asistatǎ de calculator s-au bazat pe o serie de cercetǎri efectuate anterior de matematicieni care au studiat problema aproximǎrii funcțiilor, dintre care pot fi amintiți în special A.M. Legendre, C. Hermite, K. Wierstrass, P.L. Cebasev, C. Runge, H. Lebesque, D. Jackson, S.N. Bernstein s.a.

In present, proiectarea formei autovehiculelor, avioanelor, vapoarelor, a unei game foarte largi de produse este de neconceput fǎrǎ utilizarea unor soft-uri specifice de proiectare asistatǎ pe care le folosește un inginer. Este evident faptul cǎ pentru a reuși o realizare eficientǎ a acestor produse soft este necesar sǎ se cunoscǎ în cât mai mare mǎsurǎ și fundamentarea matematica a acestora.

Câteva din lucrǎrile de referințǎ în acest domeniu sunt amintite la bibliografie.

Ecuațiile analitice ale curbelor care descriu formele diverselor obiecte reale sunt de cele mai multe ori deosebit de complicate. In aceasta situație este însǎ suficient sǎ cunoaștem un set de puncte ale curbei respective și sǎ încercǎm sǎ conectǎm aceste puncte printr-o curbǎ a cǎrei ecuație analiticǎ este în general simplǎ. Aceste puncte le vom numi în continuare puncte de control.
In grafica pe calculator, ceea ce prezintǎ interes în primul rând este aspectul vizual al curbei reprezentate și mai puțin proprietǎțile sale matematice. Pornind de la un set de puncte de control, construcția unei curbe controlate de acestea poate fi realizatǎ utilizând metode de interpolare sau de aproximare. In primul caz, curba obținutǎ va trece prin punctele de control. In cel de-al doilea caz, curba nu va mai trece în mod necesar prin punctele de control ci prin vecinatatea acestora:

4.1 Curbe Bézier

Conceptul de curbǎ Bezier a fost dezvoltat pentru prima oara de Pierre Bezier prin anii ’60 când lucra ca matematician la firma de automobile Renault și au fost printre primele utilizate in proiectarea asistatǎ de calculator.

In planul axelor Oxy se considerǎ un numǎr de n+1 puncte (noduri sau puncte de control)

ale cǎror coordonate carteziene sunt cunoscute : prin unirea cǎrora se obține așa numitul poligon de control.

Dacǎ se noteazǎ

vectorii coloana ale cǎror elemente sunt coordonatele punctelor de control, coordonatele curbei Bezier se calculeazǎ astfel

(4.1.1)

Funcțiile se numesc polinoamele fundamentale Bernstein :

și sunt definite in intervalul [0,1]. Bernstein a utilizat pentru prima datǎ aceste polinoame pentru a demonstra teorema lui Weierstrass enunțatǎ în paragraful 1.2. Câteva din proprietǎțile importante ale acestor polinoame le vom defini mai jos:

1. polinoamele Bernstein sunt positive pentru orice

2. polinoamele Bernstein formeazǎ o partiție a unitǎții . Aceastǎ proprietate este foarte importantǎ în utilizarea acestor funcții în modelarea geometrica și aplicații grafice pe calculator.

3. fiecare polinom are un singur maxim pe intervalul în punctul .

4. orice polinom poate fi scris ca o combinație liniarǎ de polinoame de grad mai mare

Pentru calculul valorilor polinoamelor lui Bernstein se poate folosi urmǎtoarea formulǎ de recurențǎ:

Formula recursivǎ de mai sus stǎ la baza algoritmului de Casteljau binecunoscut, de evaluare a unei cure Bezier într-un punct

Algoritmul de Casteljau

Input: n+1 puncte de control

parametru real

Output: un punct de pe curba

pentru executa

// memorez punctele de control

pentru executa

pentru executa

Observație Pentru faptul cǎ polinoamele Bernstein sunt simetrice în raport cu si , forma unei curbe Bézier nu se schimbǎ indiferent de direcția de parametrizare.

In conditiile în care o curba Bézier nu posedǎ suficientǎ flexibilitate în modelarea formei dorite, putem îmbunǎtǎți aceastǎ caracteristicǎ prin creșterea flexibilitǎții poligonului de control. O soluție pentru această problema ar fi adǎugarea unui nou punct de control și obținem astfel o creștere a gradului curbei. Astfel pentru algoritmi care genereazǎ suprafețe pornind de la aceste curbe, e necesar ca aceste curbe sǎ aibǎ același grad. Folosind aceastǎ tehnica putem crește gradele curbelor inițiale la gradul maxim al familiei de curbe.

O altǎ aplicație importantǎ a tehnicii este în situația transferului de date între diferite sisteme grafice. Astfel dacǎ s-a generat o parabolǎ (curba Bezier de gradul doi) și vrem s-o transferǎm într-un sistem care recunoaște doar forme cubice va trebui sǎ creștem gradul parabolei. Am amintit în paragraful 1.2 un rezultat important din teoria aproximǎrii, teorema lui Weierstrass de aproximare, care poate fi acum reformulat în contextul curbelor parametrice ([Farin]). Fie o curbǎ continuǎ pe un interval [0,1] și secvența de valori . Aceste valori pot fi intepretate ca formând poligonul de control asociat curbei Bezier de grad :

ce aproximeazǎ curba

Daca vom creste gradul acestor curbe, obtinem secvența . Teorema lui Wierstrass în acest contex afirma ca

.

Cu alte cuvinte, teoretic, se exprima ideea ca orice curbǎ poate fi aproximatǎ arbitrar cu o curbǎ polinomialǎ. Polinoamele lui Bernstein au o serie de proprietǎți care vor fi utilizate în continuare, în cadrul aproximǎrii funcțiilor cu curbe Bezier. Astfel sunt valabile urmǎtoarele relații:

1.

2.

3.

Deoarece polinoamele Bernstein satisfac prima relație rezultǎ pentru curba Bézier urmǎtoarele relații

, ,

deci curba trece prin primul și ultimul punct , , al poligonului de control, dar nu este necesar sǎ treacǎ și prin celelalte puncte.

În cazul parcurgerii complete a curbei Bézier, parametrul va lua valori între 0 și 1 și daca pe curbǎ se vor gǎsi cele n+1 puncte de control rezultǎ cǎ aceste puncte se vor obține pentru valorile ale parametrului t.

Curba Bezier nu trece în mod necesar prin alt punct de control dar are proprietatea de a se menține tot timpul în interiorul poligonului convex definit de aceste puncte.

Deoarece punctelor de control le sunt asociate valori ale parametrului care diferǎ prin rezultǎ cǎ tangenta la curba Bezier în punctul este chiar prima latura a poligonului de control, același lucru fiind valabil și în cazul ultimului punct ce corespunde la , ultima latura a poligonului fiind tangentǎ la curba Bezier.

Utilizând o secvența de curbe Bezier de grad 3 se obține o curbǎ de aproximare formatǎ din porțiuni având proprietatea de a trece prin punctele si ale setului de puncte din control, panta tangentei la curbǎ în punctul, ceea ce permite construcția unei curbe ce trece printr-un set de puncte având forma controlatǎ de punctele intermediare:

Prin alegerea punctelor coliniare se va obține aspectul neted al curbei în.

Utilizând forma parametricǎ cu pentru fiecare curbǎ componenta determinatǎ de punctele vom avea:

Atunci când parametrul parcurge intervalul formula de mai sus furnizeazǎ punctele curbei dintre punctul corespunzǎtor punctului de control Pi și cel corespunzǎtor lui. Curba finalǎ se obține prin trasarea segmentelor de pe curbǎ pentru valorile i=1, i=4,…, i=3*(j-1)+1.

Vom da în continuare un cod de program Matlab, de desenare a unei curbe Bezier, Bezier.m. Funcția de desenare, va apela o altǎ funcție de calcul a polinoamelor Bernstein, bernstein.m .

function [X] = Bezier(C,M)

if nargin<2

M = 30;

end

n = length(C)-+1; %numarul punctelor de control

X=[];

%reprezentarea punctelor de control

for i = 0:M-1

u = i/(M-1); %valorile parametrului u

P = [0;0];

%evaluarea curbei Bezier

for k = 0:n

P = P + bernstein(n,k,u) .* C(:,k+1);

end

%adaugarea lui P la lista punctelor de pe curba

X = [X P];

end

%desenarea curbei

plot(X(1,:), X(2,:));

hold on

plot(C(1,:), C(2,:), 'ro–');

și funcția bernstein.m

function [bnk]=bernstein (n,k,x)

if k == 0

C = 1;

else

C = prod ([k+1:n])/(prod([1:n-k]));

end

bnk = C .* x.^k .* (1-x).^(n-k);

return

Sǎ exemplificǎm pentru urmǎtoarea matrice a punctelor de control:

» C=[0 1 2;1 0 1];

» Bezier (C,20)

Fig. 4.1

» C=[0 1 2 3;1 2 0 1];

» Bezier(C,20)

Fig. 4.2

» C=[0 2 3 4 5 6;1 5 5 0 6 1];

» Bezier(C,20)

Fig. 4.3

O altǎ metodǎ de evaluare a punctelor de pe o curbǎ Bézier este cea recursivǎ cea bazatǎ pe algoritmul De Casteljau amintit mai sus , dupǎ numele matematicianului Paul de Casteljau, care l-a dezvoltat la Citroën la sfarșitul lui 1950.

function drawbezier(pc,npc)

%deseneaza o curba Bezier

%drawbezier(pc,npc)

%pc – punctele de control

%npc – numar de puncte (rezolutia)

h=1/npc;

t=0:h:1;

uu=decast(pc,t);

plot(pc(1,:),pc(2,:),'–',pc(1,:),pc(2,:),'o',uu(1,:),uu(2,:));

Funcția de desenare a unei curbe Bézier va apela funcția în care se evalueazǎ curba Bezier recursiv, cu ajutorul algoritmului de Casteljau.

function q=decast(pc,t)

%algoritmul de Casteljau

%pc – punctele de control

%apel q=decast(pc,t)

%t punctele in care se evalueazǎ

[mp,np]=size(pc);

x=zeros(np,np);

y=zeros(np,np);

x(1,:)=pc(1,:);

y(1,:)=pc(2,:);

lt=length(t);

for l=1:lt

for i=2:np

for j=1:np-i+1

x(i,j)=t(l)*x(i-1,j)+(1-t(l))*x(i-1,j+1);

y(i,j)=t(l)*y(i-1,j)+(1-t(l))*y(i-1,j+1);

end

end

q(1,l)=x(np,1);

q(2,l)=y(np,1);

end;

4.2 Curbe B-spline

Polinoamele lui Bernstein, folosite în “construirea” funcțiilor (curbelor) Bezier sunt definite în intervalul [0, 1], care corespunde întregului interval în care se face aproximarea, corespunzând primului punct de control iar ultimului. S-a observat cǎ în cazul în care un punct de control (un vârf al poligonului de control) își schimbǎ poziția acest lucru are influențǎ asupra tuturor coordonatelor puntelor curbei Bézier. Acesta este și motivul pentru care în cazul în care sunt date un numǎr de puncte prin care (sau pe lângǎ care) se dorește sǎ treacǎ curba de aproximare trebuie multe încercǎri privind alegerea numǎrului de puncte de control precum și a pozițiilor acestora pentru ca sǎ se rezolve problema.

De asemenea un dezavantaj al funcțiilor lui Bézier îl constituie faptul cǎ gradul polinoamelor lui Bernstein crește cu creșterea numǎrului de puncte de control.

Odata numǎrul de puncte de control fixat ar fi de dorit ca prin modificarea poziției unui punct de control sǎ nu se schimbe alura curbei de aproximare în întregime ci numai local, în zona care corespunde punctului de control deplasat. Acest lucru este realizabil doar dacǎ se va adopta o altǎ bazǎ de reprezentare decât polinoamele lui Bernstein (definite pe tot intervalul) și anume niște polinoame numite polinoamele “B-spline”, care sunt nenule și positive numai într-un interval mǎrginit de valori ce corespund la douǎ , trei sau mai multe noduri, în rest fiind nule.

Pentru prima oarǎ matematicianul Schoenberg (nǎscut în Romania) a introdus și studiat aceste polinoame și funcțiile care rezultǎ prin combinarea lor, în anul 1946.

Funcțiile B-spline sunt folosite în generarea unor intrumente importante în grafica asistatǎ de calculator : curbele și suprafețe B-spline. Clasa curbelor B-spline este cel mai mult folositǎ în modelarea geometricǎ. Dezvoltarea inițialǎ a teoriei curbelor B-spline a folosit diferențe divizate, fiind foarte laborioasǎ din punct de vedere matematic. Ulterior s-a introdus o definiție geometricǎ, datoratǎ lui Carl de Boor si M. Cox, transformatǎ în una analiticǎ, în care se dǎ parametrizarea curbei intr-o bazǎ definitǎ recursiv.

O curbǎ B-spline este formata practic din mai multe arce de curbe Bezier ce se comportǎ ca un întreg. Punctele de legaturǎ dintre aceste curbe Bézier se vor numi noduri și vor juca un rol important în problemele legate de aspectul unei curbe B-spline, netezime sau flexibilitate.

Ca și în cazul curbelor Bézier și în cazul curbelor B–spline acestea se vor exprima parametric, utilizând parametrul

Vom considera un vector al nodurilor, (“knot vector”, în engleza) având componentele

valorile acestora fiind crescatoare Dacǎ aceste noduri sunt alese în mod uniform, sau în particular 1 atunci polinoamele B-spline și funcțiile care se creeazǎ cu acestea se numesc uniforme iar în caz contrar neuniforme.

Vom nota cu o funcție B-spline de gradul k, (“blending B-spline functions”, în engleza) și calculul lor se va face pe baza urmǎtoarei relații de recurențǎ:

și pentru

.

Remarca Ordinul k este independent de numarul punctelor de control . Astfel în cazul curbelor B-spline putem obține o mai mare flexibilitate a formei prin creșterea numǎrului punctelor de control fǎrǎ creșterea gradului curbei.

Notând cu punctele de control din plan, o curbǎ B-spline de gradul ce aproximeazǎ poligonul de control se definește astfel:

.

Algoritmul de evaluare a curbelor B-spline, algoritmul Cox de Boor este o extensie a algoritmului lui de Casteljau si are câteva avantaje fațǎ de acesta din urmǎ.Utlizând datele de intrare direct și calculând doar combinațiile convexe este mult mai stabil numeric. Este apoi mult mai flexibil pentru cǎ nu are în vedere o poziționare speciala a punctelor de control. Algoritmul se prezintǎ astfel sub formǎ triunghiularǎ, fiecare element se obține prin combinația convexǎ a douǎ elemente anterioare:

Algoritm Cox de Boor

Input: n puncte de control

– punctul in care se face evaluarea

– secventa nodurilor

Output valoarea

Fie

Pentru executa

sfarsit pentru

pentru executa

pentru executa

sfarsit pentru

.

Pentru o anumitǎ valoare a parametrului , funcțiile de bazǎ B-spline de grad k și un vector al nodurilor desemnat prin knot se pot calcula în Matlab cu urmǎtoarea funcție:

function alfa=Spline(i,k,knots,u)

if (k==1)if (knots[i]<=u)&& (u<knots[i+1]))

alfa=1;

else

alfa=0;

end;

else if (knots[i+k-1]==knots[i])

alfa=(knots[i+k]-u)/(knots[i+k]-knots[i+1])*Spline(i+1,k-1,knots,v);

else

if (knots[i+k]==knots[i+1])

alfa= (v-knots[i])/(knots[i+k-1]-knots[i])*Spline(i,k-1,knots,v);

else alfa=(v-knots[i])/(knots[i+k-1]-knots[i])*Spline(i,k-1,knots,v)+(knots[i+k]-v)/(knots[i+k]-knots[i+1])*Spline(i+1,k-1,knots,v);

end;

end;

end;

end;

Sǎ reprezentǎm în continuare poligonul punctelor de control

» pc=[5 10 20 30 35 25 20;20 10 5 15 25 50 45];

și curba B-spline de grad 3, asociatǎ acestor puncte de control. Funcția de desenare a poligonului și a curbei B-spline în același sistem de referințǎ este datǎ de :

function draw_b_spli(pc,g,npd,stil)

%draw_b_spline(pc,grad,npd,stil)

%pc – punctele de control

%g – gradul

%npd – numar de puncte(rezolutia)

%stil – stilul de linie(conform sintaxei plot implicit '-')

n=length(pc);

kn=[0 0 0 0 1 2 3 4 4 4 4];

%kn=noduri(n-1,g);

ln=length(kn);

if nargin<2

error('numar ilegal de argumente');

else

if nargin==2

npd=100;

stil='-';

else

if nargin==3

stil='r-';

end

end

end

pas=kn(ln)/npd;

t=0:pas:kn(ln);

q=coxboor(g,pc,kn,t);

%plot(q(1,:),q(2,:),stil);

plot(pc(1,:),pc(2,:),'–',pc(1,:),pc(2,:),'o',q(1,:),q(2,:),stil,'LineWidth',2.5);

hold on

plot(pc(1,:),pc(2,:),'–',qnew(1,:),qnew(2,:),'r–')

Se observǎ cǎ funcția apeleazǎ funcția coxboor(g,pc,kn,t), care calculeazǎ valoarea B-spline-lui pentru o valoare datǎ a parametrului t. Calculul se face recursiv prin algoritmul lui CoxDeBoor, astfel:

function q=coxboor(grad,pc,knots,t)

%algoritmul Cox-DeBoor pentru B-spline

%q=coxboor(grad,pc,knots,t)

%grad- gradul spline-ului

%pc – poligonul de control

%knots – nodurile

%t – abscisele punctelor in care calculam spline-ul

m1=length(t);

[lin,n]=size(pc);

q=zeros(lin,m1);

for l=1:m1

%determina pozitia lui t in raport cu nodurile

for j=1:n+grad

if (t(l)>=knots(j)) & (t(l) < knots(j+1))

break

end %if

end %for j

if (j==n+grad)

j=n;

end;

qv=pc

qn=zeros(lin,n);

for r=0:grad-1

for i=j-grad+r+1:j

%sprintf('omega(%g,%g),%g',i,grad-r,r)

w=omega(i,grad-r,knots,t(l));

qn(:,i)=(1-w)*qv(:,i-1)+w*qv(:,i)

end %for i

qv=qn;

end %for r

%sprintf('t(%g)=%g',l,t(l))

q(:,l)=qv(:,j);

end %for l

return

unde

function u=omega(i,k,knots,t)

% sprintf('i=%g, k=%g\n',i,k)

% pause

if knots(i+k) ~= knots(i)

u=(t-knots(i))./(knots(i+k)-knots(i));

else

u=0;

end %if

Apelând în fereastra de comandǎ a Matlab-lui funcția

» draw_b_spli (pc,3,100,'–k')

cu punctele de control specificate mai sus, obținem urmǎtorul grafic:

Fig. 4.5

Dupǎ cum s-a vǎzut o curbǎ B-spline este determinatǎ de un poligon de control și un șir de noduri. Din aceastǎ perspectivǎ importantǎ devine problema influenței pe care o are modificarea poziției unui punct de control sau a multiplicitǎții unui nod. Aceste curbe fiind folosite în design-ul liber, deseori dupǎ alegerea poligonului de control și a șirului de noduri, forma curbei generate nu corespunde fie din punct de vedere estetic fie funcțional sau dimensional și o primǎ încercare ar fi “mutarea” unui punct de control. Aceastǎ operație va avea doar efect local , așa cum se observǎ și din formula de calcul a unui punct de pe curba B-spline, conform algoritmului Cox De Boor și acesta reprezintǎ un mare avantaj al curbelor B-spline. Dacǎ modificarea poziției unui punct de control are efecte previzibile asupra formei curbei, modificarea poziției unui nod în șirul de noduri este și imprevizibilǎ și nesatisfǎcǎtoare. Schimbând însǎ ordinul de multiplicitate al unui nod putem obține niște obiective dorite asupra formei curbei. Cel mai important aspect este acela cǎ dacǎ un nod are ordinul de multiplicitate k adicǎ atunci curba va conține punctul . Aceastǎ proprietate poate fi folositǎ in problemele de “free design” în care se impune ca viitoarea curbǎ sǎ treacǎ prin niste puncte fixate.

Bibliografie

C de Boor , A practical Guide to Splines, New-York, Springer, 1978;

W. Boehm,   G. Farin,   J. Kahmann,  A survey of curve and surface methods in CAGD, Computer Aided Geometric Design , 1  (1984)  pp. 1–60;

W. Boehm, Inserting new knots into B-spline curve, Computer Aided Design, 12(4), 199-201, 1980;

R.L. Burden, J.D. Faires, Numerical Analysis, Third Edition , Prindle, Weber&Schimdt, Boston, 1985;

I. Chiorean, T. Cǎtinaș, Gh. Coman, Advanced Course on Numerical Analysis, Ed. Presa Univ. Clujeană, 2007

Gh. Coman, Analizǎ Numericǎ, Ed. Libris, Cluj-Napoca, 1995;

S. D. Conte, Carl de Boor, Elementary Numerical Analysis, An Algorithmic Aproach, International Series in Pure and Applied Mathematics, McGraw- Hill Book Company, New-York, third edition, 1980

I. Cuculescu, Analizǎ numericǎ, Ed. Tehnicǎ, București, 1967;

G. Farin, Curves and Surfaces for Computer Aided Design. A practical Guide, Academic Press, London, 1988;

D. Foley, A. van Dam, S.K.Feiner, J.F. Hughes, and Phillips: Introduction to Computer Graphics, Addison-Wesley, 1994;

D. Foley , A. van Dam, S.K. Feiner, J.F. Hughes, Computer Graphics: Principles and Practice in C (2nd ed., Addison Wesley, 1992;

M. Ghinea, V. Fireteanu, Matlab- Calcul Numeric, grafica, aplicații, Ed. Teora,1997;

J. Gravesen, Differential Geometry and Design of Shape and Motion, Department of Mathematics, Techical University of Denmark;

G.G. Lorentz, Bernstein polynomials, Univ. Toronto Press  (1953);

G.G Lorentz, Aproximation of functions, Holt, Rinehart and Winston, New York, 1966;

J. Lowther, C. Kuang Shene, Teaching B-spline is not difficult!, Michigan Technological University, 2003;

C.V. Muraru, Metode Numerice. Seminarii Matlab, Ed. Edusoft, Bacau, 2005

Theo Pavlidis, Algorithms for Graphics and Image Processing, Spronger-Verlang, Berlin-Heidelberg, Computer Science, Inc., 1982;

T. Sederberg, BYU Bézier curves, http://www.tsplines.com/resources/class_notes/Bezier_curves.pdf;

D.D Stancu , Ghe. Coman, O. Agratini, R. Trîmbițaș, Analiza Numericǎ și Teoria Aproximǎrii, Vol. I,.Presa Universitarǎ Clujeanǎ, Cluj-Napoca, 2001;

M Talmaciu, Metode Numerice, Vol. I, Ed. Tehnicǎ Info, Chișinǎu, 2005;

http://cgg-journal.com/1999-1/01/Denis.htm;

http://math.fullerton.edu/mathews/n2003/BezierCurveMod.html;

R. Trâmbițaș, Analizǎ Numericǎ, O introducere bazatǎ pe Matlab, Ed. Presa Universitarǎ Clujeanǎ, 2005;

R. Trâmbitaș, Numerical Analysis, Ed.Presa Univ.Clujeanǎ, 2007

T. Vladislav, I. Rașa, Analizǎ Numericǎ, Ed. Tehnicǎ, București, 1997;

A. Watt, 3D Computer Graphics (3rd edition), Addison-Wesley, 2000;