Curs_11_12_Cap9_Integr_EcDif.pdf

Dumitru Dragomir Metode numerice Curs 11, Curs 12 1

Cap. 9. Integrarea ecuaiilor difereniale

Metodele numerice pentru rezolvarea ecuaiilor difereniale pot fi grupate n dou categorii principale: a) metode directe de rezolvare a ecuaiilor difereniale; b) metode indirecte de rezolvare a ecuaiilor difereniale.

8.1. Metode numerice directe pentru rezolvarea ecuaiilor difereniale 8.1.1. Metoda dezvoltrii n serie Taylor pentru rezolvarea ecuaiilor difereniale Aceast metod asigur teoretic soluia oricrei ecuaii difereniale, dar are o aplicabilitate practic redus. Importana acestei metode const n faptul c servete drept criteriu de apreciere i comparare a metodelor utilizate frecvent n practic. Fie ecuaia diferenial de ordinul nti cu condiia iniial Cauchy:

( ) ( ) [ baxyxyyxfy ,, 00 == ] (8.1.1.1) Considerm c funcia f(x,y) este continu n x i este difereniabil de un numr de ori suficient de mare. Fie o divizare echidistant a intervalului [xo,b]:

Nmhxxhmxx mmm ,010 =+=+= + (8.1.1.2)

n baza dezvoltrii n serie Taylor n vecintatea lui xm, avem relaia:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )y x y x x x f x y x x f x y x xn

f x ymm

m mm

m mm

n

nm m= +

+

+ +

1 2

0

2

1 1

!,

!, ...

!,

cu eroarea de trunchiere: ( )( )

( )( ) ( ) [ ]x xn y xn

m

n

nm m+R x x

+

++=

+

1

1

111 !

, , (8.1.1.3)

de unde pentru x=xm+1 suficient de apropiat de xm, pentru h suficient de mic, obinem: ( ) ( ) ( ) ( ) ( ) (y y h f x y h f x y h

nf x ym m m m m m

nn

m m+= + + + +1

02

1

1 2!,

!, ...

!, )1 (8.1.1.4)

cu eroarea de trunchiere: ( ) ( )( ) ( )( ) [ ]R x hn f y x xn m

nn

m m+ +

+

+= +1 1

1

11 !, , ,

unde: y y ( ) (x y y xm m m+ += =1 1 ; .)m Caz particular n=1

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

y y y h f x y m

Rh

f y x x

m m m m

T m

0 10

2

21

1

1

2

; ,

, ,

+

+

= + =

= =

N

m

, (8.1.1.5)

Caz particular n=2

( ) ( ) ( ) ( )

( ) ( )( ) ( )

( ) ( ) ( ) ( )[ ]

y y yh

f x yh

f x y m N

Rh

f y x x

y y h f x yh

f x y f x y f x y

m m m m m m

T m m

m m m m x m m y m m m m

0 10

21

3

32

1

1

2

1 21

3

2

;!

,!

, ,

!, ,

, , ,

+

+

+

= + + =

= =

= + + +

,

(8.1.1.6)

Obs.: Odat cu considerarea unui numr mai mare de termeni "n", scade eroarea de trunchiere, dar crete numrul de derivate de calculat (uneori imposibil de calculat numeric), nefiind o metod folosit n practic.


8.1.2. Metoda Euler (metoda liniilor poligonale) Fie ecuaia diferenial de ordinul nti cu condiia iniial Cauchy:

( ) ( ) = =y f x y y x y, 0 0 (8.1.2.1) unde f(x,y) este definit ntr-un domeniu D din planul (xOy).

Fig. 8.1

Interpretarea geometric Din punct de vedere geometric metoda

Euler const n a nlocui curba y=y(x), soluia ecuaiei difereniale, prin linia poligonal construit cu ajutorul segmentelor de dreapt M M0 1,..., M Mi i+1,...., respectiv numai n punctul Mo (xo,yo) dreapta M M este tangent la curba y=y(x).

0 1

Se definete un cmp de direcii astfel nct n fiecare punct M(x,y)D se ia direcia cu panta ( ) ( ) = =tg f x y , .y x

Graficul soluiei ecuaiei difereniale (8.1.2.1) trece prin punctul M(xo,yo) i este tangent n orice punct al su la cmpul direciilor definite. Descrierea metodei Euler n metoda Euler se aproximeaz soluia printr-o linie poligonal n care fiecare segment este coliniar cu direcia cmpului definit de extremitatea sa stng. Considerm divizarea echidistant: x x i h i Ni = + =0 0, (8.1.2.2) pas"1" n punctul Mo(xo,yo) se calculeaz direcia cmpului definit n Mo i se scrie ecuaia dreptei determinat de Mo i aceast direcie: ( )( )f x y x x= + 0 0 0, .y y 0

]

Aceast funcie se propune ca aproximare a soluiei problemei (8.1.2.1) pe intervalul [x0,x1], de unde valoarea aproximativ a soluiei n x=x1 este:

( )y y h f x y1 0 0 0= + , (8.1.2.3) pas"i+1" Presupunem c n xi s-a calculat valoarea aproximativ yi, de unde pentru

soluia se aproximeaz prin: [x x xi i +, 1 ( )( )f x y x x x xi i i i i i= + = ++, 1y y hi+1

)ii respectiv n x=x obinem valoarea aproximativ:

