Curs 9. Derivarea numeric

Calculul aproximativ al derivatelor

cu aplicaii n ingineria electric

METODE NUMERICE

n domeniul ingineriei electrice exist situaii practice cnd este necesar evaluarea numeric a valorilor derivatelor unor funcii pentru care derivarea analitic este greoaie sau chiar imposibil

Trebuie inut cont c n tehnic de obicei nu se cunoate expresia analitic a funciei care trebuie derivat, ci doar valorile ei n anumite puncte (determinate experimental sau prin calcule).

Se dorete determinarea aproximativ a derivatei n punctele unde se cunoate valoarea funciei, ct i n alte puncte.

Exemple:

determinarea distribuiei de sarcindeterminarea potenialelor curbelor de sarcindeterminarea consumului de energie pe baza curbelor de sarcindeterminarea valorilor maxime ale puterilor consumate de un receptorcalculul intensitii cmpurilor electromagnetice, etc.

Exemplu practic: prelucrarea curbelor de sarcin

012345

0 2 4 6 8 10 12 14 16 18 20 22 24t

P

Curba de sarcin activ zilnic a unui consumator

Se consider un receptor de energie electric pentru care se cunoate curba de sarcin zilnic referitoare la puterea activ consumat.

Derivarea numeric bazat pe polinoame de interpolare

Pentru simplitate vom considera doar derivata de ordinul I. Se pot aplica tehnici analoage i pentru derivate de ordin superior.

Vom rezolva problema prin interpolare: functia dat tabelar o determinm prin interpolare si vom deriva polinomul su de interpolare!!!

.

De regul se aproximeaz derivatele n punctele x0, x1,,xn (noduri) n care valoarea funciei este cunoscut, dar se pot utiliza i pentru alte puncte din intervalul [a,b]!!!

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

