Curs 6+7

Interpolarea şi aproximarea funcţiilor

cu aplicaţii în ingineria electrică

METODE NUMERICE

În aplicaţiile din domeniul electrotehnic nu se cunoaşte expresia analitică a funcţiei care trebuie aproximată ci doar valorile ei într-un anumit număr de puncte (tabelate - obţinute din calcule sau măsurători experimentale) urmărindu-se determinarea aproximativă a valorilor corespunzătoare unor alte puncte diferite de cele date.

Exemple: curbele experimentale caracteristice care descriu funcţionarea aparatelor electrice, curbele de sarcină, lucru mecanic, inducţia în funcţie de câmp, distribuţii de sarcină etc.

Aproximarea unei funcţii exprimată analitic sub forma unor formule explicite, implicite sau parametrice, sub forma unor serii, sau a unui algoritm se face cu scopul simplificării calculelor de evaluare a mărimii funcţiei , a derivatelor acesteia sau a integralei definite.

Evaluarea unei funcţii definită sub formă numerică (dată tabelar) în urma unor măsurători experimentale, presupune aproximarea ei (interpolarea) în intervalele dintre nodurile reţelei în orice punct al domeniului de definiţie.

Cea mai simplă metodă de interpolare a unei funcţii definită sub formă numerică prin coordonatele (xi,yi) ale unor puncte numite noduri, constă în aproximarea funcţiei cu un polinom pentru a putea fi prelucrată în continuare (interpolare, derivare, integrare etc) evaluarea funcţiei reducându-se la operaţii aritmetice elementare (adunări şi înmulţiri).