(y y h f x yi i i+ = + 1 , (8.1.2.4) n metoda Euler, rezolvarea ecuaiei (8.1.2.1) se bazeaz pe utilizarea unei formule de calcul recursiv:

( )( )

y y x h x x

y y h f x y i Ni i

i i i i

0 0 1

1 0 1

= =

= + =

+

+ , , (8.1.2.5)

Relaia (8.1.2.5) se poate obine din relaia dezvoltrii n serie Taylor, considernd doar primii 2 termeni, n=1.

Obs.: Metoda Euler este cu pai separai i algoritmul metodei este de tip explicit. Pentru fiecare pas n metoda Euler, eroarea de trunchiere este:

( ) [ ] T i ih

y x x= +2

12

2, , T h (8.1.2.6)

obinut pe baza relaia erorii de trunchiere a metodei Taylor la n=1. Metoda Euler are o precizie de calcul redus. Pentru obinerea unei precizii mai bune de calcul se recomand folosirea unei divizri

mai fine (h mai mic) sau utilizarea metodei de tip predictor-corector.


Algoritm Start Introducere date: x , funcia f(x,y) y h N0 0, , , Cicleaz i=1,N { x xi i= +1 h Calculul funciei f x( )yi i 1 1, ( )y y h f x yi i i i= + 1 1, 1 } Scrie: ( )x y i Ni i, , ,= 0 Stop 8.1.3. Metoda Euler perfecionat (predictor-corector) Interpretarea geometric

Fig. 8.2

Fie ecuaia diferenial de ordinul nti cu condiia iniial:

( ) ( ) ( ) =y x f x y y x y, 0 = 0 (8.1.3.1) n metoda Euler perfecionat se folosete pentru pasul urmtor media pantelor din punctele: ( ) ( ) ( )( )x y x h y h y x ym m m m m m m, ; , ,+ + + +1 10 Pentru a gsi geometric punctul ( )x h y h ym m+ + , m

)m

)1

folosim metoda obinuit Euler, acesta gsindu-se pe dreapta (L1) din fig. 8.2.

( ) ( ) ( ) ( ) (L y x y y x x y y h f x ym m mx x

m m m

m

11

10: ,= + = +

=

+

+

-formul predictoare (8.1.3.2) Calculm n acest punct panta curbei soluiei ecuaiei (8.1.3.2) care trece prin punctul

, evalund n acest punct funcia f(x,y) se obine dreapta (L2) din fig. 8.2. ( )( x ym m+ +1 0, Apoi construim dreapta (L) a crei pant este egal cu media aritmetic a pantelor dreptelor (L1) i (L2). n final trasm prin punctul (xm,ym) o dreapt L paralel cu L i notm cu ( )punctul de intersecie al dreptei L cu verticala x=x

( )x ym m+ +1 11,

m+1.

Panta dreptei L (i a lui L) rezult: ( ) ( )([ ]+ + +12 1 10f x y f x ym m m m, , ) = , de unde:

( ) ( ) ( ) ( ) ( ) ( )( )[ ]L y y x x y y h y y h f x y f x ym m m m m m m m m m: ,= + = + = + ++ + 11 11 1 102 ,+ + x -formul corectoare (8.1.3.3) x hm m+ = +1

Descrierea metodei Euler predictor-corector Pe baza relaiilor (8.1.3.2) i (8.1.3.3) obinem algoritmul corector: iter"0" y y (predictor) ( ) ( )h f x y x xm m m m m m+ += + = +10 1, h


iter"1" ( ) ( ) ( )([ ]h f x y f x ym m m m m m+ += + +11 1 102 , , )+y y (corector) ..... (8.1.3.4)

iter"i" ( ) ( ) ( )([ ]h f x y f x ymi m m m m i++ + += + +11 1 12 , , )my y (corector) Criteriul de oprire: ( ) ( ) ( )y y dat y ym

imi

m mi

++

+ + > 11

1 10 +1 Algoritmul (8.1.3.4) se va utiliza pentru fiecare punct (xm,ym),m>0. Obs.: Se demonstreaz c algoritmul (8.1.3.4) este sigur convergent dac funcia f(x,y) satisface condiia lui Lipschitz n "y" i este continu n "x". Eroarea de trunchiere a metodei Euler mbuntit (predictor-corector) satisface relaia de proporionalitate: 3 (8.1.3.5) T h Algoritm Start Introducere date: ( )x y h N functia f x y iter0 0, , , , , , ,max Cicleaz i=1,N { x xi i= +1 h

i-1 i-1

1

Calculul funciei: f(x ,y ) ( )y y h f x yi i i i= + 1 1, iter=1 Repet { ts t yi= Calculul funciei: f(xi,yi)

( ) ([ ]+ 12 1 1f x y f x yi i i i, , ) = y y hi i= + 1 ts t tst yi= iter=iter+1 } pn cnd: ( ) ( )tst sau iter iter< > max } Scrie: ( )x y i Ni i, , ,= 0 Stop 8.1.4. Metoda Runge-Kutta Metodele Runge-Kutta au 3 proprieti:

1) sunt metode directe, au nevoie doar de informaiile existente la punctul precedent (xk,yk) pentru a gsi (xk+1,yk+1);

2) sunt identice cu seriile Taylor pn la termenii hp, unde p este diferit pentru metode distincte (ordinul metodei);

3) nu necesit evaluarea nici unei derivate a funciei f(x,y), ci numai a valorii acesteia ntr-un punct.

Obs. Ca urmare a proprietii (3), metodele Runge-Kutta sunt cele mai folosite n practic. Fie ecuaia diferenial de ordinul nti cu condiia iniial Cauchy:

( ) ( )( ) ( ) = =y x f x y x y x y, 0 0 (8.1.4.1)


Fie o diviziune echidistant a intervalului [xo,b]: xm=xo+mh, m=0,N (8.1.4.2) Obs.: Pe baza metodei seriei Taylor de ordinul "p" a funciei y(x), pentru diviziunea (8.1.4.2) avem:

( ) ( )( ) ( )( )

( )( ) ( )( ) ( ) ( )( ) ( )( ) ( )( ) ( )yxfyxfyxfyxfxxyf

phR

Nmyxfphyxfhyxfhyy

py

px

pmm

pp

p

mmp

p

mmmmmm

,,,,,,,,!1

1,0,!

...,!2

,!1

111

1

1

112

1

+=+

=

=++++=

+

+

+

+

(8.1.4.3) Relaia (8.1.4.3) prin metoda Runge-Kutta va fi nlocuit cu o formul de recuren de tipul:

( )( )( )

( )

y daty y k k k

k h f x y

k h f x h y k

k h f x h y k k

k h f x h y k k k

m m p p

m m

m m

m m

p m p m p p p p

0

1 0 0 1 1

0

1 1 10 0

2 2 20 0 21 1

0 0 1 1 1 1

+

= + + + +

=

= + +

= + + +

= + + + + +

...

,

,

,.....

, ... , p

(8.1.4.4)