( ) ( ) ( )( ) ))((!2

''

))((!2

''

000

00

0

10110011

hxxxxxfhxxhxf

hhxxxf

xxxxxfxlxfxlxfxRxLxf

+++=

=++=+=

( ) [ ]b,ax

Derivarea numeric bazat pe polinomul de interpolare LagrangeInterpolarea lui f s-a facut prin polinomul de interpolare Lagrange i vom construi funcia pe baza polinomului Lagrange de ordin I i a erorii (funcia este cunoscut experimental n nodurile x0 i x1):

( ) ( ) ( )( )( )

( ) ( ) ( ) ( )( ) ( )( )( )dx

x''fd2

)hxx)(xx(x''f2

hxx2h

xfhxfdx

)hxx)(xx(!2x''fd

hxfhxfx'f

00000

0000

+++=

=

++=

( ) ( ) ( )h

xfhxfx'f 00 +

Derivnd se obine:

Deci:

( ) [ ] hxxbaCfbax += 0120 ,,,,

( )( )( )dx

x''fd

( ) ( ) ( ) ( )++= ''f2h

hxfhxfx'f 000

Un inconvenient al formulei este c nu avem informaii despre termenul

Cnd x=x0 coeficientul termenului problem este zero deci formula se simplific :

deci eroarea de trunchiere nu poate fi estimat.

( ) ( )h

xfhxf 00 +

( ) Mx''f,M2h

Er [ ]b,ax

0h >0h calculul derivatei cu operatorul definit in Mathcad;

Dv4 2.048 104= x4 5.1= y4 1.711 10

4=

j 0 N 1..:= Def j f' xj( ):=xT 0 1 2 3 4 5 6 7

0 5 5.025 5.05 5.075 5.1 5.125 5.15 5.175=

yT 0 1 2 3 4 50 1.51610 4 1.56310 4 1.61110 4 1.6610 4 1.71110 4 1.76210 4

=

DvT 0 1 2 3 4 50 1.85310 4 1.910 4 1.94810 4 1.99810 4 2.04810 4 2.09910 4

=

DT 0 1 2 3 4 50 1.85310 4 1.910 4 1.94810 4 1.99810 4 2.04810 4 2.09910 4

=

DefT 0 1 2 3 4 50 1.85310 4 1.910 4 1.94810 4 1.99810 4 2.04810 4 2.09910 4

=

eri f' xi( ) Dvi( ):= eroare max er( ):=eroare 5.815 10 6= \\ referinta considerat operatorul

Mathcad.

er_reli f' xi( )( ):= er_rel min er_rel( ):= eroare_rel eroareer_rel:=\\ eroarea relativa arata pana la a catea zecimala rezultatul este corect.

eroare_rel 3.139 10 8 %=

Derivarea numerica utilizand dezvoltarea in serie Taylor (derivata de ordinul I, pentru functii date in puncte)

Fie un set de puncte N, cu pasul h intre ele, in care se cunosc valorile unei functii, y:

N 102:= h 10 2:= i 0 N..:=

x0 1:= xi 1+ xi i h+:= y0 20:= yi 1+ yi i2 h:=

Aproximarea numerica a derivatei

der1 i( ) Ai1

i 1+( ) h yi 1+ yi 1( )Ai

:= \\ elementele unui vector se incarca cu relatia dependenta de variabila i, care schimba nodurilesi y

der

Bj der1 j( )j 1 N..for

B

:= \\ elementele individuale ale vectorului anterior se stocheaza in alt vector, numit der

length der( ) 101=length y( ) 102=

Valorile numerice obtinute sunt:

x

0

01

2

3

4

5

6

7

8

9

10

11

11

1.01

1.03

1.06

1.1

1.15

1.21

1.28

1.36

1.45

1.55

= y

0

01

2

3

4

5

6

7

8

9

10

11

2020

19.99

19.95

19.86

19.7

19.45

19.09

18.6

17.96

17.15

16.15

= der

0

01

2

3

4

5

6

7

8

9

10

11

0-0.5

-1.667

-3.25

-5

-6.833

-8.714

-10.625

-12.556

-14.5

-16.455

-18.417

=

Sa se determine numeric campul electric pe suprafata unui cilindru, dupa o directie in situatia in care se cunoaste expresia analitica a fluxului electric care strabate suprafata si coordonatele cilindrice care o definesc. Date numerice:

-- campul scalar flux electric: x y, z,( ) 1012 x2 z 3x( ) y 3z:= Wb

-- raza cilindrului: r 5:= cm 2:=

-- inaltimea cilindrului: ab 6:= cm a 0:= b 6:= cm

(coordonatele cilindrice cunoscute: r, , z=ab)

A. Avand data o forma analitica pentru fluxul electric, se poate determina campul electric care produce acel flux prin derivare numerica dupa directia z, cand valorile x si y ale suprafetei si, implicit, ale campului de flux sunt impuse de raza cilindrului , r, si unghiul din coordonatele cilindrice. -- fluxul electric el -- campul electric (camp vectorial) Ev-- elementul de suprafata scris in coordonate cilindrice dS r dz

el SEv

d (se aplica operatorul de derivare)

EvS eld

d(care devine) Ev

1r z el

dd

B. Formula matematica care sta la baza rezolvarii problemei se exprima prin teorema lui Gauss:

-- dependenta variabilelor x si y de coordonatele cilindrice

x r cos ( ):= y r sin ( ):=

Pentru rezolvarea problemei se alege o metoda numerica de derivare utilizand dezvoltarea in serie Taylor. Aceasta metoda aproximeaza valorile derivatelor functiei intr-un sir de puncte in care valoarea functiei este cunoscuta.

-- sirul de puncte in care se evalueaza derivata se defineste in felul urmator:

m' 50:= (numarul de elemente din sir) i 0 m' 1..:= (numarul de iteratii)

a 0:= b 6:= (limitele intervalului din care face parte sirul)

hb a

m':= h 0.12= (pasul sirului, distanta dintre doua valori consecutive)

z0 a:= zm' b:= (valorile initiala si finala ale sirului)

zi 1+ zi h+:= (relatia generala intre elementele sirului)

Observatie: Derivarea numerica prin aceasta metoda poate conduce la erori semnificative, in anumite conditii, datorita imposibilitatii reducerii pasului h sub o anumita valoare limita.

-- functia matematica de derivat, dupa modelul completat este urmatoarea:

f z( ) x y, z,( )

r :=-- iar derivarea numerica are loc in cadrul formulei deduse din seria Taylor:

f' z i,( )f zi 1+( ) f zi 1( )

2h:=

-- versorul directiei campului electric dupa axa z: rz

0

0

1

:=

-- identificarea completa a punctelor in care se analizeaza campul electric, dupa axa z se face tinand cont de parametrii r si :

r 5= cm 7

:= x 5= cm y r sin ( ):= y 2.169= cm

-- expresia campului electric scrisa vectorial dupa directia z:

Ez z i,( ) f' z i,( ) rz:=

-- punctele sirului considerat in care se doreste aflarea campului (derivatei numerice):

zi z i,( ) (echivalenta)

z18 2.16= mm Ez z 18,( )0

0

0.003422280488564

= V

m

(valoarea campului electric)

Observatie: Metoda numerica expusa, pentru problema de fata se dovedeste eficace la un numar restrans de argumente ale functiei, adica pentru o valoare limitata de coordonate z pentru care se doreste cunoasterea marimii campului electric.

D. In continuare, se va apela operatorul predefinit de derivare din utilitarul Mathcad. Dupa aplicarea acestuia functiei, la valorile campului electric se ajunge prin inlocuirea coordonatei z cu punctele de interes.

z 2.16:= mm f' z( )zf z( )d

d:= f' z( ) 0.028=

-- valoarea campului electric calculata in acest fel: Ez z( ) f' z( ) rz:=

Ez z( )

0

0

0.028376055249965

= V

m

E. Raportarea celor doua valori obtinute, prin metoda numerica, respectiv prin functia definita in toolbox-ul Mathcad, evidentiaza acuratetea metodei numerice. Problema, prin referirea metodei numerice, face posibila exprimarea valorilor campului electric la suprafata de separatie a doua medii, dupa directia tangentiala. Prin alegerea corespunzatoare a valorilor r, si z se poate acoperi intreg domeniul suprafetei cilindrului, ceea ce inseamna determinarea valorilor campului electric pe toata suprafata considerata.

Un acumulator cu tensiune electromotoare si rezistenta interna date alimenteaza "n" rezistente egale, R. Cate rezistente R trebuie legate in serie si cate in paralel pentru ca in grupul de "n" rezistente astfel montate sa se dezvolte puterea maxima?

n 150:= R 20:= Ri 30:= Ue 100:= V

A. Puterea maxima obtinuta in rezistoare se determina efectuand derivata expresiei puterii dupa rezistenta echivalenta a circuitului de disipare. Matematic, in punctul in care derivata se anuleaza, exista un maxim pentru puterea disipata. B. Expresia analitica a puterii disipate are urmatoarea forma:

P n R Ip2 Ip --> curentul printr-un grup de rezistente conectate in serie.

--> se noteaza cu ns numarul de rezistente legate in serie si cu np

numarul de grupuri serie legate in paralel.

--> prin aplicarea teoremelor lui Kirchhoff pe unul din noduri si pe un ochi se deduce formula curentului electric printr-un grup de rezistente conectate in serie in functie de curentul electric total debitat de acumulator:

I np Ip 0 I Ri ns R Ip+ Ue --> dupa inlocuiri in formule se ajunge la expresia finala a puterii electrice debitate pe circuitul rezistiv:

Pnsnp

RUe

2

Rinsnp

R+

2 x

nsnp

--> notatie a raportului de numere intregi;

Rech x R --> rezistenta echivalenta a circuitului din rezistente R;

--> asadar, puterea electrica poate fi considerata ca o functie cu variabila rezistenta echivalenta direct dependenta de raportul x:

P x( ) x RUe

2

Ri x R+( )2:=

C. In vederea solutionarii problemei, se considera o metoda numerica pentru derivare, care foloseste aproximarea derivatei de ordinul I a functiei putere electrica cu o formula recurenta dupa 4 puncte vecine, cunoscute prin calcularea valorilor functiei intr-un set de noduri, care se atribuie a fi chiar valorile raportului x. Pentru aceasta se alege un interval posibil de situare a raportului, devenit variabila, x, un numar de noduri de calcul si o formula de salt in noduri:

a 0:= b 5:= --> capetele intervalului de situare pentru x;

N 102:= i 0 N..:= --> numarul de noduri de calcul;

hb a2 N:= xi a h i+:=

--> formula de salt in noduri dupa pasul h calculat;

--> formula de aproximare a derivatei:

Dvi

Ai

P xi h( ) 8 P xi h2 8 P xi h2+ + P xi h+( )6 h

n 1 N 1..for

Ai

:=

xP x( )

=>

--> se reprezinta grafic ajustat la intervalul de interes valorile calculate cu formula de aproximare a derivatei puterii:

--> cu optiunea Trace x-y din paleta grafica a utilitarului Mathcad se citeste valoarea x in care derivata se anuleaza:

x 1.5:= D. Utilitarul Mathcad ofera un operator pentru derivarea functiilor analitice, care va fi folosit in continuare spre a se evalua posibilele erori introduse de efectuarea derivarii numerice.

DP x( )x

P x( )

:= si reprezentarea grafica:

--> citirea graficului cu optiunea indicata ofera aceeasi informatie, adica x 1.5:=

Ramane de rezolvat sistemul de doua ecuatii cu doua necunoscute cu blocul Given - Find indicat intr-o problema anterioara:

ns 1:= np 1:=Given

ns np n

nsnp

x

ns

np

Find ns np,( ):= nsnp

15

10

= --> solutia problemei

E. Se propune o varianta grafic-numerica pentru stabilirea configuratiei de montare a rezistentelor intr-un circuit, astfel incat puterea sa fie maxima, fara a se continua mersul analitic al rezolvarii, pornit in paragraful B, dedicat modelului matematic al problemei. Varianta se justifica din punct de vedere al expresivitatii numerice directe a rezultatului.

Reprezentarea grafica a gradientului

Introducem functia f(x,y): f x y,( ) 2 sin x( )2 cos y 2( ):=

Introducem valorile de capat (punctele finale) pentru x si y:xmin 2:= xmax 2:= ymin 4:= ymax 4:=

Introducem numarul valorilor lui x, y in sir: xn 20:= y n 20:=

Gradientul functiei:

i 0 xn 1..:=grad x y,( ) x

f x y,( )dd

yf x y,( )d

d

:=

j 0 yn 1..:=

xindi xmin ixmax xmin

xn 1+:=

yind j ymin jymax ymin

yn 1+:=

Vi j, grad xindi yind j,( ):= \\ se stocheaza intr-o matrice valorile gradientului in punctele de calcul ale functiei;Mi j, Vi j,( )0:= Ni j, Vi j,( )1:= Fi j, f xindi yind j,( ):=

M N,( )Campul vectorial:

F

Reprezentarea grafica pe suprafata:

Jacobianul unei functii vectoriale

Fie o functie vectoriala de valori x, y si z: A x y, z,( )x z 2y+

y2 x+2 x3 y

:=

J A x, y, z,( )

xtr A x y, z,( )T( ) 0 d

d

xtr A x y, z,( )T( ) 1 d

d

xtr A x y, z,( )T( ) 2 d

d

ytr A x y, z,( )T( ) 0 d

d

ytr A x y, z,( )T( ) 1 d

d

ytr A x y, z,( )T( ) 2 d

d

ztr A x y, z,( )T( ) 0 d

d

ztr A x y, z,( )T( ) 1 d

d

ztr A x y, z,( )T( ) 2 d

d

T

:=

Jdet A x, y, z,( ) J A x, y, z,( ):= \\ determinantul matricei jacobian;

Evaluare simbolica din paleta Evaluation

J A x, y, z,( )z

2

x

1

2 y0

6 x2 y2 x3

0

Jdet A x, y, z,( ) x 2 x3 12 x2 y2+( )Evaluare numerica

J A 1, 3, 2,( )2

2

1

1

60

18

2

0

=

Jdet A 1, 3, 2,( ) 110=

Evaluarea numerica a derivatelor unei functii intr-un punct

Introducem functia pe care dorim sa o derivam:

f x( ) x3 ln x( ) 3 x2+ cos x( )2+:=Evaluarea simbolica a derivatei din paleta Evaluation:

xf x( )d

d3 x2 ln x( ) x2+ 6 x 2 cos x( ) sin x( )+

n 3:= nxf x( )d

d

n6 ln x( ) 11+ 8 cos x( ) sin x( )+

Punctul in care calculam derivata: x 4:=

Derivata de ordin n:Prima derivata:

n 3:= nxf x( )d

d

n23.275=

xf x( )d

d105.553=

Dupa copierea functiei de derivat, prin selectarea variabilei in functie de care se face derivarea, daca se trece in meniul Symbolics-Variable-Differentiatese genereaza expresia derivatei:

x3 ln x( ) 3 x2+ cos x( )2+ by differentiation

3 x2 ln x( ) x2+ 6 x 2 cos x( ) sin x( )+

Fie functia: f x( ) 3 x4 x2+ x+ 1+:= sa se evalueze derivatele in punctele (ca un sir):

n 1 5..:= D f w,( )w

f w( )dd

:= f n( )6

55

256

789

1.90610 3

= D f n,( )15

101

331

777

1.51110 3

=

D f w,( ) 12 w3 2 w+ 1+

Introducem punctele in care dorim calculul ca si componentele unui vector :

v

i

3 2i+sin 3 ( )

e2 0.123

ln e3( )1.123

:=

vT i 3 2i+ 0 535.492 0.123 3 1.123( )=n 0 rows v( ) 1..:= rows v( ) 7=

f vn( )3+1i

-348+374i

1

2.46710 11

1.139

256

8.155

= D f vn,( )1-10i

-101+556i

1

1.84310 9

1.268

331

20.241

=

Daca avem o functie de doua variabile si dorim calcularea derivatei partiale:

f u w,( ) 2 u3 w 5 sin u( ) cos 2 w( )+:=

Du f u, w,( )u

f u w,( )dd

:= Dw f u, w,( )w

f u w,( )dd

:=

Du f u, w,( ) 6 u2 w 5 cos u( ) cos 2 w( )+Dw f u, w,( ) 2 u3 10 sin u( ) sin 2 w( )

u 1 2..:= w 4 5..:= Du f u, w,( )23.60796.303

27.733

121.746

= Dw f u, w,( )-6.3257.004

6.578

20.947

=

Consideram o functie scalara: f x y, z,( ) x y x2 z3+ 3 y3 z:=

Introducem o functie sub forma unui vector: A x y, z,( )2 x2 y3

y2 z3 x4 z2

:=

Dorim calcularea simbolica si numerica a gradientului, divergentei si rotorului.

Gradientul: Grad f x', y', z',( )

x'f x' y', z',( )d

d

y'f x' y', z',( )d

d

z'f x' y', z',( )d

d

:=

Evaluare simbolica: Grad f x', y', z',( )y' 2 x' z'3+x' 9 y'2 z'

3 x'2 z'2 3 y'3

Evaluare numerica: Grad f 2, 1, 4,( )257

34189

= Daca Grad ,

atunci Div A .

Div A x', y', z',( )x'

tr A x' y', z',( )T( ) 0 dd y'

tr A x' y', z',( )T( ) 1 dd

+

z'tr A x' y', z',( )T( ) 2 d

d+

...

:=

Evaluare simbolica: Div A x', y', z',( ) 4 x' y'3 2 y' z' 6 x'4 z'+Evaluare numerica: Div A 2, 1, 4,( ) 384=

Daca Grad , atunci Curl A . (rotor)

Curl A x', y', z',( )

y'tr A x' y', z',( )T( ) 2 d

d z'tr A x' y', z',( )T( ) 1 d

d

z'tr A x' y', z',( )T( ) 0 d

d x'tr A x' y', z',( )T( ) 2 d

d

x'tr A x' y', z',( )T( ) 1 d

d y'tr A x' y', z',( )T( ) 0 d

d

:=

Evaluare simbolica: Curl A x', y', z',( )y'2

12 x'3 z'26 x'2 y'2

Evaluare numerica: Curl A 2, 2, 3,( )4

86496

=

B x' y', z',( ) f x' y', z',( ) A x' y', z',( ):= B x' y', z',( )2 x'2 y'3 x' y' x'2 z'3 3 y'3 z'+( )

y'2 z' x' y' x'2 z'3 3 y'3 z'+( )3 x'4 z'2 x' y' x'2 z'3 3 y'3 z'+( )

fA x' y', z',( )x'3 z'4 y'y'3 x'2 z'3

2 x'4 y'2 z'3

:= Div fA x', y', z',( ) 3 x'2 z'4 y' 3 y'2 x'2 z'3 6 x'4 y'2 z'2+

Curl fA x', y', z',( )4 x'4 y' z'3 3 y'3 x'2 z'2+4 x'3 z'3 y' 8 x'3 y'2 z'3

2 y'3 x' z'3 x'3 z'4