Se măsoară la momente discrete x0, x1...,xn (noduri), valorile unor funcţii f(x) şi se pune problema de a găsi valorile sale în alte puncte diferite de noduri.

)(...)()(),...,,( 11000 xgaxgaxgaaaxg nnn +++=

nnn xaxaxaaaaaxg ++++= ...),...,,;( 2

21010- interpolare polinomială

1 ,...),...,,;( 221010 −=++++= ieaeaeaaaaaxg nxi

nxixi

n

- interpolare trigonometrică

mm

nn

mn xbxbbxaxaa

bbaaxg++++++

=......

),...,,,...,;(10

1000 - interpolare raţională

- interpolare liniară

Funcţia de aproximare este de forma:

),...,,;()( 10 naaaxgxf ≅ - model matematic.

Dacă nu există informaţii asupra problemei tehnice care a generat

modelul matematic, atunci cel mai des se utilizează pentru interpolare

polinoame!

Avantaje:

valoarea polinoamelor se calculează uşor;

sumele, diferenţele, produsele de polinoame au ca rezultat polinoame;

prin derivare şi integrare (care se fac uşor), rezultă tot polinoame;

teoria interpolării polinomiale este simplă şi bine pusă la punct.

Aproximarea prin polinoame a funcţiilor

[ ] [ ],Cf,b,a b,a∈ ( )xp,0 n∃>ε∀

[ ]b,ax∈∀ ( ) ( ) ε<− xpxf n

polinom de gradul nFie f o funcţie definită şi continuă pe intervalul

Se porneşte de la teorema lui Weierstrass care are următorul enunţ:

În aplicaţiile electrotehnice alegerea funcţiei de aproximare se bazează şi pe cunoaşterea formei funcţiei care trebuie aproximată ţinând cont de informaţiile asupra aplicaţiei practice care a generat modelul matematic.

Problema care se pune este determinarea polinomului pn(x) care satisface relaţia de mai sus.

Interpolarea polinomială

[ ] ( ) ( ) ( ) ji,xx,xf...,,xf,xfb,ax...,,x,x jin10n10 ≠≠→∈

( ) ( ) n,0i,xfxp iin ==

Cunoscând valorile funcţiei f determinate experimental prin măsurători în punctele

se pune problema determinării valorilor funcţiei în alte puncte adică să se găsească un polinom pn astfel încât:

- polinomul pn(x) trebuie să coincidă cu funcţia f(x) pe n+1 puncte.

Unicitatea polinomului de interpolare

Se va demonstra că există un polinom şi numai unul de grad mai mic ca n care îndeplineşte această condiţie.

( ) nn

2210n xa...xaxaaxp ++++=

( ) ( ) n,0i,xfxp iin ==

Demonstraţie: Dacă avem polinomul de grad n:

Condiţia devine:

( )( )

( )⎪⎪

⎩

⎪⎪

⎨

⎧

=++++

=++++

=++++

nnnn

2n2n10

1n1n

212110

0n0n

202010

xfxa...xaxaa.....................................................xfxa...xaxaa

xfxa...xaxaa

n10 a,...,a,a

( ) 0xx

x...x1............x...x1x...x1

dn

jxixji

0j,iji

nnn

n11

n00

s ≠−== ∏≠

≠=

- necunoscute

Determinantul sistemului este Vandermonde:

Regula lui Cramer de rezolvare a sistemelor de ecuaţii, sistemul poate fi compatibil şi unic determinat adică să aibă soluţie unică: coeficienţii unic determinaţi:

( )( )

( )

( )( )

( )

( )( )

( )s

nn

11

00

ns

nnn

n11

n00

1s

nnnn

n111

n000

0 dxf...x1............xf...x1xf...x1

a...,d

x...xf1............x...xf1x...xf1

a,d

x...xxf............x...xxfx...xxf

a ===

( ) ( ) n,0i,xfxp iin ==

Există un singur polinom care interpolează funcţia în nodurile x0, x1,...,xn pentru care:

n10 a...,,a,a

Problema interpolării presupune parcurgerea etapelor:

determinarea coeficienţilor polinomului de interpolare prin rezolvarea unui sistem liniar de ecuaţii algebriceevaluarea polinomului de interpolat

Această variantă de interpolare poate fi aplicată doar pentru valori miciale gradului polinomului (n<5) deoarece are două mari dezavantaje:

Efort de calcul mare pentru determinarea coeficienţilor (Cramer).Erorile soluţiei sunt mari deoarece sistemul poate fi rău condiţionat pentru valori mari a gradului polinomului

Metoda de interpolare bazată pe polinomul de interpolare Lagrange elimină dezavantajele metodei clasice de interpolare polinomială, în schimb timpul necesar evaluării polinomului de interpolare creşte de la ordinul liniar O(n) la cel pătratic O(n2).

n][0,i )( , ∈= iii xfyx

Fie o funcţie f(x) definită pe intervalul [a,b], ale cărei valori yi sunt cunoscute numai în nodurile

- interpolare liniară

ix - noduri de interpolare

Interpolarea Lagrange

Tabel de valori - din măsurători experimentale

ii y)x(p =

]b,a[x )x(p)x(f ∈∀≅

]n,0[i ,y)x(p)x(f iii ∈==

Se numeşte polinom de interpolare asociat tabelului de valori un polinom p de grad mai mic sau egal cu n, cu coeficienţi reali astfel încât:

are loc formula aproximativă:

unde p(x) este unic pentru un tabel dat iar f şi p au aceleaşi valori în nodurile fixate.

ixx ≠Observaţie: Polinomul p(x) permite calculul valorilor sale şi în punctele

deci între noduri, ceea ce justifică denumirea de interpolare.

pk0 ≤≤ ])n,0[i( li ∈∀

( ) ik,0xl ki ≠=( ) 1xl ii =

⎩⎨⎧

=≠

=δ=ik ,1ik ,0

)x(l kiki

Pentru orice k, polinomul ce verifică condiţiile

condiţii ce pot fi scrise cu ajutorul simbolului lui Kroneker astfel:

,se notează cu

n1i1i10 x,...,x,x,...,x,x +− ( ) ( ) ( )( ) ( )n1i1i0ii xx...xxxx...xxcxl −−−−= +−

ixx = ( ) 1xl ii =

( ) ( )( ) ( )ni1ii1ii0ii xx...xxxx...xx

1c−−−−

=+−

Deoarece se cere ca li să se anuleze în n puncte rezultă că acesta are forma:

- ci este o constantă

rezultă constanta:Punând în această relaţie şi ţinând cont că

Înlocuind constanta ci în relaţia de mai sus rezultă – polinomul bază:

( ) ( ) ( )( ) ( )( ) ( )( ) ( ) ∏

≠=+−

+−

−−

=−−−−

−−−−=

n

ij0j ji

j

ni1ii1ii0i

n1i1i0i xx

xxxx...xxxx...xx

xx...xxxx...xxxl

( )}

)x(L)x(ly)x(lxf)x(P n

notn

0iii

n

0iiin =⋅=⋅= ∑∑

==

)(xLn

∑= +−

+−

−−−−−−−−

⋅=n

0i ni1ii1ii0i

n1i1i0in )xx)...(xx)(xx)...(xx(

)xx)...(xx)(xx)...(xx(y)x(L

( ) ( ) i

n

0ikii

n

0ikiiin yyxlyxL =δ⋅=⋅= ∑∑

==( ) ( )iiin xfyxL ==

Se construieşte polinomul de interpolare Lagrange de grad cel mult n:

- se numeşte polinomul Lagrange de interpolare:

( ) ( )xLxf n≅ ( ) ( ) ni0,yxLxf iini ≤≤==

Deci Ln(x) este polinomul Lagrange de interpolare asociat tabelului de valori

Formula este sugerată de faptul că

adică funcţia f şi polinomul Ln au aceleaşi valori în nodurile fixate.

( ) ( )1010 xf,xfx,x →

( ) ( )01

01

10

10 xx

xxxl,xxxxxl

−−

=−−

=

1. Exemplificăm o interpolare polinomială de gradul I

Dacă n=1 atunci se cunoaşte funcţia în două noduri

Vom avea polinom de gradul 1:

( ) ( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( )001

010

01

01

10

1011001

xxxx

xfxfxf

xxxx

xfxxxxxfxlxfxlxfxL

−⋅−−

+=

=−−

⋅+−−

⋅=⋅+⋅=

( ) ( ) ( ) ( ) 0xl,1xl,0xl,1xl 01111000 ====

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 111011110011

001001100001

yxf1xf0xfxlxfxlxfxLyxf0xf1xfxlxfxlxfxL

==⋅+⋅=⋅+⋅===⋅+⋅=⋅+⋅=

Deci formula de interpolare liniară este:

Astfel L1(x) este polinomul unic care trece prin punctele (x0,y0) şi (x1,y1) şi aproximează funcţia f(x) pe intervalul [x0, x1].

Se trece la generalizarea conceptului de interpolare liniară prin determinarea unui polinom de grad cel mult n care trece prin n+1 puncte în care funcţia se cunoaşte:

( )( ) ( )( ) ( )( )nn1100 xf,x,...,xf,x,xf,x

Fig. 4 Interpolarea Lagrange de ordinul n

Fig. 3 Interpolarea Lagrange de ordinul I

( ) ( ) ( ) ( ) ( )001

0101 xx

xxxfxf

xfxL −⋅−−

+=

)x(R)x(L)x(f 11 += ( ) ( )xLxf 1≈

x

( ) ( )xLxf 1≈ ( ) ( ) hxfxL 01 +=

Considerăm formula de interpolare liniară Lagrange:

unde R1(x) este restul în formula de interpolare liniară Lagrange.

avem:

Fig. 5 Interpolarea liniară – demonstraţie

În punctul intermediar

( ) ( )( ) ( ) ( )0

01

01

01

0

01

xxxx

xfxfhxxxx

xfxfh

−⋅−−

=⇒−−

=−

( ) ( ) ( ) ( ) ( ) ( ) ( )001

01001 xx

xxxfxfxfhxfxfxL −⋅

−−

+=+==

Din asemănarea celor două triunghiuri:

)()()( xRxLxf nn +=

)xx)...(xx)(xx()!1n()(f)x(R n10

)1n(

n −−−+ξ

=+

- restul Lagrange de ordin n

Formula de interpolare a lui Lagrange de ordin n:

ξ - este un punct din cel mai mic interval care conţine nodurile x0, x1,...xn si variabila x.

.

Determinarea restului in dezvoltarea in serie Taylor

Demonstraţie pe tabla! .

Determinarea restului Lagrange

Demonstratie pe tabla!

xe)x(f =⎭⎬⎫

⎩⎨⎧ −− 1,

21,0,

21,1 - noduriExemplu:

Se interpolează pe intervalul [-1,1] funcţia dată prin polinomul Lagrange Ln(x):

( )( ) ( )

( ) ( )

( ) ( ) ( )

( ) ( )

( ) ( )

( ) ( ) ( )

( ) ( )

( ) ( )1

21

0

21

1n

e

21101

21111

21x0x

21x1x

e1

210

21

21

211

21

1x0x21x1x

e10

210

21010

1x21x

21x1x

e1

21

21

210

211

21

1x21x0x1x

e11

21101

211

1x21x0x

21x

xL

⋅⎟⎠⎞

⎜⎝⎛ −⋅+⋅⎟

⎠⎞

⎜⎝⎛ +⋅+

⎟⎠⎞

⎜⎝⎛ −⋅−⋅⎟

⎠⎞

⎜⎝⎛ +⋅+

+

+⋅⎟⎠⎞

⎜⎝⎛ −⋅⎟

⎠⎞

⎜⎝⎛ −⋅⎟

⎠⎞

⎜⎝⎛ +⋅⎟

⎠⎞

⎜⎝⎛ +

−⋅−⋅⎟⎠⎞

⎜⎝⎛ +⋅+

+⋅−⋅⎟

⎠⎞

⎜⎝⎛ −⋅⎟

⎠⎞

⎜⎝⎛ +⋅+

−⋅⎟⎠⎞

⎜⎝⎛ −⋅⎟

⎠⎞

⎜⎝⎛ +⋅+

+

+⋅⎟⎠⎞

⎜⎝⎛ −−⋅⎟

⎠⎞

⎜⎝⎛ −−⋅⎟

⎠⎞

⎜⎝⎛ −−⋅⎟

⎠⎞

⎜⎝⎛ +−

−⋅⎟⎠⎞

⎜⎝⎛ −⋅−⋅+

+⋅−−⋅⎟

⎠⎞

⎜⎝⎛ −−⋅−−⋅⎟

⎠⎞

⎜⎝⎛ +−

−⋅⎟⎠⎞

⎜⎝⎛ −⋅−⋅⎟

⎠⎞

⎜⎝⎛ +

=−−

( ) 432n x04.0x17.0x49.0x99.00.1xL ⋅+⋅+⋅+⋅+=

x1)x(f = { }4x,5.2x,0x 210 ===

( ) ( ) ( )( ) ( ) ( )

( ) ( ) ( )( ) ( )

( )

( ) ( ) ( )( ) ( )

( )3

5x5.4x5.24245.2x2xxl

332x24x4

45.225.24x2xxl

10x5.6x425.224x5.2xxl

2

1

0

+−=

−⋅−−⋅−

=

−+−=

−⋅−−⋅−

=

+−=−⋅−−⋅−

=

( ) ( ) ( ) ( ) ( ) 15.1x42.0x05.0xfxlxLxP2

0iii2 +⋅−=⋅== ∑

=

( ) ( ) 325.03L3f 2 =≅

Exemplu

Se interpolează pe intervalul [1,5] funcţia dată prin polinomul Lagrange de ordinul II:

O aproximaţie pentru

este

Interpolare Lagrange de ordinul II

( ) ( ) ( ) ( )( ) ( ) 25.04fxf

;4.05.2fxf;5.02fxf

2

10

======

( ) 333.0313f ==

Alte exemple:

1. Să se determine o aproximare polinomială pentru caracteristica unei diode semiconductoare:

( ) ( ) ( )V027.0v,A10I,1eIui 6s

vu

s ==⎟⎟⎠

⎞⎜⎜⎝

⎛−= −

şi analiza erorii introdusă prin interpolare.

2. Să se determine o aproximare polinomială pentru caracteristica de magnetizare:

( ) 04

s7

0s

s0 10,T7.1B,104,BHthBHB μ=μ=π=μ⎟⎟⎠

⎞⎜⎜⎝

⎛ μ+μ= −

şi analiza erorii introdusă prin interpolare.

Interpolarea Newton cu diferenţe divizatePresupunem că Ln(x) este polinomul Lagrange de grad n care interpolează funcţia f(x) înnodurile distincte x0, x1,…,xn

Diferenţele divizate ale funcţiei f(x) sunt utilizate pentru a exprima polinomul Ln(x) sub forma lui Newton Nn(x), polinomul se va scrie sub o formă generalizată sub forma lui Taylor:

( ) ( ) ( ) ( )( ) ( )( ) ( )1n10n102010nn xx...xxxxa...xxxxaxxaaxNxL −−−−++−−+−+==

Pentru determinarea constantelor a0, a1,...an, evaluăm polinomul în nodurile:

( ) ( )

( ) ( ) ( ) ( ) ( ) ( )01

01111n01101

00n00

xxxfxfaxfxNxxaxfxx

xfxNaxx

−−

=⇒==−+⇒=

==⇒=

Tabelul de valori rezultat din măsurători experimentale:

n,...,1i xx

)x(f)x(f

1ii

1ii =−−

−

−

[ ] ]f;x,x[sau x,xf i1ii1i −−

- diferenţe divizate de ordinul I ale lui f pe nodurile xi

- se notează

.

[ ] [ ] [ ]1iki

1kii1iki1iikii1i xx

x,...,x,xfx,...,x,xfx,...,x,xf−+

−+−+++− −

−=

[ ] ( )00 xfxf =

[ ] [ ] [ ]01

0110 xx

xfxfx,xf−−

=

[ ] [ ] [ ]02

1021210 xx

x,xfx,xfx,x,xf−−

=

[ ] [ ] [ ]0n

1n10n21n210 xx

x,...x,xfx,...x,xfx,...x,x,xf−−

= −

Diferenţele divizate de ordinul n pe n+1 noduri:

Diferenţele divizate de ordinul (k+1) se definesc:

Diferenţele divizate de ordinul 0 pe un nod

Diferenţele divizate liniare de ordinul 1 pe două noduri:

Diferenţele divizate de ordinul 2 pe trei noduri:

.

( ) ( ) [ ]1001

011 x,xf

xxxfxfa =

−−

= [ ]k10k x,...,x,xfa =Se notează astfel:

[ ] [ ]0110 x,xfx,xf = [ ] [ ]n10n10 x,...,x,xfx,...,x,xf ααα=

[ ] [ ] [ ]10

0

01

110 xx

xfxx

xfx,xf−

+−

=

[ ] [ ]( )( )

[ ]( )( )

[ ]( )( )2010

0

2101

1

1202

2210 xxxx

xfxxxx

xfxxxx

xfx,x,xf−−

+−−

+−−

=

[ ] [ ]( )( ) ( )∑

= −−−=

n

i niii

in xxxxxx

xfxxxf0 10

10 ...,...,,

Proprietăţi:

[ ] [ ]100 x,xfx,xf =

[ ] [ ] [ ] [ ] [ ] [ ] ( )0000

00 xxx,xfxfxf

xxxfxfx,xf −⋅+=⇒

−−

=

a. Pentru f funcţie liniară

- o funcţie liniară pentru care valoarea derivatei este constantă.

- funcţia şi polinomul coincid

[ ] [ ] [ ] ( ) ( )xRxxx,xfxfxf 10100 +−⋅+=

( ) ( ) [ ] [ ]( ) ( )( ) [ ] [ ]( )

[ ]444 3444 21

1x,0x,xf

1

1001010001 xx

x,xfx,xfxxxxx,xfx,xfxxxR−−

⋅−−=−⋅−=

b. Pentru f funcţie neliniară

( ) [ ] [ ] [ ] ( ) ( ) ( ) [ ] [ ]( ) [ ]⎟⎟

⎠

⎞⎜⎜⎝

⎛−

−−

⋅−=⇒−⋅−−= 100

00101001 x,xf

xxxfxf

xxxRxxx,xfxfxfxR

( ) ( )( ) [ ]10101 x,x,xfxxxxxR ⋅−−=Polinomul construit este polinomul de gradul 2 de interpolare sub forma Newton cu diferenţedivizate considerând nodul x2 nod intermediar:

( ) [ ] ( ) [ ]( ) ( ) ( ) ( )

( )( ) [ ]( )

4444 34444 214444 34444 21x1R

10210

xfx0Rxfx1N

10002 x,x,xfxxxxx,xfxxxfxN −−+−+==−=

( ) [ ] ( ) [ ]( ) ( ) ( ) ( )

( )( ) [ ]( )

( )

( ) ( ) [ ]n01n0

xN

xR

21010

xfxRxfxN

1000n x,...,xfxx...xx...x,x,xfxxxxx,xfxxxfxN

2

101

−

=−=

−−++−−+−+=

4444444444 34444444444 2144444 344444 214444 34444 21

( ) ( ) ( ) [ ]n10n0n x,...,x,x,xfxx...xxxR −−=

Observaţii:S-a arătat că în cazul în care se adaugă un nod de interpolare, faţă de cazul interpolării Lagrange, interpolarea Newton nu necesită întreaga recalculare.Metoda de interpolare Newton realizează un compromis între evaluare şi construcţia algoritmului care este stabil din punct de vedere al erorilor numerice.Permite mărirea gradului polinomului de interpolare prin adăugarea unui nod nou în reţeaua de interpolare cu reutilizarea coeficienţilor de la gradul anterior care nu se modifică.Coeficienţii polinomului Newton reprezintă diferenţele divizate ale funcţiei de interpolat ceea ce uşurează calculul numeric al derivatelor polinomului de interpolare.Timpul de calcul este dependent de eroarea impusă având valori mari doar în cazurile în care se doreşte o precizie ridicată.

Exemplu: Se cunosc valorile funcţiei în punctele

( ) ( ) ( ) ( ) 554f,253f,11f,50f ===−=

Să se calculeze valoarea funcţiei în 0.5 utilizând un polinom de interpolare Newton de gradul I.

( ) [ ] ( ) [ ]10001 x,xfxxxfxN ⋅−+= ( ) ( ) 2605.055.0N1 −=⋅−+−=

( ) ( )( ) [ ]10101 x,x,xfxxxxxR ⋅−−= ( ) ( )( ) 5.0215.005.05.0R1 −=⋅−−=

( ) ( )xsinxf π=21x,

61x,0x 210 ===

Exemplu: Să se scrie polinomul de interpolare Newton cu diferenţe divizate de gradul II pentru funcţia

şi nodurile

.

Aproximarea cu abatere medie pătratică minimă-metoda celor mai mici pătrate- least squares data fitting

Se utilizează pentru aproximarea funcţiilor definite prin noduri ale căror coordonate prezintă un oarecare grad de incertitudine. Este cazul funcţiilor ce exprimă dependenţe obţinute experimental prin măsurători sau ca urmare a unor calcule care folosesc rezultatele măsurătorilor, nu se recomandă întrebuinţarea metodelor de aproximare Lagrange care uneori pot amplifica aceste erori.

Prin folosirea unui polinom de aproximare al cărui grad m este mai mic decât cel al polinomului Lagrange corespunzător se obţine o aproximare cu abatere pătratică minimă.

Metoda este cunoscută şi sub numele metoda celor mai mici pătrate a lui Gauss.

Polinomul de aproximare P(x) nu trece prin toate cele n noduri sau chiar prin nici unul dar este cel mai aproape de toate acestea.

Metoda urmăreşte minimizarea erorii calculate de norma euclidiana (suma pătratelor abaterilor dintre datele experimentale şi cele determinate teoretic).

Determinarea aproximării în sensul celor mai mici pătrate se reduce la rezolvarea unui sistem de ecuaţii algebrice liniare cu un număr de ecuaţii mai mare decât numărul de necunoscute, care este supradeterminat. Acest sistem este de obicei incompatibil deci graficul polinomului de interpolare nu trece prin punctele (xk,yk) ci printre ele.

În general datele experimentale sunt afectate de erori de măsurare, astfel încât reconstrucţia unei caracteristici corespunzătoare unui fenomen fizic sau proces tehnic, în loc să treacă prin fiecare punct de măsură trebuie să urmărească tendinţa rezultatelor experimentale cu o eroare generală cât mai mică.

Exemple: Aproximaţi prin polinoame de grad: 3, 5, 10, 15, 20-curba de magnetizare B(H)-caracteristica unei diode semiconductoare i(u)-variaţia puterii disipate de un rezistor P(R) într-un circuit electric de curent continuu;

Metoda celor mai mici pătrate se aplică nu numai la aproximarea caracteristicilor neliniare dar şi la rezolvarea problemelor inverse, unde conceptul este echivalent cu pseudosoluţia.

Aplicaţia 1

n \\ numarul de receptoare;

n 3 4 5 6 9 10 12 15 17 20 22 24 25 27 30 33 35 40( ):=

ka \\ coeficientul de influenta;

ka 1 1.4 1.8 2.35 2.5 3 3.4 4 4.5 5.1 5.8 6 6.35 6.9 7.2 8 8.5 9.2( ):=

F x( )

log x( )

x2

x

1

⎛⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎠

:= S linfit nT kaT, F,⎛

⎝⎞⎠:= S

0.565

4.201− 10 4−×

0.222

0.216

⎛⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎠

=

f x( ) F x( ) S⋅:= x 3 3.01, 40..:=

0 10 20 30 40

5

10

f x( )

k a

x n,

Metoda II. Funcţia linfit

i 0 last nT( )..:= Abatere

i

f nT( )i

⎡⎣

⎤⎦ ka

T⎛⎝

⎞⎠i

−⎡⎢⎣⎤⎥⎦

2

∑:= Abatere 0.385= Abatere 0.5< 1=

(verificare)

dispersie max f nT( ) kaT−

→⎯⎯⎯⎯⎯⎛⎜⎝

⎞⎟⎠ 1.1⋅:= dispersie 0.415=

10 20 30 40

0

Dispersia punctelor fata de o referinta

0

Aproximarea analitica liniara prin minimizarea abaterii medii patratice

\\ datele tabelare provenite din masuratori experimentale pentru care trebuie realizata aproximarea analitica liniara, cu conditia de minimizare a abaterii medii patratice;

x

0

1

2

3

4

5

⎛⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎠

:= y

10

25

51

66

97

118

⎛⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎠

:= n last x( ):= n 5= k 0 n..:=

f x( ) m x⋅ b+ \\ functia analitica de identificat

F m b,( )

0

n

k

yk m xk⋅ b+( )−⎡⎣ ⎤⎦2∑

=

:= \\ formula abaterii medii patratice

\\ conditia de minimizare a abaterii medii patratice;m

F m b,( )dd

0b

F m b,( )dd

0

Sistemul de ecuatii care determina coeficientii functie de gradul I, de aproximare:

m 1:= b 1:= \\ initializarea solutiilorGiven

m

0

n

k

xk( )2∑=

⎡⎢⎢⎣

⎤⎥⎥⎦

⋅ b

0

n

k

xk∑=

⎛⎜⎜⎝

⎞⎟⎟⎠

⋅+

0

n

k

xk yk⋅∑= \\ sistemul de ecuatii

m

0

n

k

xk∑=

⎛⎜⎜⎝

⎞⎟⎟⎠

⋅ b

0

n

k

1∑=

⎛⎜⎜⎝

⎞⎟⎟⎠

⋅+

0

n

k

yk∑=

m

b⎛⎜⎝

⎞⎟⎠

MinErr m b,( ):= m 22.029= b 6.095= \\ rezultate numerice

⎝ ⎠f x( ) m x⋅ b+:= \\ forma analitica si grafica (mai jos) a functiei de aproximare

f x( ) 22.029x⋅ 6.095+

0 1 2 3 4 5 6

50

100

150

f x( )

y k

x xk,

Functie de aproximare generalizata

f x( ) a0 f0 x( )⋅ a1 f1 x( )⋅+ ..+ an fn x( )⋅+= \\ combinatie liniara de functii de aproximare;

\\ functiile de aproximare adoptate;

f t( )

t2

t

1

0

⎛⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎠

:=

x

0

1

2

3

4

5

⎛⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎠

:= y

0

4

9

12

25

33

⎛⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎠

:=

\\ perechile de puncte de aproximare;

coef linfit x y, f,( ):= \\ procedura de aproximare a coeficientilor combinatiei liniare de functii alese;

g t( ) coef f t( )⋅:= coef

0.929

1.957

0.429

0

⎛⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎠

=

0 1 2 3 4 5

10

20

30

40

y

g t( )

x t,

Observatie: functia de aproximare generalizata se preteaza la calcule si reprezentari grafice cu dificultate sporita.

Aplicaţia 2

1. Se cunoaste un sir de valori ale sarcinii electrice de pe un fir intr-un plan cu sistem de coordonate adoptat, in puncte de abscise cunoscute. Care este, in aceasta situatie, functia numerica de distributie de sarcina electrica pe fir in planul considerat? Datele problemei:

-- coeficient de multiplicare pentru afisare grafica optima m 106:=

x

1

1.2

2

2.5

3.4

6

7.6

8

9

10

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠

:= mm Q x( )

0.5

0.6

0.8

1.4

2.4

3.5

5.7

6

8

20

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠

10 6−⋅ m⋅:= Coulomb

-- lungimea firului: L 15:= mm

In aceasta problema se urmareste exprimarea functiei numerice dupa care variaza sarcina unui corp bidimensional, cand sunt determinate, posibil experimental, valori ale sarcinii in puncte diferite de pe corp (fir). Astfel, se ajunge la densitatea lineica de sarcina pentru corpul considerat, pe baza careia se pot calcula intensitatea campului electric, potentialul electric in orice punct din plan.

. Modelul matematic la care se reduce problema:

Q x( ) Functie x( ) dQ ρL x( ) dx⋅ ρL x( ) Functie1 x( )

-- intervalul ales: L 0 15..:= mm

-- numarul de valori stabilite: n 10:= i 0 n 1−..:= j 0 n 1−..:=

In vederea calcularii functiei de repartitie se apeleaza la metoda numerica de interpolare Lagrange (interpolare polinomiala).-- metoda consta in gasirea unui polinom, cu coeficienti reali, a carui valori in abscisele cunoscute sa fie egale cu valorile sarcinii in acele puncte, iar pe intervalul de lungime a firului functia de repartitie sa fie egala cu polinomul respectiv;-- intervalul considerat este lungimea firului cu sistemul de coordonate fixat la inceputul acestuia; -- polinomul este unic pentru setul de valori prescris;

-- valorile absciselor dispuse orizontal (se apeleaza din meniul Vector and Matrix Palette, obiect Matrix Transpose, ori Ctrl+1):

xT 0 1 2 3 4 5 6 7 8 9

0 1 1.2 2 2.5 3.4 6 7.6 8 9 10=

-- valorile ordonatelor (sarcinii): y x( ) Q x( ):=

y x( )T 0 1 2 3 4 5 6 7 8 9

0 0.5 0.6 0.8 1.4 2.4 3.5 5.7 6 8 20=

-- coeficienti de calcul cu instructiunea de conditionare 'if' (daca i=j, valoarea coeficientului este 1):

coef ij

if i j( ) 1, xi xj−( ),⎡⎣ ⎤⎦⎡⎣ ⎤⎦∏:=

-- functie intermediara necesara rularii metodei (polinom de ordinul n-1, obtinut prin inmultire):

intermediar z( )

j

z xj−( )∏:= z este variabila functiilor intermediare si finala, de interpolare;

-- formula definita a polinomului initial Lagrange:

l i z,( ) if z xi 1,intermediar z( )

z xi−( ) coef i⋅,⎡

⎢⎣

⎤⎥⎦

:=

L z( )

i

l i z,( ) y x( )i⋅∑:=-- polinomul de interpolare, numit Lagrange:

-- forma functiei de interpolare pentru lungimea corpului incarcat cu sarcina:

z 0 0.1, 10..:=

0 2 4 6 8 10

10

20

L z( )

y x( )

z x,

D. Determinarea functiei de interpolare si reprezentarea ei grafica este posibilasi prin utilizarea unor functii predefinite in Mathcad, care metode au la baza algoritmi interni precum cel prezentat anterior. In continuare, se exemplifica o functie predefinita de interpolare din butonul 'Insert Function', linterp(x,y,z).

B z( ) linterp x y x( ), z,( ):= z 0 0.001, 10..:=

0 2 4 6 8 10

10

20

B z( )

y x( )

z x,

Observatie: Acuratetea de reprezentare este evidenta, in cazul graficelor de mai sus. Astfel, cand se cere o aproximatie cat mai apropiata de adevar metoda numerica expusa se dovedeste mai avantajoasa, atat pentru figurarea variatiecat si pentru calculul in puncte aflate in domeniul celor cunoscute.

-- in vederea exemplificarii ultimului argument se stabileste un set de valori de pe lungimea corpului pentru care se calculeaza valorile sarcinii:

k 0 3..:=

X

1.5

3

7

9.5

⎛⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎠

:= L Xk( )0.5912.051

5.044

11.776

= B Xk( )0.6751.956

4.875

14

=

metoda numerica functia predefinita

Obtinerea cat mai exacta a functiei de interpolare depinde si de gradul polinomului de interpolare; daca acesta este inferior si precizia va fi de nivel scazut; cresterea gradului este insa limitata si superior.

Aplicaţia 3

Intre doi electrozi cilindrici coaxiali de lungime egala, perfect conductori, este dispus un material de conductivitate σ data. Pentru lungimea cilindrilor si pentru raza cilindrului exterior se dau cateun set de valori in functie de raza cilindrului interior. Se cere determinarea valorii rezistentei dispozitivului coaxial pentru orice valoare a razei cilindrului interior pe baza variatiei seturilor de date furnizate.

σ 5.7 107⋅:=

1Ω m⋅

(conductivitatea materialului)

a

0.1

0.2

0.35

0.65

1

⎛⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎠

:= m L

10.05

9.40

8.98

9.96

14.04

⎛⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎠

:= m b

2.20

3.42

5.26

9.06

13.72

⎛⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎠

:= m

A-B. Rezistenta electrica a dispozitivului poate fi privita ca o suma de rezistente infinitezimale de forma cilindrica cu peretii de grosime dr, unde r este raza unui contur cerc situat intre cei doi cilindri supraconductori. Astfel, redarea analitica a dependentei de marimile geometrice a rezistentei totale rezulta prin integrarea rezistentelor infinitezimale considerate a fi in serie:

R

a

b

r1

σ S⋅

⌠⎮⎮⌡

d R

a

b

r1

2 π⋅ r⋅ L⋅ σ⋅

⌠⎮⎮⌡

d R1

2 π⋅ L⋅ σ⋅ln

ba

⎛⎜⎝

⎞⎟⎠

⋅

C. Deoarece parametrii L si b sunt determinati de parametrul a, variatia rezistentei se va face dupa o functie compusa din combinatia dependentelor celor doi parametrii de raza cilindrului mic, deci va avea o forma analitica necunoscuta. Corespondenta rezistenta electrica - raza cilindru interior se va exprima prin interpolare, pe baza datelor din cele 3 seturi: a, L si b.

L f a( ) b g a( ) => R1

2 π⋅ σ⋅ f a( )⋅ln

g a( )a

⎛⎜⎝

⎞⎟⎠

⋅

--> setul de valori ale rezistentei dupa parametrul geometric a, in conditiile in care s-au inglobat si ceilalti parametri in formula arata in modul urmator:

i 0 last b( )..:= --> numarul de indici ale valorilor din set

--> calculul individual al fiecarei valori Ri

12 π⋅ σ⋅ Li⋅

lnbi

ai

⎛⎜⎜⎝

⎞⎟⎟⎠

⋅:= R σ⋅

0.04895

0.04807

0.04803

0.0421

0.02969

⎛⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎠

=

--> in acest stadiu se aplica metoda lui Newton cu diferente divizate de interpolare a functiei data in puncte .

--> algoritmul de scriere a metodei face apel la paleta de programare din Mathcad, se prezinta mai jos si este explicat pentru fiecare linie:

--> se incarca intr-un vector intern sirul valorilor lui a; --> se incarca in alt vector sirul valorilor R--> se stabileste numarul de linii si coloane ale matricei cu diferente divizate, n;

--> prima coloana a matricei primeste valorile vectorului Y incarcat anterior;

--> printr-un calcul recursiv se ocupa toate elementele superioare diagonalei principale cu diferentele divizate de ordin 1 pe a doua coloana, de ordin 2 pe a 3, s.a.m.d.

--> celorlalte elemente, de sub diagonala principala, li se atribuie valoarea 0.

--> algoritmul returneaza matricea Matrice_dif_div completata.

Matrice_dif_div X a←

Y R←

n last X( )←

A j 0, Y j←

j 0 n..∈for

A k j,

A k 1+ j 1−, A k j 1−,−

Xk j+ Xk−←

k 0 n j−..∈for

A n j− i+ j, 0←

i 1 j..∈for

A

j 1 n..∈for

A

:=

A Matrice_dif_div:= --> redenumirea matricei

A

8.588 10 10−×

8.433 10 10−×

8.426 10 10−×

7.386 10 10−×

5.208 10 10−×

1.546− 10 10−×

4.711− 10 12−×

3.467− 10 10−×

6.222− 10 10−×

0

5.995 10 10−×

7.6− 10 10−×

4.238− 10 10−×

0

0

2.472− 10 9−×

4.202 10 10−×

0

0

0

3.214 10 9−×

0

0

0

0

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠

=

--> functia de interpolare Newton apeleaza elementele matricei, care constituie diferedivizate, construind astfel aproximarea:

j 1 last a( )..:=

N t( )

j

A0 j,0

j 1−

i

t a( )i−⎡⎣ ⎤⎦∏=

⋅⎡⎢⎢⎣

⎤⎥⎥⎦

∑ A0 0,+:=last a( ) 4=

--> reprezentarea grafica a variatiei rezistentei dupa parametrul a se figureaza mai josdomeniul limitat inferior si superior de minimul, respectiv maximul sirului dimensiunii a

0.2 0.4 0.6 0.8 10.02

0.03

0.04

0.05

N t( ) σ⋅

Rσ⋅

t a,

D. Varianta oferita de Mathcad pentru identificarea numerica si grafica a variatiei rezistentei electrice in functie de dimensiunea parametrului a, implica apelarea functiei de interpolare linterp cu argumente: sirul valorilor lui a, sirul valorilor lui R inmultit cu σ si o variabila notata:

f z( ) linterp a R σ⋅, z,( ):= --> denumirea functiei de interpolare cu variabila notata z;

--> in continuare se reprezinta grafic dependenta rezistentei 'R' de 'a':

0.2 0.4 0.6 0.8 10.02

0.03

0.04

Variatia rezistentei dupa parametrul a

Rσ⋅

f z( )

N z( ) σ⋅

a z,

R dim_a( ) N dim_a( ):=

E. In grafic se compara caracteristica trasata de functia de interpolare Newton cu diferente divizate si caracteristica trasata de functia linterp de interpolare liniara; ultima functie uneste punctele prin segmente de dreapta, ceea ce reduce precizia de calcul. Avand exprimata functia de interpolare Newton cu diferente divizate este posibila calcularea rezistentei electrice pentru orice valoare numerica a dimensiunii 'a'. Se dau cateva exemple:

a1 0.25:= m R a1( ) 8.429 10 10−×= Ω

a2 0.47:= m R a2( ) 8.249 10 10−×= Ω

Aplicaţia 4

x

1

5

9

12

20

⎛⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎠

:= U x( )

6−

1

4

13

28

⎛⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎠

:= n 4:= i 0 n..:= j 0 n..:=

dij

if i j( ) 1, xi xj−( ),⎡⎣ ⎤⎦⎡⎣ ⎤⎦∏:= ω z( )

i

z xi−( )∏:= l i z,( ) if z xi 1,ω z( )

z xi−( ) di⋅,⎡

⎢⎣

⎤⎥⎦

:= L z( )

i

l i z,( ) U x( )i⋅∑:=

z 0 10..:= L 10( ) 6.277= L 8( ) 2.468=z01

2

3

4

5

6

7

8

9

10

= L z( )-12.158

-6

-2.262

-0.243

0.668

1

1.19

1.591

2.468

4

6.277

=

z 1 20..:= \\ pentru reprezentarea grafica;

0 5 10 15 2020

0

20

40

U x( ) i

L z( )

xi z,

1. Se dau urmatoarele rezultate ale unor masuratori, in functie de distanta, ale valorilor tensiunii electrice. Se cer valorile polinomului de interpolare Lagrange, de un grad fixat in jurul unui punct z in m+1 puncte, acest polinom aproximand functia. Sa se reprezinte grafic punctele masurate si polinomul de interpolare si sa se calculeze valoarea lui in z=1100 cu un polinom de gradul doi.

Aplicaţia 5

2. Intr-un mediu de dimensiuni foarte mari, se masoara in puncte spatiale egal distantate consecutiv, pe o linie dreapta pornind de la o origine bine fixata, un set de valori ale potentialului electric, creat de un camp electric static. Pe baza acestor masuratori si a faptului ca starea electrica a mediului se considera a avea simetrie, se urmareste determinarea distributiei de potential in intreg mediul. Prin ce modalitate este posibila realizarea acestei optiuni? Datele problemei:

-- sirul de valori spatiale cunoscute:

0 0, 0,( ) 1 1, 1,( ) 2 2, 2,( ) 3 3, 3,( ) 4 4, 4,( ) 5 5, 5,( )

-- starea numerica a potentialului in punctele de masura in ordinea definirii lor:

Vp 0 0.5 0.7 0.86 1 0.86( )V:=

D. Solutionarea problemei se poate realiza si prin apelul la o functie predefinita in Mathcad; aceasta functie, numita 'cspline' are ca algoritm intern o metoda de interpolare polinomiala, la care precizia de calcul este, dupa cum va fi aratat, diferita de metoda utilizata mai sus.

t

0

1

2

3

4

5

⎛⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎠

:= w

0

0.5

0.7

0.86

1

0.86

⎛⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎠

= S cspline t w,( ):= g r( ) interp S t, w, r,( ):=

-- reprezentarea grafica a functiei:

0 1 2 3 4 5

0.5

1

1.5

g r( )

r

var 0 0.5, 5..:= --> esantion de valori ale variabilei spatiale

f var( )0

0.312

0.5

0.616

0.7

0.778

0.86

0.942

1

0.993

0.86

= g var( )0

0.307

0.5

0.618

0.7

0.777

0.86

0.945

1

0.985

0.86

=

--> valorile functiei definita dupa metoda

valorile functiei definita intern in Mathcad <--

Reprezentarea grafica tridimensionala a distribuitiei de potential electric:

Ai j, f xi( ) sin rot j( )⋅:= Bi j, f yi( ) cos rot j( )⋅:= Ci j, f zi( ):=

A B, C,( )

Aplicaţia 6

Functii Mathcad pentru interpolarea unei suprafete in 2D

Se introduce o matrice patratica care genereaza o suprafata:

Mz

2

4.7−

12−

2

9

0

0.9−

2

5

0.1

8

0.8

5

8.9

17

0.8−

9.2

9.6−

7

9

12−

3.4

1−

3

0.9

7.1

1

5

9

2

3

12

0

7

3.2

0

⎛⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎠

:=

rows Mz( ) 6=

cols Mz( ) 6= n rows Mz( ):= n 6=

Se scriu vectorii X si Y cu n-linii care determina reteaua:

X

0

1

2

3

4

5

⎛⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎠

:= Y

0

1

2

3

4

5

⎛⎜⎜⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎟⎟⎠

:=

Mxy augment sort X( ) sort Y( ),( ):= rows Mxy( ) 6=

Coeficientii functiei spline

S cspline Mxy Mz,( ):= fit x y,( ) interp S Mxy, Mz,x

y⎛⎜⎝

⎞⎟⎠

,⎡⎢⎣

⎤⎥⎦

:=

Suprafata originala realizata doar prin unirea cu segmente a punctelor matricei:

2D Spline- Suprafata interpolata in care punctele intermediare sunt calculate si afisate:

Mz FIT