n funcie de ordinul "p" vom evidenia 4 dintre cele mai importante variante ale

metodelor Runge-Kutta. 8.1.4.1. Formula Runge-Kutta de ordinul I (p=1)

Din dezvoltarea n serie Taylor (8.1.4.3) avem:

( )y y h f x y p y datm m m m+ = + =1 1 1! , , 0 (8.1.4.5)

Din forma general a metodei Runge-Kutta (8.1.4.4) avem:

( ) ( )y y kk h f x y y y h f x y y d

m m

m mm m m m

++

= +=

= + 1 0 0

01 0 0

, , ; at (8.1.4.6)

Punem condiia de identitate a relaiilor (8.1.4.5),(8.1.4.6) i din identificarea

termenilor rezult: 0 1= , de unde formula Runge-Kutta de ordinul I este:

( )y y k y d

k h f x ym m

m m

+ = +

= 1 0 0

0

;

,

at (8.1.4.7)

Obs.: Formula Runge-Kutta de ordinul I este identic cu formula liniilor poligonale a

lui Euler. Eroarea de trunchiere este:

( )( )( )

Th

f y hI

= 2

2

2!, T (8.1.4.8)


Algoritm Start Introducere date: y x h ( )N f x y0 0, , , , , Cicleaz m=1,N xm=xm-1+h Cicleaz m=0,N-1 {

( )k h f x y

y y km m

m m

0

1 0

=

= ++

,

} Scrie: ( )x y m Nm m, , ,= 0 Stop 8.1.4.2. Formula Runge-Kutta de ordinul II (p=2)


( ) ( ) ( )y y h f x y h f x y p y dam m m m m m+ = + + =12

101 2

2!

,!

, , t (8.1.4.9)


( )( )01

0

011001

,,

;

kyhxfhkyxfhk

datkkkyy

mm

mm

mm

++==

++=+ (8.1.4.10)

Dezvoltm n serie Taylor termenul k1, din care reinem doar primii 2 termeni:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )mmymmmxmmmmmmm yxfyyyxfxxyxfyxdfyxfyxf ,,,,!11,, ++=+

( ) ( ) ( ) ( ) ( )mmmmymmxmmmm yxfyxfhyxfhyxfkyhxf ,,,,, 0 ++=++ (8.1.4.11)

Din relaiile: (8.1.4.10),(8.1.4.11) avem: ( ) ( ) ( ) ( ) ( )[ ]y y h f x y h f x y f x y f x ym m m m x m m y m m m m+ = + + + + 1 0 1 1 2 , , , ,

iar din relaia (8.1.4.9) avem:

( ) ( ) ( ) ([ ]y y h f x y h f x y f x y f x ym m m m x m m y m m m m+ = + + + 12

2, , , ),

i din identificarea celor dou relaii avem:

21211

1

1

10

==

=+

, considerm 0 11 2 1 2 1= = = = (8.1.4.12)

Din relaiile: (8.1.4.10),(8.1.4.12), formula Runge-Kutta de ordinul II este: ( )

( )( )01

0

0101

,,

;2

kyhxfhkyxfhk

datykkyy

mm

mm

mm

++==

++=+ (8.1.4.13)

i este echivalent cu metoda Euler perfecionat. Obs.: Se demonstreaz c eroarea de trunchiere satisface relaia:

( )T II h

3 (8.1.4.14)


Algoritm Start Introducere date: y x h ( )N f x y0 0, , , , , Cicleaz m=1,N xm=xm-1+h Cicleaz m=0,N-1

( )( )

( ) 2,

,

101

01

0

kkyykyhxfhk

yxfhk

mm

mm

mm

++=++=

=

+

Scrie: ( )x y m Nm m, , ,= 0 Stop 8.1.4.3. Formula Runge-Kutta de ordinul III (p=3)


( ) ( )( ) ( )( ) datypyxfhyxfhyxfhyy mmmmmmmm 023

12

1 ,3,!3,

!2,

!1=+++=+ (8.1.4.15)


( )( )( )

y y k k k k dat

k h f x y

k h f x h y k

k h f x h y k k

m m

m m

m m

m m

+ = + + +

=

= + +

= + + +

1 0 0 1 1 2 2 0

0

1 1 10 0

2 2 20 0 21 1

;

,

,

,

(8.1.4.16)

Vom dezvolta funcia f(x,y) n serie Taylor n punctul (xm,ym), reinnd termenii de

ordinul I i II: ( ) ( ) ( ) ( ) ( ) ( )

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

f x y f x y x x f x y y y f x y

x xf x y x x y y f x y

y yf x y

m m m x m m m y m m

mxx m m m m xy m m

myy m m

, , , ,

, ,

= + + +

+

+ +

2 2

2 2,

(8.1.4.17)

Facem notaiile: ( ) ( ) ( )

( ) ( ) (f f x y f f x y f f x y

f f x y f f x y f f x y

m m x x m m y y m m

xx xx m m xy xy m m yy yy m m

= = =

= = =

, ; , ; , ;

, ; , ; , )

f

(8.1.4.18)

Pe baza relaiilor (8.1.4.17),(8.1.4.18) obinem:

k h0 =

( ) ( )k h f h f f f h f f f fx y xx xy yy12

1 10

3

12

1 10 102 2

22 2

63 6 3= + + + + + f (8.1.4.19)

( )

( ) ( )( )k h f

hf f f f f

hf f f f f f f

x y y

xx xy yy x y y

2

2

2 20 21

3

22

2 20 21 20 21

2 21 21 10 21

2

22 2 2

63 6 3 6 6

= + + + +

+ + + + + + +

f f

Din relaiile (8.1.4.16),(8.1.4.19) pentru formula Runge-Kutta de ordinul III avem:


( ) ( ) +

++++

+++=+22

1010121

3

101

2

101 363622

2fffffhfffhfhfhyy yyxyxxyxmm

( )

+++++ fffffhfh yyx 212022

2 2222 (8.1.4.20)

( ) ( )( )

+++++++ fffffffffh yyxyyxyxx2

211021122

21202120222

3

663636

i respectiv relaia seriei Taylor (8.1.4.15), descompus este:

