Metode Numerice de Rezolvarea Ecuaiilor Difereniale

.l. Dr. ing. Levente CZUMBIL E-mail: Levente.Czumbil@ethm.utcluj.ro WebPage: http://users.utcluj.ro/~czumbil

mailto:Levente.Czumbil@ethm.utcluj.rohttp://users.utcluj.ro/%7Eczumbil

Generaliti

O ecuaie diferenial este o ecuaie care conine pe lng variabilele independente ifunciile necunoscute i derivatele acestor funcii. Ordinul ecuaiei difereniale, n, este datde ordinul maxim al derivatelor funciei necunoscute din cadrul aceste ecuaii.

( ) 0,...,'',',, )( =nyyyyxfEcuaiile difereniale cu derivate pariale conin mai multe variabile independente

i derivatele pariale ale funciilor necunoscute.

0,,,,,,, 222

2

2

=

yz

yxz

xz

yz

xzzyxy

Rezolvarea unei ecuaii difereniale de ordin n implic impunerea a n condiii iniiale.

Se consider un circuit format dintr-un rezistor de rezisten R i o bobin de inductivitate L,alimentate n serie la o tensiune electromotoare .tEe cos=

Circuitul R-L serie n regim tranzitoriu

Se studiaz variaia curentului n circuit la nchiderea ntreruptorului K.

tcosEiRdtdiL

dtdiLiReee LR =++=+=

)t(iR)t(udt

)t(diL)t(u

R

L

=

=

tEtiRdt

tdiL cos)()( =+

Se scriu teoremele lui Kirchhoff i rezult o ecuaie diferenial de ordinul I:

innd cont de dependena de timp (t) : unde i(t) - curentul din circuit la momentul t

uL,R(t) - tensiunea la bornele bobinei respectiv rezistenei la momentul t ecuaia diferenial se va rescrie :

Circuitul R-L Serie

Fie funcia, f(x):I x RR unde I este un interval real, f fiind o funcie continu dat, iary0 fiind valoarea iniial a acesteia.

Evaluam funcia y:I R n nodurile intervalului de definiie, funcie care satisfaceproblema cu condiia iniial Cauchy impus:

bxxxa n

iar: ),(),(),(),( )1()1()( yxfyxfyxfyxf nyn

xn +=

))x(y,x(f)x(y )n()n( =+1

Pentru n=2 rezult:

( ) ( ) ( ) ( )( )

+++= 0000000001 ,,,2

, yxfyxfyxfhyxfhyy yx

nyy ,...,2Prin recuren

))(2

(1 yxii fffhfhyy +++=

unde:yxff

yff

xff

yff

xff xyyyxxyx

=

=

=

=

=2

2

2

2

2

; ; ; ;

Fie circuitul R-L serie din cadrul aplicaiei prezentate pentru care avem cunoscuteparametrii electrici: E = 12V, R = 4 i L = 3,2uH. S se determine curentul prin bobina deinductivitate L dup nchiderea ntreruptorului K (t ia valori pe intervalul [0;40ms])

Pasul 1. Se definesc parametri electric ai circuitului R-L serie:

Pasul 2. Se scrie ecuaia diferenial ce descrie funcionarea circuitului R-L serie:

Pasul 3. Se extrage derivata curentului din ecuaia diferenial corespunztoare circuitului:

Pasul 4. Se definete funcia asociat membrului drept a ecuaiei difereniale:

F t i, ( )1

LE cos t( ) R i( ):=

E 12:= R 4:= L 32 10 6:= f 50:= 2 f:= 314.159=

Pasul 5. Se definesc capetele intervalului, numrul de puncte de calcul i se determin pasul de parcurgere al intervalului de definiie:

Pasul 6. Se determin irul de puncte intermediare tk n care se evalueaz valoareacurentului:

ti 0:= tf 40 103

:= N 500:= htf ti

N:= h 8 10 5=

Pasul 7. Se definesc derivatele pariale ale funciei ataate, F, ecuaiei difereniale:

k 0 N..:= tk ti h k+:=

Pasul 8. Din condiia iniial Cauchy a problemei (ntreruptorul K deschis), reiese cvaloarea curentului n momentul t = 0s este egal cu 0A:

Ft t i, ( ) tF t i, ( )d

d:= Fi t i, ( ) i

F t i, ( )dd

:=

I0 0:=

Pasul 9. Se implementeaz formula recursiv de calcul a valorilor funciei, pe bazadescompunerii n serie Taylor pn la elementul de gradul al II-lea:

Ik 1+ Ik h F tk Ik, ( )h2

Ft tk Ik, ( ) F tk Ik, ( ) Fi tk Ik, ( )+( )+

+:=

Pasul 10. Se vizualizeaz valoarea curentului la momentele de timp tk:

0 0.01 0.02 0.03 0.043

2

1

0

1

2

3

I

t

)884.1848.1811.1185.0094.00( =TI

Se obine din metoda Taylor pentru n=1, adic se rein numai primii doi termeni dindezvoltare rezultnd forma explicit a metodei lui Euler:

),( 111 += iiii yxfhyy ,...2,1=i

Avem bineneles aceeai problem de rezolvare a ecuaiilor difereniale cu condiiiiniiale:

),(' yxfy = 00 )( yxy =

Pentru variante ale metodei lui Euler cu precizie mai mare se folosesc relaii derecuren de forma:

),,( 111 hyxhyy iiii +=unde n:

a) Metoda lui Euler mbuntit (Euler Heun):

[ ])',(),(21),,( 1111111 +++= iiiiiii yhyhxfyxfhyx

),(' 111 = iii yxfycu

b) Metoda lui Euler modificat (versiunea Cauchy)