( ) ( )fffffffffhfffhfhyy yyxyyxyxxyxmm ++++++++=+ 2232

1 262 (8.1.4.21)

Se identific termenii din relaiile (8.1.4.20), (8.1.4.21) i rezult urmtorul sistem de ecuaii:

( )

( )( )

6161

31

3131

2121

1

21102

2112

221202

2101

2120221011

222

211

21202101

2211

210

==

=++

=++=+

=++=+=++

21220

1221

110

222

211

2211

201

61

3121

1

===

=+=+

=

(8.1.4.22)

unde s-au considerat ca fiind variabilele principale. 0 , 2Varianta 1 Se consider 0 2 1 6= = i din sistemul de ecuaii (8.1.4.22) rezult urmtoarea formul Runge-Kutta de ordinul III:

( )

( )

( )102

01

0

02101

2,2

,2

,

;461

kkyhxfhk

kyhxfhk

yxfhk

datykkkyy

mm

mm

mm

mm

++=

++=

=

+++=+

(8.1.4.23)

Varianta 2 Se consider 0 21 4 3 4= =; i din sistemul de ecuaii (8.1.4.22) rezult urmtoarea formul Runge-Kutta de ordinul III:

( )

( )

++=

++=

=

++=+

32,

32

3,

3

,

;341

12

01

0

0201

kyhxfhk

kyhxfhk

yxfhk

datykkyy

mm

mm

mm

mm

(8.1.4.24)

Obs.: Se demonstreaz c eroarea de trunchiere a metodei Runge-Kutta de ordinul III satisface relaia de proporionalitate:

( )T III h

4 (8.1.4.25)


Algoritm Start Introducere date: y x h ( )N f x y0 0, , , , , Cicleaz m=1,N xm=xm-1+h Cicleaz m=0,N-1 {

( )

( )

( )

k h f x y

k h f xh

yk

k h f x h y k k

y y k k k

m m

m m

m m

m m

0

10

2 0

1 0 1

2 2

216

4

=

= + +

= + +

= + + ++

,

,

, 1

2

} Scrie: ( )x y m Nm m, , ,= 0 Stop 8.1.4.4. Formula Runge-Kutta de ordinul IV (p=4)

Din dezvoltarea n serie Taylor (8.1.4.3) avem (p=4): yo dat ;

( ) ( ) ( ) ( ) ( ) ( ) (y y h f x y h f x y h f x y h f x ym m m m m m m m m m+ = + + + +12

13

24

3

1 2 3 4!,

!,

!,

!, ) (8.1.4.26)


( )( )( )( )

y y k k k k k dat

k h f x y

k h f x h y k

k h f x h y k k

k h f x h y k k k

m m

m m

m m

m m

m m

+ = + + + +

=

= + +

= + + +

= + + + +

1 0 0 1 1 2 2 3 3 0

0

1 1 10 0

2 2 20 0 21 1

3 3 30 0 31 1 32 2

;

,

,

,

,

(8.1.4.27)

Facem notaiile:

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

( ) ( ) ( ) ( )



f f x y f f x y f f x y f f x y

m m x x m m y y m m

xx xx m m xy xy m m yy yy m m

xxx xxx m m xxy xxy m m xyy xyy m m yyy yyy m m

= = =

= = =

= = = =

, ; , ; , ;

, ; , ; ,

, ; , ; , ; ,

(8.1.4.28)

Descompus, relaia (8.1.4.26) cu notaiile (8.1.4.28), are expresia:

( ) ( )

()

y y h fh

f f fh

f f f f f f f f f

hf f f f f f f f f f f f

f f f f f f f f f f f f

m m x y xx xy yy x y y

xxx xxy xy x xyy yy x xx y

xy y y x yyy yy y y

+ = + + + + + + + + +

+ + + + + + +

+ + + + +

1

2 32 2

42

2 2 2 3

2 62

243 3 3 3

5 4

(8.1.4.29)


Dezvoltm n serie Taylor funcia f(x,y) i reinem primii 4 termeni:

( ) ( ) ( ) ( ) (f x y f x y df x y d f x y d f x ym m m m m m m m, , ! , ! , ! ,= + + +11

12

13

2 3 ) (8.1.4.30)

( ) ( ) ( )[ ] ( ) ( )( ) ( )[ ] = + + + + + f x y f x x f y y f x x f x x y y f y y fm x m y m xx m m xy m yy, 12 22 2

( ) ( ) ( ) ( )( ) ( )[ ]+ + + 16 3 33 2 2x x f x x y y f x x y y f y y fm xxx m m xxy m m xyy m yyy+

3

Folosind relaia (8.1.4.30), se calculeaz termenii k1,k2,k3 din relaia (8.1.4.27) i se compar cu relaia (8.1.4.29). Din identificarea termenilor asemenea rezult sistemul ecuaiilor coeficienilor necunoscui (analog metodei Runge-Kutta de ordinul III). Varianta 1 (formula treimii pentru metoda Runge-Kutta de ordinul IV)

( )( )

( )

y y k k k k y d

k h f x y

k h f xh

yk

k h f xh

yk

k h f x h y k

m m

m m

m m

m m

m m

+ = + + + +

=

= + +

= + +

= + +

1 0 1 2 3 0

0

10

21

3 2

16

2 2

2 2

2 2

;

,

,

,

,

at

(8.1.4.31)

Varianta 2 (formula optimii pentru metoda Runge-Kutta de ordinul IV)

( )( )

( )

y y k k k k y d

k h f x y

k h f xh

yk

k h f xh

y kk

k h f x h y k k k

m m

m m

m m

m m

m m

+ = + + + +

=

= + +

= + +

= + + +

1 0 1 2 3 0

0

10

2 10

3 2 1 0

18

3 3

3 323 3

;

,

,

,

,

at

(8.1.4.32)

Obs.: Se demonstreaz c eroarea de trunchiere a metodei Runge-Kutta de ordinul IV

satisface relaia de proporionalitate: ( )

T III h5 (8.1.4.33)

Algoritm Start Introducere date: y x h ( )N f x y0 0, , , , , Cicleaz m=1,N xm=xm-1+h Cicleaz m=0,N-1 ( )mm yxfhk ,0 =

++=

2,

20

1kyhxfhk mm


++=

2,

21

2kyhxfhk mm

( )23 , kyhxfhk mm ++=

( )32101 2261 kkkkyy mm ++++=+

Scrie: ( )x y m Nm m, , ,= 0 Stop 8.2. Metode numerice indirecte pentru rezolvarea ecuaiilor difereniale 8.2.1. Metoda Adams-Bashforth Descrierea metodei

Fie ecuaia diferenial de ordinul nti cu condiia iniial Cauchy: ( ) ( ) ( ) = =y x f x y y x y, 0 0 (8.2.1.1)

Fie o diviziune echidistant a intervalului [xo,xN] de pas h>0:

x x k h k Nk = + =0 0, (8.2.1.2)

Presupunem c printr-o anumit metod numeric direct s-a construit tabelul de valori care aproximeaz soluia y(x) a problemei (8.2.1.1) n nodurile diviziunii:

x0 x1 xn (8.2.1.3)y0 y1 yn

Fie mn , mN . Notm fk=f(xk,yk), k=0,n i cu pm(x) polinomul de interpolare tip

Lagrange generat de tabelul de valori care aproximeaz funcia f(x,y):

xn-m xn-m+1 xn (8.2.1.4)fn-m fn-m+1 fn

Formula exact de calcul a soluiei problemei (8.2.1.1):

( ) ( ) ( )( )y x y x f x y x dxn nx

x

n

n

+ =+

11

, (8.2.1.5)

se nlocuiete cu formula aproximativ:

( )y y p x dn n mx

x

n

n

+ =+

11

x

m

formula Adams-Bashforth (8.2.1.6)

Polinomul p (x) de tip Lagrange are relaia:

( )( ) ( )( ) ( )

( ) ( )( ) ( )p xx x x x x x x x

x x x x x x x xfm

n m i i n

i n m i i i i i ni n m

n

=

+

+=

... ...... ...

1 1

1 1i (8.2.1.7)

de unde obinem:

( )p x dx A fmx

x

ii n m

n

n

n+

= =

1

i ; ( ) ( )( ) ( )

( ) ( )( ) ( )A (8.2.1.8) x x x x x x x x

x x x x x x x xdxi

n m i i n

i n m i i i i i nx

x

n

n

=

+

+

+

... ...

... ...1 1

1 1

1

Calculul coeficienilor Ai, i=n-m,...,n se face considernd schimbarea de variabil:

( ) [x x t x t x h t tn n n= + = + +1 01 , ]1 de unde rezult:


( )

x x hx x h

x x n i h

i i

i i

i n

= =

=

+

+

1

2 2.....

,

( )

x x hx x h

x x i n m

i i

i i

i n m

= =

= +

1

2 2.....

h

, ( )

( )

x x h tx x h t

x x h t m

n

n

n m

= = +

= +

1 1.....

(8.2.1.9)

Din relaiile (8.2.1.8),(8.2.1.9) obinem: ( )

( ) ( )( ) ( )

Ah

i n m n it t t m

t n idt i n m ni

n i

=

+ + +

+ =

1 10

1

! !.....

,..., (8.2.1.10)

Observaii: - Pe baza relaiei (8.2.1.10) rezult c coeficienii Ai nu depind de ecuaia diferenial

iniial, deci se pot calcula i folosi apoi pentru orice problem (8.2.1.1). n baza relaiei (8.2.1.8) formula Adams-Bashforth devine:

i (8.2.1.11) y y An n ii n m

n

+=

= + 1 f- Se demonstreaz c eroarea de trunchiere a metodei Adams-Bashforth satisface

relaia de proporionalitate: (8.2.1.12) Tmh +2

- Metoda Adams-Bashforth datorit modului cum se calculeaz aproximaia la un anumit pas "n+1", n funcie de paii anteriori, poart denumirea i de "metoda cu pai legai". Cazuri particulare

(m y y h f fn n n n= = + +1 2 31 )1 (8.2.1.13) Pentru a determina condiiile iniiale pentru (8.2.1.13), folosim condiia Cauchy (xo,yo)

i cu metoda Runge-Kutta de ordinul II calculm (x1,y1).

(m y y h f f fn n n n n= = + ++2 12 23 16 51 ) 1 2 (8.2.1.14) Pentru a determina condiiile iniiale pentru (8.2.1.14), folosim condiia Cauchy (xo,yo)

i cu metoda Runge-Kutta de ordinul III calculm (x1,y1), (x2,y2).

(m y y h f f f fn n n n n n= = + + + 3 24 55 59 37 91 1 ) 2 3 (8.2.1.15) Pentru a determina condiiile iniiale pentru (8.2.1.15), folosim condiia Cauchy (xo,yo)

i cu metoda Runge-Kutta de ordinul IV calculm (x1,y1), (x2,y2), (x3,y3). Obs.: Am grupat metoda Runge-Kutta cu metoda Adams-Bashforth de un anumit