++= 11111 '2

,2

),,( iiiii yhyhxfhyx ( )111 ,' = iii yxfy

c) Metoda Euler modificat predictor corector (metoda implicit).

Cu metoda lui Euler clasic se calculeaz o prim aproximaie (valoarea prezis a soluiei n punctul urmtor) adic se iniializeaz valoarea lui yi cu o relaie:

( ) ( )1110 , += iiii yxfhyyLa un pas oarecare al procesului iterativ de calcul noua valoare a lui rezult prin aplicarea

unei relaii:( ) ( ) ( )( )

2,, 111

1

+

+=k

iiiii

ki

yxfyxfhyy

Calculul se consider terminat cnd se atinge o precizie impus aprioric:

( ) ( )

Fie ecuaia diferenial de ordinul nti cu condiia

iniial Cauchy y(7)=6, unde x ia valori pe intervalul [7,15]. S se determine valorilefunciei y(x) folosindu-se metoda lui Euler mbuntit (Euler-Heun), respectiv variantamodificat (versiunea Cauchy).

y ' x( ) 6y x( ) co s x

9

+ x2 5

Pasul 1. Se scrie ecuaia diferenial ce urmeaz a fi rezolvat:

Pasul 2. Se extrage derivate funciei necunoscute:

y' x( ) x2 5 6 y x( )+ cos x

9

Pasul 3. Se definete funcia asociat ecuaiei difereniale:

f x y, ( ) x2 5 6y+ cos x

9

:=

Pasul 4. Definirea funciei caracteristice metodei mbuntite Euler-Huen:

EH x y, h, ( )12

f x y, ( ) f x h+ y h f x y, ( )+, ( )+( ):=

Pasul 5. Definirea funciei caracteristice metodei mbuntite Euler-Huen:

EC x y, h, ( ) f xh2

+ yh2

f x y, ( )+,

:=

Pasul 6. Se definesc capetele intervalului, numrul de puncte de calcul i se determin pasul de parcurgere al intervalului de definiie:

a 7:= b 15:= N 100:= hb a

N:= h 0.08=

Pasul 7. Se determin irul de puncte intermediare xi n care se evalueaz valoareafunciei necunoscute:

i 0 N..:= xi a h i+:=

Pasul 8. Se impune condiia iniial Cauchy y(7)=5:

Pasul 9. Se evalueaz valorile funciei necunoscute conform metodei lui Eulermbuntite (Euler-Heun) :

Pasul 11. Se vizualizeaz valorile funciei necunoscute determinate n punctele xi:

yEH05:= yEC0

5:=

yEHi 1+yEHi

h EH xi yEHi, h,

+:=

Pasul 10. Se evalueaz valorile funciei necunoscute conform metodei lui Eulermodificat (versiunea Cauchy) :

yEC i 1+yEC i

h EC xi yEC i, h,

+:=

yEHT 0 1 2 3 4 5 6

0 5 12.295 24.046 42.911 73.122 121.436 ...=

yECT 0 1 2 3 4 5 6

0 5 12.294 24.046 42.91 73.12 121.433 ...=

Pasul 12. Se reprezint grafic alura funciei determinate cu cele dou metode:

6 8 10 12 14 160

1 1021

2 1021

3 1021

yEH

yEC

x

Pasul 13. Se evalueaz abaterea procentual dintre cele dou metode:

Err1N

i

yEHiyECi

yEHi:= Err 2.412 10 3 %=

Metodele Runge Kutta de integrare numeric a unei ecuaii difereniale, nlocuiesccalculul derivatelor funciei f(x,y) prin evaluri ale sale n diverse puncte.

Fie ecuaia diferenial cu condiii iniiale de forma: ( )( )

==

00

,'yxyyxfy

unde: f(x):DR, D R2, este o funcie cu derivatele pariale continue pe D, undei+j=k, respectiv . ji

k

yxf

mk ,1=

Soluia calculndu-se cu o relaie unipas de forma:

nnii kakakayy ++++=+ ...11001Din condiia ca dezvoltarea n serie Taylor a membrului drept (n funcie de h) s coincidcu membrul drept al formulei lui Taylor de ordinul (n+1) avem:

( )( )( )

( )11,1100

12102022

01011

0

...,...............................................................

,,

,

+++++=

+++=++=

=

nnnnninin

ii

ii

ii

kkkyhxfhk

kkyhxfhkkyhxfhk

yxfhk

Particulariznd parametrul n determinm diverse formule de tip Runge Kutta:

a) n=0 (Runge-Kutta de ordin II), y0 dat (formula lui Euler clasica):

( )iiii yxfhyy ,1 +=+b) n=1 (Runge-Kutta de ordin II), y0 dat:

( )101 21 kkyy ii ++=+ ( )ii yxfhk ,0 = ( )01 , kyhxfhk ii ++=

( ) ( )( )[ ]iiiiiiii yxfhyhxfyxfhyy ,,,21

++++=+

c) n=2 (Runge-Kutta de ordin III), y0 dat:

( )2101 461 kkkyy ii +++=+ ( )ii yxfhk ,0 =

++=

2,

20

1kyhxfhk ii

( )012 2, kkyhxfhk ii ++=

S se rezolve ecuaia diferenial de ordinul I cu condiia iniial Cauchy

y(0)=0 pe intervalul [0,2] cu pasul h=0.01 folosind metoda Runge-Kutta de ordinul III.yxx

yxy++

= 2

22

14'

Pasul 1. Se definete funcia asociat membrului drept a ecuaiei difereniale:

f x y, ( )4 x2 y2

1 x2+ x y+:=

Pasul 2. Se definesc capetele intervalului de definiie i se determin numrul de puncte de calcul:

a 0:= b 2:= Nb a

h:=

200N =

Pasul 3. Se determin vectorul al punctelor din intervalul n care se calculeaz valorilefunciei :

i 0 N..:= xi a h i+:=

h 0.01:=

N

b

a

-

h

:=

200

N

=

a

0

:=

b

2

:=

i

0

N

..

:=

x

i

a

h

i

+

:=

y0 0:=

Pasul 5. Se definesc coeficienii Runge-Kutta de ordinul III:

Pasul 4. Din condiia iniial Cauchy a problemei, reiese c funcia y(x) ia valoarea 1 n punctul 0 (n captul a al intervalului):

k0 x y, ( ) h f x y, ( ):=

k1 x y, ( ) h f xh2

+ yk0 x y, ( )

2+,

:=

k2 x y, ( ) h f x h+ y 2 k1 x y, ( )+ k0 x y, ( ), ( ):=Pasul 6. Se implementeaz formula recursiv Runge-Kutta pentru calculul valorilor funciei n punctele :

yi 1+ yi16

k0 xi yi, ( ) 4 k1 xi yi, ( )+ k2 xi yi, ( )+( )+:=

Pasul 7. Se vizualizeaz valorile funciei y(x) pe intervalul [0,2] i se reprezint grafic:

( )45.1453.1237.0198.0159.012.008.004.00 =Ty

Metode Numerice de Rezolvarea Ecuaiilor Difereniale

.l. Dr. ing. Levente CZUMBIL E-mail: Levente.Czumbil@ethm.utcluj.roWebPage: http://users.utcluj.ro/~czumbil

mailto:Levente.Czumbil@ethm.utcluj.rohttp://users.utcluj.ro/%7Eczumbil

Slide Number 1Slide Number 2Slide Number 3Slide Number 4Slide Number 5Slide Number 6Slide Number 7Slide Number 8Slide Number 9Slide Number 10Slide Number 11Slide Number 12Slide Number 13Slide Number 14Slide Number 15Slide Number 16Slide Number 17Slide Number 18Slide Number 19Slide Number 20