ordin, astfel nct s avem acelai ordin de mrime pentru eroarea de trunchiere. Algoritm Start Introducere date: y x N ( )h f x y0 0, , , , , Cicleaz i=1,N x xi i= +1 h

) ( )(

( )

k h f x y

k h f x h y k

y y k k

0 0 0

1 0 0

1 0 0 1

12

=

= + +

= + +

,

, 0 {Runge-Kutta II}

Cicleaz i=1,N-1

( ) ([y y h f x y f x yi i i i i i+ = + 1 2 3 , , )] 1 1 {Adams-Bashforth m=1}


Scrie: ( )x y i Ni i, , ,= 0 Stop 8.2.2. Metoda Adams-Moulton Descrierea metodei

Fie ecuaia diferenial de ordinul nti cu condiia iniial Cauchy: ( ) ( ) ( ) = =y x f x y y x y, 0 0 (8.2.2.1)

Fie o diviziune echidistant a intervalului [xo,xN] de pas h>0: x x k h k Nk = + =0 0, (8.2.2.2)

Presupunem c printr-o anumit metod numeric am construit tabelul de valori care aproximeaz soluia y(x) a problemei (8.2.2.1) n nodurile diviziunii: x x x x xy y y y y

n n

n n

0 1 2 1

0 1 2 1

.....

.....+

+

(8.2.2.3)

Fie mn , mN . Notm fk=f(xk,yk), k=0,n+1 i cu qm+1(x) polinomul de interpolare tip

Lagrange generat de tabelul de valori care aproximeaz funcia f(x,y(x)): x x x xf f f f

n m n m n n

n m n m n n

+

+

1

1 1

..........

+

+

1

x

(8.2.2.4)

Formula exact de calcul a soluiei problemei (8.2.2.1):

( ) ( ) ( )( )y x y x f x y x dxn nx

x

n

n

+ =+

11

, (8.2.2.5)

se nlocuiete cu formula aproximativ:

( )y y q x dn n mx

x

n

n

+ + =+

1 11

formula Adams-Moulton (8.2.2.6)

Formula (8.2.2.6) a lui Adams-Moulton genereaz o ecuaie cu necunoscuta yn+1, care

se rezolv printr-o metod iterativ.

Avem ecuaia cu necunoscuta yn+1 de tipul ( )yn+ =1y F i care se rezolv iterativ conform: y F (8.2.2.7)

n+1

N[ ] [ ]( )y knk nk++ += 11 1 ,

Observaii: - Pentru a determina aproximaia iniial [ ]yn+1 se folosete o metod direct de calcul.

Metoda Adams-Moulton are rolul de a mbunti precizia de calcul asigurat de alte metode.

0

- Pentru o funcie f(x,y) ce satisface proprietatea Lipschitz n raport cu "y", convergena metodei Adams-Moulton (8.2.2.7) este asigurat.

Polinomul de interpolare Lagrange, pe baza tabelului (8.2.2.4) are expresia:

( )( ) ( )( ) ( )

( ) ( )( ) ( )q xx x x x x x x x

x x x x x x x xfm

n m i i n

i n m i i i i i ni n m

n

i+ + +

+=

+

=

1

1 1 1

1 1 1

1 ... ...... ... +

i

(8.2.2.8)

de unde obinem:

( )q x dx B fmx

x

ii n m

n

n

n

+=

++

= 111

; ( ) ( )( ) ( )

( ) ( )( ) ( )B (8.2.2.9) x x x x x x x x

x x x x x x x xdxi

n m i i n

i n m i i i i i nx

x

n

n

=

+ +

+

+

... ...

... ...1 1 1

1 1 1

1

+


Analog metodei Adams-Bashforth vom considera schimbarea de variabil:

( ) [x x t x t x h t tn n n= + = + +1 01 , ]1 de unde rezult:

( )

x x hx x h

x x n i

i i

i i

i n

= =

= +

+

+

+

1

2

1

2

1.....

h

,

( )

x x hx x h

x x i n m h

i i

i i

i n m

= =

= +

1

2 2.....

,

( )

( )

x x h tx x h t

x x h t m

n

n

n m

= =

= +

+

1 1

..... (8.2.2.10)

Din relaiile (8.2.2.9),(8.2.2.10) obinem:

( )( ) ( )

( ) ( ) ( )B

hi n m n i

t t t t mt n i

dt i n m nin i

=

+ + + +

+ = +

+

11

1 11

1

0

1

! !.....

,..., (8.2.2.11)

- Coeficienii Bi nu depind de ecuaia diferenial iniial, se pot tabela i folosi apoi

pentru orice problem (8.2.2.1). Din relaiile (8.2.2.6),(8.2.2.9) formula Adams-Moulton devine:

[ ] [ ]( )y y B f y y B f B f x yn n i ii n m

n

nk

n i ii n m

n

n n nk

+=

+

++

= + += + = + + 1

1

11

1 1, +1 k>0 (8.2.2.12)

cu criteriul de oprire: [ ] [ ]y y dnk

nk

++

+ < >11

1 0 at (8.2.2.13) - Se demonstreaz c eroarea de trunchiere a metodei Adams-Moulton satisface relaia

de proporionalitate: Tmh +3 (8.2.2.14)

- Metoda Adams-Moulton ca i metoda Adams-Bashforth este tot "o metod cu pai legai". Cazuri particulare

[ ] ( ) [ ]([m y y h f x y f x ynk n n n n k= = + +++ + +0 211

1 1, , )]n n0 (8.2.2.15) Expresia de mai sus reprezint relaia de corecie de la metoda Euler mbuntit

(predictor-corector). Folosind condiia Cauchy (xo,yo) i cu metoda Runge-Kutta de ordinul II obinem

aproximrile iniiale ( ) . [ ]x y i Ni i, , ,0 1=[ ] [ ]([m y y h f f f x ynk n n n n k= = + + +++ +1 12 8 51

11 1 , )]n+1

o o

n1 (8.2.2.16)

Pentru a determina condiiile iniiale pentru (8.2.2.16), folosim condiia Cauchy (x ,y ) i cu metoda Runge-Kutta de ordinul III obinem aproximrile iniiale

. [ ]( )x y i Ni i, , ,0 1=Pentru m=1 , soluia ( )[ ] ( )x y x y1 10 1 1, , nu se mai prelucreaz suplimentar.

[ ] [ ]([m y y h f f f f x ynk n n n n n nk= = + + +++ + +2 24 5 19 911

2 1 1 , )]1 n2 (8.2.2.17) Pentru a determina condiiile iniiale pentru (8.2.2.16), folosim condiia Cauchy (xo,yo)

i cu metoda Runge-Kutta de ordinul IV obinem aproximrile iniiale [ ]( )x y i Ni i, , ,0 1= . Pentru m=2 , soluia [ ]( ) ( )x y x y1 10 1 1, , ; [ ]( ) ( )x y x y2 20 2 2, , nu se mai prelucreaz

suplimentar. Algoritm Start Introducere date: y x N ,iter(h f x y0 0, , , , , ) max,


Cicleaz i=1,N x xi i= +1 h Cicleaz i=0,N-1 {

( )

( )

( )

k h f x y

k h f xh

yk

k h f x h y k k

y y k k k

i i

i i

i i

i i

0

10

2 0

1 0 1

2 2

216

4

=

= + +

= + +

= + + ++

,

,

, 1

2

{Runge-Kutta III}

} Cicleaz i=1,N-1 { iter=1 Repet { tst=yi+1

( ) ( ) ([ ]h f x y f x y f x yi i i i i i i i+ = + + +1 1 112 8 5, , , )+ +1 1y y ts {Metoda Adams-Moulton m=1} t tst yi= +1 iter=iter+1 } pn cnd: ( ) ( )tst sau iter iter< > max } Scrie: ( )x y i Ni i, , ,= 0 Stop

8.3. Sisteme de ecuaii difereniale 8.3.1. Sisteme de ecuaii difereniale de ordinul nti. Ecuaii difereniale de ordin superior

Fie sistemul de ecuaii difereniale de ordinul nti:

( )

( )

( )

dydx

f x y y ydydx

f x y y y

dydx

f x y y y

m

m

mm m

11 1 2

22 1 2

1 2

=

=

=

, , ,...,

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

, , ,...,

(8.3.1.1)

cu condiiile iniiale Cauchy: [ ] ( ) ( ) ( )x a b y x y y x y y x ym m = = =, ; ; ; ...1 0 1 2 0 2 00 0 0 (8.3.1.2)

Definiie: Funcia f(x,y ,y ,...,y ) satisface condiia Lipschitz n variabilele y ,y ,...,y

pe mulimea 1 2 m

( ) [ ]1 2 m

{ }D x y y y x a b y R j mm j= , , ,..., , , , ,1 2 1= dac exist o constant L>0 astfel nct:


( ) ( )f x y y y f x z z z L y zm mj

m

, , ,..., , , ,...,1 2 1 21

=

(x y y y x z z zm m, , ,..., ; , , ,...,1 2 1 2 j j

) (8.3.1.3)

pentru orice ( ) . D

Obs.: Se demonstreaz c dac funciile din sistemul de ecuaii difereniale (8.3.1.1) sunt continue i satisfac condiia Lipschitz (8.3.1.3), atunci sistemul de ecuaii difereniale (8.3.1.1), (8.3.1.2) are soluie i este unic.

Fie ecuaia diferenial de ordinul "m" de forma: ( ) ( )( ) [ ]y f x y y y x a bm m= , , ,..., ,(1) 1 (8.3.1.4)

cu condiiile iniiale: ( ) ( ) ( ) ( ) ( ) ( ) ( )y y x y y x ym0 0 1 0 01 1 0 0 1= = ; .....y x (8.3.1.5) m=

Pentru rezolvarea ecuaiei difereniale de ordinul "m", aceasta va fi transformat ntr-un sistem de ecuaii difereniale de ordinul nti echivalent.

Notm: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )y x y x y x y x y xm m1 2 1=y x 1= = ; ..... (8.3.1.6) i obinem sistemul de ecuaii difereniale de ordinul nti echivalent:

( )

( )

( )( ) ( )

dydx

dydx

y

dydx

dydx

y

dydx

dydx

y

dydx

dydx

y f x y y y

mm

m

mm

mm

12

21

3

12

1

1 2

= =

= =

= =

= = =

.....

, , ,...,

(8.3.1.7)

cu condiiile iniiale: ( ) ( ) ( ) ( ) (y y x y y x ym m1 0 0 2 0 01 0 0 1= = ; ..... )y x (8.3.1.8) = Din relaiile (8.3.1.4)-(8.3.1.8) rezult c problema ecuaiei difereniale de ordin superior se reduce la a rezolva tot un sistem de ecuaii difereniale de ordinul nti echivalent de forma (8.3.1.1),(8.3.1.2). Pentru rezolvarea sistemului de ecuaii (8.3.1.1),(8.3.1.2) se vor aplica generalizri ale metodelor numerice folosite n cazul ecuaiilor difereniale de ordinul nti (vezi capitolele 8.1 i 8.2). Spre exemplificare vom prezenta generalizarea metodei Runge-Kutta IV de la ecuaii difereniale de ordinul nti la sisteme de ecuaii difereniale de ordinul nti. Fie o divizare echidistant a intervalului [a,b]:

x a hb a

Nx x i h ii0 0 0= =

= + =, , , ,

T

N

}

)}

(8.3.1.9)

Notm cu { } vectorul necunoscutelor, astfel nct sistemul de ecuaii (8.3.1.1),(8.3.1.2) poate fi adus la forma matriceal:

{Y y y ym= 1 2, ,...,

{ } { }({ =Y F x Y, cu condiia iniial: ( ){ } { }Y x Y0 0= (8.3.1.10) unde: { }{ }( ) ( ) ( ) ( ){ }F x Y f x y y y f x y y y f x y y ym m m

T

, , , ,..., ; , , ,..., ;.....; , , ,...,= 1 1 2 2 1 2 1 2 m

{ } { }Y y y ymT

0 1 20 0 0= , ,.....,

Relaia de calcul a vectorului necunoscutelor { }Yi+1 n punctul divizrii xi+1=xi+h prin generalizarea metodei Runge-Kutta IV este:


( ){ } { }{ } { } { } { } { } { }( )Y x Y

Y Y K K K K x x h i Ni i i i

0 0

1 1 2 3 4 1

16

2 2 0

=

= + + + + = + = + +; ; 1, (8.3.1.11)

unde: { } { }( ){ }{ } { } { }

{ } { } { }

{ } { } { }( ){ }

K h F x Y

K h F xh

Y K

K h F xh

Y K

K h F x h Y K

i i

i i

i i

i i

1

2 1

3 2

4 3

212

212

=

= + +

= + +

= + +

,

,

,

,

Sub form algebric relaiile (8.3.1.11) devin:

( )y x y j mj j0 0 1= = , ( )y y k k k k j m x x h i Nj i j i j j j j i i, , ; , ,+ += + + + + = = + = 1 1 2 3 4 116 2 2 1 0 1 (8.3.1.12)

unde: ( )k h f x y y y jj j i i i m i1 1, 2 1= =, , ,..., ,, , m

k h f xh

yk

yk

yk

j mj j i i i m im

2 1,11

212 1

2 2 2 21= + + + +

=, , ,..., ,, ,

k h f xh

yk

yk

yk

j mj j i i i m im

3 1,21

222 2

2 2 2 21= + + + +

=, , ,..., ,, ,

( )k h f x h y k y k y k jj j i i i m i m4 1, 31 2 32 3 1= + + + + =, , ,..., ,, , m Algoritm Start Introducere date: ( ) ( )( )y j m h N f x y y y j mj j m, , , , , , , , , ,..., , ,0 1 20 1 1= =m x Cicleaz i=1,N x xi i+ = +1 h

m i

Cicleaz i=0,N-1 { Cicleaz j=1,m k h ( )f x y y yj j i i i1 1, 2= , , ,...,, , Cicleaz j=1,m

k h f xh

yk

yk

yk

j j i i i m im

2 1,11

212 1

2 2 2= + + + +

, , ,...,, , 2

Cicleaz j=1,m

k h f xh

yk

yk

yk

j j i i i m im

3 1,21

222 2

2 2 2= + + + +

, , ,...,, , 2

m3

Cicleaz j=1,m k h ( )f x h y k y k y kj j i i i m i4 1, 31 2 32= + + + +, , ,...,, , Cicleaz j=1,m


( )k k k kj i j i j j j j, ,+ = + + + +1 1 2 316 2 2y y 4 } Scrie: ( ) x y y y i Ni i i m i, , ,..., , ,, ,1, 2 0= Stop 8.3.2. Metoda -Newmark pentru rezolvarea sistemelor de ecuaii difereniale de ordinul doi

Fie sistemul de ecuaii difereniale de ordinul doi de forma: [ ] ( ){ } [ ] ( ){ } [ ]{ } ( ){ } [ ]A Y B Y C Y F x x a b2 1+ + = , (8.3.2.1) cu condiiile iniiale: { }( ) { } ( ) ( ){ } ( ){ }Y x Y Y x Y0 0 1 0 01= =; (8.3.2.2)

unde:

{ } { }[ ] [ ] [ ]

( ){ } ( ) ( ) ( ){ }

Y y y y

A B C M

F x f x f x f x

m

T

mxm

m

T

=

=

1 2

1 2

, ,...,

, ,

, ,...,

Obs.: Problema (8.3.2.1),(8.3.2.2) are o aplicaie direct n studiul dinamicii sistemelor mecanice.

Fie divizarea echidistant: NihixxN

abhax i ,0;; 00 =+=

== (8.3.2.3)

Algoritmul -Newmark pentru rezolvarea problemei (8.3.2.1),(8.3.2.2) este urmtorul:

pas"0" Condiia iniial: { } ( ){ }Y0 0 01, ,x Y , ( ){ } [ ] ( ){ } [ ] ( ){ } [ ]{ }( )Y A F x B Y C Y0 2 1 0 01 0= unde: det [ ]A 1 0..... pas"i" Se determin: x Y { } ( ){ } ( ){ }Y Yi i i i, , ,1 2pas"i+1" (i=0,N-1) { } ( ){ }h F F xi i i i+ += + =1 1; x x +1

( ){ } ( ){ } ( ){ } ( ){ }[ ]

{ } { } ( ){ } ( ){ } ( ){ }[ ]

Y Yh

Y Y

Y Y h Yh

Y Y

i i i i

i i i i i

+ +

+ +

= + +

= + + +

11 1 2

12

11

22

12

2

4

(8.3.2.4)

[ ] ( ){ } [ ] ( ){ } [ ]{ } { }A Y B Y C Y Fi i i+ + ++ + =12 11 1 i+1 i nlocuind { } ( ){ }Y Yi i+ +1 1, 1 obinem sistemul de ecuaii liniare cu necunoscutele { } : ( )Yi+12

[ ] [ ] [ ] ( ){ } { } [ ] [ ] ( ){ } [ ] [ ]( ) ( ){ } [ ]{ }A h B h C Y F h B h C Y B h C Y C Yi i i i+ + i

= +

+ + +2 4 2 2

2

12

12 1

i care este echivalent cu forma restrns: [ ] { }{ } { }D Y Ei+ =12 , sistem algebric liniar rezolvabil prin metoda Gauss.

Obs.: De la pasul anterior pas"i" avem: [ ] ( ){ } [ ] ( ){ } [ ]{ } { }A Y B Y C Y Fi i i2 1+ + = i , de unde se poate obine la pas"i+1" i sistemul algebric liniar echivalent de forma:


[ ] [ ] [ ] ( ){ } { } { } [ ] [ ] [ ] ( ){ } [ ] ( ){A h B h C Y F F A h B h C Y h C Yi i i i+ + }i

= +

+ +2 4 2 4

2

12

1

22 1 (8.3.2.5)

Algoritm Start Introducere date: { } ( ){ } [ ] [ ] [ ] ( ){ }h Y Y A B C F x0 0 01, , , , , , , ,x N Cicleaz i=1,N x xi i+ = +1 h

Calculeaz: { }( ) [ ] ( ){ } [ ] ( ){ } [ ]{ }( )Y A F x B Y C Y0 2 1 0 01 0= 2

Calculeaz: [ ] [ ] [ ] [ ]D Ah

Bh

C= + +2 4

Cicleaz i=0,N-1 { { } ( ){ }F F xi i+ +=1 1 { } { } [ ] [ ] ( ){ } [ ] [ ]( ) ( ){ } [ ]{ }E F h B h C Y B h C Y C Yi i i= +

+ +1

2 1

2 2 i

Rezolvare prin metoda Gauss: [ ] { }{ } { }D Y Ei+ =12 ( ){ }Yi+12

( ){ } ( ){ } ( ){ } ( ){ }[ ]

{ } { } ( ){ } ( ){ } ( ){ }[ ]

Y Yh

Y Y

Y Y h Yh

Y Y

i i i i

i i i i i

+ +

+ +

= + +

= + + +

11 1 2

12

11

22

12

2

4

} Scrie: x Y { } ( ){ } ( ){ }Y Y i Ni i i i, , , , ,1 2 0= Stop

Cap. 9. Integrarea ecuatiilor diferentialeCaz particular n=1Caz particular n=2Algoritm

Curs_11_12_Cap9_Integr_EcDif.pdf

Documents

Transcript of Curs_11_12_Cap9_Integr_EcDif.pdf