REZOLVAREA NUMERIC A ECUAIILOR DIFERENIALE Introducere Acest tip de probleme provine din cadrul vast al analizei funcionale. Ecuaiile difereniale sau cu derivate

pariale constituie modelele matematice pentru majoritatea problemelor inginereti: studiul eforturilor la care sunt supuse elementele de rezisten: bare, grinzi, plci subiri, groase, conducte; studiul problemelor de cmp electric n dielectrici, cmp magnetic, cmp termic, propagarea undelor, curgerea fluidelor etc.

Odat stabilit fenomenul fizico-tehnic i ecuaiile difereniale care l guverneaz, ca form, coeficieni, condiii la limit (pe frontier) rmne de rezolvat ultima problem: rezolvarea acestui model matematic. Din diverse motive: neomogenitile fizice, frontiere cu geometrie dificil, numr de necunoscute, etc., rezolvarea o vom face cutnd o soluie aproximativ cu ajutorul unui cod numeric, folosind calculatorul.

O ecuaie diferenial este o ecuaie n care necunoscuta este o funce i n care intervine funcia

necunoscut, derivatele ei de diverse ordine i variabile independente de care depind aceste funcii. In cazul n care funcia necunoscut depinde de o singur variabil independent, ecuaia se numete

ecuaie diferenial ordinar, iar n situaia n care funcia necunoscut depinde de mai multe variabile independente, ecuaia se numete cu derivate pariale.

Ordinul unei ecuaii difereniale este cel mai nalt ordin a derivatei funciei necunoscute ce figureaz n

ecuaia respectiv. Expresia general a unei ecuaii difereniale, sub form implicit este:

0,...,,,, 22

=

n

n

dxyd

dxyd

dxdyyxF

unde . )(xyy = Forma explicit se poate scrie:

=

1

1

2

2

,...,,,, nn

n

n

dxyd

dxyd

dxdyyxf

dxyd

Privitor la condiiile la limit exist dou tipuri: (i). Condiii Cauchy: se cunosc ntr-un punct att valoarea funciei necunoscute ct i valorile derivatelor, pn la ordinul cel mai mare ce figureaz n ecuaie;

0x

(ii). Condiii la limit: se cunosc valorile funciei necunoscut n puncte diferite. Rezolvarea numeric a unei probleme asociate unei ecuaii difereniale poate fi privit sub dou aspecte:

(a). determinarea unei funcii )(~~ xyy = , aparinnd unei anumite clase de funcii (n general polinoame, dat fiind importana lor teoretic fundamental), i care aproximeaz suficient de bine soluia exact )(xyy = ; (b). determinarea valorilor aproximative ale soluiei exacte )(xyy = , ntr-o mulime de puncte date

. ,...1,0 , =ixi Vom expune n continuare principalele metode numerice al cror algoritm are un cost de calcul redus i se

preteaz la implementarea pe calculator pentru rezolvarea numeric a ecuaiilor difereniale. Pentru ecuaii difereniale ordinare acestea se pot clasifica n dou mari tipuri: - metode unipas (de tip Euler, Runge-Kutta) n care determinarea valorii aproximative a soluiei n fiecare

punct se va obine direct, pe baza informaiilor din punctul precedent;

- metode multipas, n care valoarea soluiei exacte n fiecare punct este aproximat folosind informaiile din mai multe puncte anterioare.

Evident este vorba de soluii aproximative pe care nu avem cum s le comparm cu o soluie exact, deoarece practic aceasta este imposibil de gsit.

De aceea n practic trebuie s procedm cu atenie pentru alegerea algoritmilor cei mai potrivii pentru problema concret de rezolvat.

1. METODE DE TIP EULER

1.1. Metoda Euler

Se consider ecuaia diferenial:

y = (x, y) (1) cu condiia iniial:

y(x0) = y0, (2) unde funcia este definit ntr-un domeniu D din planul xOy. Perechea (1),(2) constituie o problem Cauchy. Presupunem asigurate existena i uniciitatea soluiei. Se definete un cmp de direcii n D dac n fiecare punct M(x, y) D se ia direcia = arctg (x, y) ( fiind unghiul format de direcie cu sensul pozitiv al axei Ox).

Acest cmp de direcii are urmtoarea interpretare: graficul soluiei ecuaiei (1) cu condiia (2) trece prin punctul M(x0, y0) i este tangent n orice punct al su direciilor cmpului.

Metoda lui Euler propune aproximarea soluiei printr-o linie poligonal n care fiecare segment este coliniar cu direcia cmpului definit de extremitatea sa stng. Astfel se consider nodurile echidistante xi = x0 + ih, n0i ,= . n punctul M0(x0, y0) se calculeaz direcia cmpului definit de M0 i se scrie ecuaia dreptei determinate de M0 i de aceast direcie:

y = y0 + (x0, y0)(x x0) (3) Funcia (3) se propune ca aproximant a soluie problemei (1)+(2) pe [x0,x1]. Valoarea aproximativ a

soluiei n x1 este dat de: y1 = y0 + (x0, y0)(x1 x0) = y0 + oh Repetm procedeul i presupunnd c n xi s-a calculat valoarea aproximativ yi, atunci pe intervalul [xi,

xi+1] se aproximeaz soluia cu: y = yi + (xi, yi)(x xi) = yi + i(x xi),

iar n punctul xi+1 se obine valoarea aproximativ: yi+1 = yi + hi.

Aproximarea este justificat i de urmtoarea teorem:

Teorem.

a) Dac: y C2[a, b], atunci 1)( (x, x + h) cu proprietatea )(''y2h)x('y

h)x(y)hx(y

1+=+

b) Dac y C2[a, b], atunci 2)( (x h, x) cu proprietatea )(''y2h)x('y

h)hx(y)x(y

2=

Demonstraie. Din formula lui Taylor avem:

y(x + h) = y(x) + )(''y2

h)x('hy 12

+ , 1 (x, x + h)

y(x h) = y(x) )(''y2

h)x('hy 22

+ , 2 (x h, x)

)(''y2h)x('y

h)x(y)hx(y

1+=+

)(''y2h)x('y

h)hx(y)x(y

2= Deci

=)x('y )h(h

)x(y)hx(y1++ ; =)x('y )h(h

)hx(y)x(y2+

Considernd cunoscut aproximarea yi a soluiei problemei (1)+(2) n xi, procedeul de aproximare Euler, poate fi acum rezumat astfel:

(4) ,...1,0),(

1

1 =

+=+=

=

+

+ ihfyyhxxyxff

iii

ii

iii

Observaii: (i).Neglijarea termenilor de ordin superior n (4) face ca metoda s fie comod n calcul, dar puin precis, erorile cumulndu-se la fiecare pas. (ii).Metoda se poate aplica i dac nodurile xi nu sunt echidistante, avnd la fiecare iteraie alt pas h n acest caz. In general o metod unipas poate fi scris sub forma:

);,(1 hyxhyy iiii +=+Pentru soluie a problemei Cauchy (1)+(2) (i pentru orice )(xyy = }1,...,2,1,0{ ni ) considerm:

( ) ));(,()()(1)( 1 hxyxxyxyhh iiiii = + Cantitatea de mai sus se numete eroare de consisten a metodei n i reprezint o msur a calitii metodei de aproximare.

ix

Pentru metoda Euler, innd cont de algoritmul acesteia, de expresia erorii de consisten i de faptul c

obinem c ),();,( iiii yxfhyx ),( ),(''

21)( hxxhyhi +=

In consecin eroarea de consisten este de ordinul . )(hO

Exemplul 1 Se consider problema Cauchy:

[ ]

==

0,1 x,1)0(yyx2y'y

Ne propunem s determinm o soluie aproximativ a acestei probleme folosind metoda Euler cu pasul h = 0,2.

Rezolvare Folosind formulele (4) pentru x0 = 0; x1 = 0,2; x2 = 0,4;

x3 = 0,6; x4 = 0,8; x5 = 1 i y0 = 1, obinem:

y1 = y0 + hf(x0, y0); f(x0, y0) = y0 0

0

yx2

= 1.

Deci y1 = 1 + 0,2 1 = 1,2.

y2 = y1+0,2f(x1, y1) = y1+0,2

1

11 y

x2y = 1,2 + 0,20,8667 =

= 1,3733. Obinem n final urmtorul tabel, ultima coloan reprezentnd valorile exacte ale soluiei problemei

propuse (y = 1x2 + ).

xi yi(Euler) yi (exact) 0 0,2 0,4 0,6 0,8 1

1 1,2 1,3733 1,5294 1,6786 1,8237

1 1,1832 1,3416 1,4832 1,6124 1,7320

n continuare, vor fi prezentate cteva variante mbuntite ale algoritmului Euler. 1.2. Metoda Euler modificat Considerm pe [xi, xi+1] ca direcie a segmentului MiMi+1 direcia definit de punctul de la mijlocul

segmentului (nu de extremitatea stng ca n formula iniial) se obine metoda Euler modificat. Dac xi, yi sunt valori calculate, procesul iterativ este urmtorul:

+=

=+=

+==+=

++++++

+

;hfyy ;y,xff ;f2hyy

;2hxx );y,x(ff ;ihxx

21ii1i

21i

21i

21iii

21i

i21iiii0i

Pentru aceast metod )(hi este de ordinul . )( 2hO Exemplul 2 Vom aplica metoda Euler modificat n rezolvarea aceleiai probleme de la exemplul 1 pentru a putea face

o comparaie rapid privind eficiena sporit a metodei. Deci:

[ ]1,0x ,1)0(y

yx2y'y

==

, h = 0,2

Rezolvare x0 = 0; x1 = 0,2; x2 = 0,4; x3 = 0,6; x4 = 0,8; x5 = 1;

f(x,y) = y yx2

x1/2 = x0 + 2h = 0,1; y1/2 = y0 + 2

h f(x0, y0) = 1,1

f(x1/2, y1/2) = 0,9182 = f1/2 y1 = y0 + hf1/2 = 1 + 0,2 0,9182 = 1,1836

Continund n acelai mod obinem tabelul:

xi yi (Euler modificat) yi (exact) 0, 0,2

1 1,1836

1 1,1832

0,4 0,6 0,8 1

1,3426 1,4850 1,6152 1,7362

1,3416 1,4832 1,6124 1,7320

1.3. Metoda EulerCauchy Introducnd y ca o nou variabil, se poate nlocui ecuaia diferenial (1) cu urmtorul sistem de ecuaii:

==

).y,x(f'y

'ydxdy

Integrnd prima ecuaie a sistemului de la xi la xi+1, unde xi = x0 + ih se obine:

=

+= +

+

)y,x(fy

;dx'yyy

ii'i

1x

xi1i

i

i

Folosind metoda dreptunghiului pentru aproximarea integralei din prima ecuaie obinem: yi+1 = yi + hyi = yi + h(xi, yi). (5) Dac aplicm metoda trapezului pentru aproximarea aceleiai integrale se obine:

yi+1 = yi + )(2'

1'

++ ii yyh = yi + ))y,x(f)y,x(f(2h

1i1iii +++ =

= yi + )(2 1++ ii ffh . (6)

Ecuaia (6) poate fi rezolvat iterativ, relativ la necunoscuta yi+1. Dac se alege ca aproximaie iniial dat de (5), obinem pentru yi+1 valoarea:

)0(1+iy

)),(),((2

)0(11

)1(1 +++ ++= iiiiii yxfyxfhyy . (7)

Deci pentru construirea liniei poligonale care aproximeaz soluia problemei (1) + (2) este dat algoritmul EulerCauchy:

,...1,0

).(2

);,(

;

;

)0(1

)1(11

)0(11

)0(1

)0(1

1

=

++===

+=+=

+++

+++

+

+

i

ffhyyy

yxffhfyyhxx

iiiii

iii

iii

ii

Pentru aceast metod )(hi este de ordinul . )( 2hO

Observaie: In scopul ameliorrii aproximaiilor se pot utiliza metodele prezentate de o manier iterativ. Astfel, metode euler-Cauchy ne conduce la urmtorul algoritm de calcul:

;)0( 1 iii hfyy +=+ )(

2)1(

1)(1

++ ++= kiiiki ffhyy , k = 1, 2,

cu ( ))1k( 1i1i)1k( 1i y,xff +++ = Se demonstreaz c dac funcia este lipsichitzian n y, de constant L i dac h este suficient de mic cu

hL < 2 atunci , pentru k . 1)( 1 ++ iki yy

n practic, se continu iteraiile dup indicele k pn cnd un criteriu de eroare impus de exemplu:

)k(1i

)1k(1i

)k(1i

yyy

+

++ < , impus, 0, este satisfcut. + )k( 1iy

Exemplul 3 S se gseasc soluia aproximativ a problemei Cauchy:

[ ]

=+=

1,0x,1)0(yyx'y

folosind metoda Euler-Cauchy i o precizie de 10-4. Rezolvare Vom folosi pasul h = 0,05; x0 = 0; x1 = 0,05; y0 = 1;

f(x, y) = x + y; )0(

1y = y0 + hf(x0, y0) = 1 + 0,05 1 = 1,05 )1(

1y = y0 + 2h ( )( ))0(1100 y,xf)y,x(f + =

= 1 + 205,0 (1 + 1,1) = 1,0525

( ) ( )( ))1(11000)2(1 y,xfy,xf205,0yy +++ = = 1 +

205,0 (1 + 1,1025) = 1,0525625

Deci fiind stabile primele 4 zecimale rezult: y1 = 1,0525. Similar cu x2 = 0,1 obinem:

)y,x(hfyy 111)0(

2 += =1,1077 )1(

2y = 1,11036; = 1,11042. )2(

2yDeci y2 = 1,1104. Soluia exact fiind y = 2ex x 1, obinem n x = 0,1 valoarea: y(0,1) = 1,1103 eroarea fiind de 0,0001.

1.4. Metoda EulerHeun

Vom ncheia trecerea n revist a variantelor Euler pentru rezolvarea ecuaiilor de tipul (1) prezentnd o

ultim variant de ordinul doi a metodei lui Euler i anume aceea propus de Heun, [26]. Presupunnd calculat valoarea yi la pasul xi, se propune pentru calcularea soluiei la pasul xi+1 expresia:

41hyy ii +=+

+++ iiiiii yxhfyhxfyxf ,(3

2,323),( .

2. Metode de tip Runge-Kutta

Metodele de tip Euler prezentate sunt explicite i nu necesit valori de start. Faptul c au un ordin sczut

al erorii de consisten conduce la o aplicabilitate limitat. In scopul obinerii unor metode de ordin ridicat trebuie renunat fie la proprietatea de a fi unipas i pastrat liniaritatea, fie viceversa. Metodele de tip Runge-Kutta sunt neliniare i conserv caracteristicile metodelor unipas, avnd un ordin ridicat. Ele au trei proprieti principale:

1. Sunt metode directe, adic pentru determinarea aproximrii soluiei la pasul i+1 avem nevoie de

informaiile existente n punctul precedent xi, yi. 2. Sunt identice cu seriile Taylor pn la termenii hn, unde h este pasul curent iar n este diferit pentru

metode diferite din aceast familie i definete ordinul metodei. 3. n procesul de calcul nu necesit dect evaluarea funciei din memebrul drept pentru diverse valori x i y . Nu este nevoie de calculul derivatelor acesteia.

Metodele de tip Euler pot fi i ele incluse n familia Runge-Kutta i putem astfel observa c metoda Euler

este o metod R-K de ordinul nti iar metodele Euler-Cauchy i Euler-Heun sunt metode R-K de ordinul 2.

2.1. Construirea formulelor Runge-Kutta Ne ocupm n special de rezolvarea ecuaiilor difereniale ordinare cu condiii iniiale. Adic pentru

ecuaii difereniale ordinare de ordinul nti ne intereseaz rezolvarea problemei Cauchy:

==

00 y)x(y))x(y,x(f)x('y

, x [a, b], x0 = a (8) Metodele Runge-Kutta constau n aproximarea soluiei problemei (8) astfel:

=

+=+m

iii hkcxyhzhxy

1)()()()(

unde

+++++=

++==

11,2211

12122

1

...,()(........................

),()(),()(

mmmmmmm kbkbkbyhaxhfhk

kbyhaxhfhkyxhfhk

i mibai

jiji ,...,3,2 ,

1

1==

=

c1, ci, ai, bij urmnd a fi determinai (i = m,2 ; j = 1i,1 )

Observaie: Pentru a fi consisente, metodele Runge-Kutta trenuie s satisfac condiia: . =

=m

iic

11

Notm (h) = y(x+h) z(h) (eroarea de aproximare). Vom determina parametrii c1, ci, ai, bij din condiiile:

====

+ f funcie anumit opentru ,0)0(

i f funcie oricepentru ,0)0(...)0(')0()1p(

)p(

Din formula lui Taylor n 0 avem:

(h) = =

++ ++p

0k

)1p(1p

)k(k

)h()!1p(

h)0(!k

h

(h) = )h()!1p(

h )1p(1p +++ , 0 < < 1 (9)

expresie ce indic ordinul de mrime al erorii de aproximare (deci O(hp)); Cazuri particulare

1) Dac m = 1, atunci: (h) = y(x+h) y(x) c1k1(h) = y(x+h) y(x) c1hf(x,y) Se verific uor c (0) = 0 pentru orice funcie f. (h) = y(x+h) c1f(x,y) Obinem (0) = y(x) c1f(x,y) = (1 c1)f(x,y). Se observ c pentru c1 = 1, (0) = 0 pentru orice funcie f. (h) = y(x+h), deci (0) = y(x) = fx(x,y) + f(x,y) fy(x,y) Se observ c (0) 0 pentru o anumit funcie f. De exemplu, pentru f(x,y) = y se obine (0) = y. Deci p = 1.

Astfel, pentru c1 = 1, se obine formula: y(x+h) = y(x) + hf(x,y)

care, este formula din metoda Euler. Eroarea este de ordinul lui h.

2) Dac m = 2, atunci: (h) = y(x+h) y(x) c1k1(h) c2k2(h) = y(x+h) y(x)

c1hf(x,y) c2hf(x + a2h, y + b21k1) Se verific uor c (0) = 0 pentru orice funcie f. (h) = y(x+h) c1f(x,y) c2f(x + a2h, y + b21k1)

-c2h [ + )kby,hax(fa 1212hax2 2 +++ )kby,hax(f)y,x(fb 1212kby21 121 +++ ] Notm cu: expresia inclus n parantezele drepte. /hf )kby,hax( 1212 ++

Obinem: (0) = y(x) c1f(x,y) c2f(x,y) = (1 c1 c2)f(x,y). (h) = y(x+h) 2c2[ (x + a2h, y + b21k1)] hf

c2h[ (x + a2h, y + b21k1)] . hhf Obinem: (0) = y(x) 2c2[a2fx (x,y) + b21f(x,y) fy(x,y)] =

=(fx + ffy) 2c2(a2fx + b21f fy) = (1 2c a )f + (1 2c b )ff 2 2 x 2 21 y[ ])kby,hax(f 1212hh(h) = y(x+h) 3c2 ++ c2h [ ])kby,hax(f 1212hhh ++

Obinem: (0) = y(x)3c2( fxx+a2b21ffxy+a2b21ffyx+ f2fyy) = 22a 221b= (fxx + ffxy + ffyx + fxfy + ff2y + f2fyy) 3c2( fxx + a2b21ffxy+ a2b21ffyx + f2fyy) = (1 3c2 ) fxx + (1 3c2a2b21) ffxy+

22a

221b

22a

+ (1 3c2a2b21)ffyx + (1 3c2 ) f2 fyy + fxfy + f 221b2yf

Dac:

===

0bc210ac210cc1

212

22

21

sau

===+

2212

21

c21ba

1cc (10)

atunci (0) = (0) = (0) = 0 pentru orice funcie f. Sistemul (10) este compatibil nedeterminat. De asemenea, (0) nu este identic nul pentru orice funcie f. De exemplu, pentru funcia f(x,y) = y, se

obine (0) = y. Deci p= 2. Pentru m = p = 2 se pot obine oricte formule dorim, alegnd parametrii c1, c2, a2, b21 astfel nct s

verifice (10). Din (9) rezult c eroarea n aceste formule este de ordinul lui h2.

Astfel, o soluie a sistemului (10) este: c1 = c2 = 21 , a2 = b21= 1, i pentru aceast soluie se obine formula

Runge-Kutta simpl:

++==

++=+

)ky,hx(hfk)y,x(hfk

)kk(21)x(y)hx(y

12

1

21

(10)

Observaie. Dac f nu depinde de y obinem:

[ ])hx(f)x(f2h)x(y)hx(y ++=+

Pe de alt parte:

+ +==+ hxx

hx

x

dt)t(fdt)t('y)x(y)hx(y

Din cele dou relaii de mai sus rezult formula trapezului (vezi 5.1.).

[ ]+ ++=hxx

)hx(f)x(f2hdt)t(f

O alt formul Runge-Kutta de ordinul 2 se obine pentru

c1 = 41 , c2 = 4

3 , a2 = b21 = 32 i anume:

++==

++=+

12

1

21

k32y,h

32xhfk

)y,x(hfk

)k3k(41)x(y)hx(y

(11)

cunoscut ca formula Euler-Heun. Pentru calculul aproximativ al soluiei unei probleme Cauchy (8) pe un interval [a, b] de noduri

echidistante de pas h procedm astfel : - presupunem cunoscute valorile xi, yi (yi y(xi)) de exemplu, pentru i = 0 cunoatem din (8) x0, y0; - calculm, folosind de exemplu (10) xi+1 = xi + h

yi+1 = yi + 21 (k1 + k2), unde

k1 = hf(xi, yi) k2 = hf(xi + h, yi + k1) Exemplul 4

Se consider problema Cauchy:

==

1)0(yxy'y

Folosind o formul Runge-Kutta de ordinul 2, s se aproximeze soluia problemei date n xi = 0.1i, 1 i 3.

Rezolvare f(x,y) = xy; x0 = 0; y0 = 1; h = 0.1 Se observ c soluia exact este y(x) = . Deci 2/x

2e

y(0,1) 1,005013 y(0,2) 1,020201 y(0,3) 1,046027 Vom ntocmi urmtorul tabel:

I 0 1 2 3 xi 0 0.1 0.2 0.3 yi 1 1.005 1.0201755 1.0459859f(xi,yi) 0 0.1005 0.2040351 k1 = hf(xi,yi) 0 0.01005 0.0204035 xi + h 0.1 0.2 0.3 yi + k1 1 1.01505 1.040579 f(xi + h, yi + k1) 0.1 0.20301 0.3121737 k2= hf(xi+h, yi+k1) 0.01 0.020301 0.0312174 (k1 + k2) /2 0.005 0.0151755 0.0258104

3) Pentru m = 3, avem: (h) = y(x + h) y(x) c1k1 c2k2 c3k3, unde:

k1 = hf(x,y); k2 = hf(x+a2h, y+b21k1); k3 = hf(x+a3h, y+b31k1 + b32k2) Efectund calcule asemntoare ca n cazul 2) se constat c (0) = (0) = (0) = (0) = 0, f, dac:

==+=

=+=+=++

61bac;ba;bba

31acac;

21acac;1ccc

322321232313

233

2223322321

care este de asemenea un sistem compatibil nedeterminat cu o infinitate de soluii. O soluie a acestui sistem este:

c1 = 61 ; c2 = 3

2 ; c3 = 61 ; a2 = 2

1 , a3 = 1

b21 = 21 , b31 = -1, b32 = 2

obinndu-se urmtoarea formul Runge-Kutta de ordinul 3:

++=

++==

+++=+

)2,(2

,2

),(

)4(61)()(

213

12

1

321

kkyhxhfk

kyhxhfk

yxhfk

kkkxyhxy

(12)

O alt formul tot de ordinul 3 este urmtoarea:

++=

++==

++=+

23

12

1

31

32,

32

3,

3

),(

)3(41)()(

kyhxfhk

kyhxfhk

yxfhk

kkxyhxy

(13)

i n acest caz se constat c nu exist formule Runge-Kutta cu m = 3 i p = 4, deoarece IV(0) nu este identic nul pentru orice funcie f. ntr-adevr pentru f(x,y) = y se obine IV(0) = y.

Observaii. (i). Din (9) rezult c eroarea n aceste formule este de ordinul lui h3. (ii).Dac funcia f nu depinde de y, atunci

+hxx

dt)t(f = +hxx

dt)t('y =y(x+h)y(x) = ( )

++

++ hxf2hxf4)x(f

6h

adic formula de integrare numeric Simpson. 4) m = 4. Raionamente i calcule asemntoare celor anterioare conduc la construirea unei mulimi de

formule Runge-Kutta de ordinul 4, din care mai cunoscute sunt:

++=

++=

++==

++++=+

),( 21,

21

2,

2

),(

)22(61)()(

34

23

12

1

4321

kyhxfhk

kyhxfhk

kyhxfhk

yxfhk

kkkkxyhxy

(14)

formula propriu-zis Runge-Kutta i

+++=

++=

++==

++++=+

),( 31,

32

3,

3

),(

)33(81)()(

3214

213

12

1

4321

kkkyhxfhk

kkyhxfhk

kyhxfhk

yxfhk

kkkkxyhxy

(15)

formula Kutta-Simpson Observaie. Aceste formule sunt de ordinul patru. .

Pentru calculul aproximativ al soluiei unei probleme Cauchy (8), pe un interval [a, b], cu pasul h procedm astfel:

x0 = a, xi = x0 + ih, i = N,0 cu xN = b , deci h = Nab

cunoscnd pentru xi, valorile yi (yi y(xi)) de exemplu pentru i = 0, cunoatem x0, y0. Calculm: xi+1 = xi + h

yi+1 = yi + 61 (k1 +2k2 + 2k3 + k4)

k1 = hf(xi, yi)

k2 = hf

++2ky,

2hx 1ii

k3 = hf

++2

ky,2hx 2ii

k4 = hf(xi + h, yi +k3) Exemplul 5

Se consider problema Cauchy: .

==

1)0(yxy'y

S aproximm soluia acestei ecuaii n punctele xi = ih, i = 9,1 , h = 0.1, folosind o metod Runge-Kutta de ordin 4. xi yi valoarea soluiei exacte 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1.005013 1.020202 1.046028 1.083287 1.133148 1.197217 1.277621 1.377128 1.499303

1.005013 1.020201 1.046027 1.083286 1.133148 1.197217 1.277620 1.377128 1.499303

y(x) = 2/x

2e

2.2. Metod de tip Runge-Kutta implicit Metodele de tip Runge-Kutta, expuse anterior, sunt explicite. Pentru a ameliora proprietile de

stabilitate ale acestor metode se consider cele de tip implicit.

Prin proprieti de stabilitate ne referim la restriciile impuse asupra pasului de integrare n situaia

utilizrii metodei respective.

Forma general a unei astfel de metode Runge-Kutta implicit de tip este: m

=

+

=+=

m

jjjii

iiii

kchyx

hyxhyy

1

1

);,(

);,(

mjba

mjkbhyhaxfk

m

sjsj

m

ssjsjj

,...,2,1 ,

,...,2,1 ,,(

1

1

==

=++=

=

=

Comparativ cu metodele explicite, funciile kj numai sunt definite explicit, ci printr-un set de m

ecuaii implicite, n general neliniare. n practic, se folosete cazul metodelor Runge-Kutta implicite cu

= 2

m

m

Astfel:

kj = f(x + haj, y + bj1hk1 + bj2hk2), j = 2,1

Considernd k1, k2 dezvoltabile sub forma

kj = Aj + hBj + h2Cj + h3Dj + O(h4), j = 2,1

pe baza dezvoltrii n serie Taylor n raport cu (x, y) a lui kj, j = 2,1 se obine prin identificarea puterilor lui

h:

==

=

;41bb

;63

21a

2211

1

41ab;

41ab

63

21a

221112

2

==

= m (16)

Observaie. Prin simetrie se observ faptul c valorile date prin alternarea semnelor, n relaiile

(16), conduc la aceeai metod.

n concluzie, o metod Runge-Kutta implicit de ordinul patru, este dat de formulele:

+

+

+=

+++

++=

++=+

212

211

211

41

63

41,

63

21

63

41

41,

63

21

)(2

kkyhxfk

kkyhxfk

kkhyy

ii

ii

ii

(17)

Remarc. Este important precizarea existenei unei metode Runge-Kutta de tip semi-explicit de ordinul

patru, descris de relaiile de mai jos:

( )( )

( )

++=

+++==

+++=+

23

212

1

3211

,41

41,

21

,

46

kyhxfk

kkyhxfk

yxfk

kkkhyy

ii

ii

ii

ii

(18)

3. Metode numerice multipas

Fie problema Cauchy

(19) ],[ ],,[ )(

),()('0

00

baxbaxyxy

yxfxy

==

i diviziunea intervalului dat de: ],[ ba bxxxa n

Cele mai cunoscute metode implicite, multipas, sunt cele de tip Adams-Moulton: (i). de ordinul trei, cu doi pai:

,...2,1

)),(),(8),(5(12 11111

=++= +++

i

yxfyxfyxfhyy iiiiiiii

(24) (ii). de ordinul patru, cu trei pai:

,...3,2

)),(),(5),(19),(9(24 2211111

=+++= +++

i

yxfyxfyxfyxfhyy iiiiiiiiii

(25)

4. Metode numerice de tip predictor-corector

O metod nmeric de tip predictor-corector este o combinaie ntre o metod numeric explicit i o

metod numeric implicit. Metoda explicit permite predicia unei valori aproximative i cea implicit corecteaz (odat sau de mai multe ori) aceast predicie.

Metodele predictor-corector, ofer o precizie superioar fa de metodele prezentate n paragrafele anterioare, fr a implica sporiri ale numrului de operaii aritmetice.

n ceea ce privete comportarea metodelor de tip predictor-corector, acestea au eroarea de procedeu mai mic, dar sunt puternic afectate de eventualele erori ale valorilor de pornire necesare algoritmului.

Cea mai simpl metod de acest tip este cea prezentat sub numele de Euler-Cauchy. Astfel pentru predictorul ne d valoarea ihxxi += 0 ),( ,)0( 1 yxfffyy iiii =+=+ i corectorul o corecteaz cu formula

,...2,1 ),(2

)1(1

)(1 =++= ++ kffhyy kiiiki

Criteriul de oprire este:

++ )1( 1)( 1 kiki yy sau

+

++

)(1

)1(1

)(1

ki

ki

ki

y

yy cu i 0)( 1 +kiy precizia de calcul impus.

4.1 Metoda lui Milne

Algoritmul de calcul este urmtorul:

( )( )

+++=

++=+=

++

+

+

)0(1111

123)0(1

1

43

223

4

iiiii

iiiii

ii

fffhyy

fffhyy

hxx

(26)

unde . ),( )0( 11)0(

1 +++ = iii yxff Observatii:

(i). este o metod de ordinul ; )( 4hO(ii). In scopul ameliorrii rezultatelor putem aplica formulele de calcul de manier iterativ:

( ) ,...2,1 ,43

)1(111

)(1 =+++= ++ kfffhyy ikiiiki

4.2 Metoda Adams-Bashforth-Moulton de ordinul patru

Pentru problema Cauchy

],[ ],,[ )(

),()('0

00

baxbaxyxy

yxfxy

==

predictorul este representat de metoda Adams-Bashforth de ordinul patru:

,...4,3 )),(9),(37

),(59),(55(24

3322

11)0(1

=+++=

+

iyxfyxf

yxfyxfhyy

iiii

iiiiii

(27)

Remarc: Pentru determinarea valorilor utilizm o metod de ordinul patru. iiii yyyy ,, 1,23 Corectm aproximarea genert mai sus prin intermediul corectorului dat de metoda Adams-Moulton de

ordinul patru:

,...2,1 ,...,3,2

)),(),(5),(19),(9(24 2211

)1(11

)(1

==+++= +++

ki

yxfyxfyxfyxfhyy iiiiiik

iiik

i

BIBLIOGRAFIE 1. R. Burden, J. Faires, Numerical Analysis, PWS-Kent, 2001.

2. C. Carasso, Analyse Numrique, Lidec, Canada, 1970. 3. B. Demidovitch. I. Maron, lments de Calcul Numrique, Mir, Moscow, 1973. 4. D. Ebnc, Metode de Calcul Numeric, Editura Sitech, Craiova, 1994. 5. R. Militaru Mthodes Numriques. Thorie et Applications Ed. Sitech, 2008 6. J.P. Nougier, Mthodes de Calcul Numrique, Hermes Sciences Publication, Paris, 2001. 7. M. Popa, R. Militaru, Metode Numerice algoritmi si aplicaii, Ed. Sitech, Craiova, 2007. 8. M. Popa, R. Militaru, Analiz numeric note de curs, Editura Sitech, Craiova, 2003.

1. REZOLVAREA NUMERICA A SISTEMELOR DE ECUATII DIFERENTIALE

Metodele numerice uilizate pentru o singur ecuaie diferenial, anterior prezentate, se pot extinde i n cazul sistemelor de ecuaii difereniale. Considerm situaia metodelor numerice de tip Runge-Kutta. Metodele Runge-Kutta, prezentate pentru ecuaii difereniale, pot fi aplicate cu uurin i la sisteme de ecuaii difereniale. Fie sistemul de ecuaii difereniale:

(1)

=

=

)y...,y,x(fy................................),y,...,y,x(fy

n1nn

n111

i urmrim s determinm soluia care satisface condiiile iniiale: yi(x0) = yi,0 , i = n,1 (2)

Presupunnd c dispunem de soluia problemei (1) + (2) la pasul i: y1,i , ..., yn,i metoda Runge-Kutta de ordinul 4 calculeaz soluia n pasul i+1 cu formulele:

y1,i+1 = y1,i + y1,i ...........................

(3) yn,i+1 = yn,i + yn,i Coreciile: y1,i , y2,i , ... , yn,i care intervin n formulele precedente se calculeaz cu ajutorul relaiilor: y1,i = 14131211 k6

1k31k

31k

61 +++ ;

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

yn,i = n4n3n2n1 k61k

31k

31k

61 +++ ;

iar coeficienii au urmtoarea form: n4n1

14

11 k,...,k,...,k,...,k

)y....,y,x(hfk i,ni,1ijj1 = , j = n,1 ;

+++=

2ky,...,

2ky,

2hxhfk

n1

i,n

11

i,1ijj2 , j = n,1 ;

+++=

2ky,...,

2ky,

2hxhfk

n2

i,n

12

i,1ijj3 , j = n,1 ;

( )n3i,n13i,1ijj4 ky,...,ky,hxhfk +++= , j = n,1 ;

Exemple Folosind o formul Runge-Kutta de ordinul 4, s se rezolve problema Cauchy:

==

+==

1)1(z1)1(y

)xy/()yx(zzy/xy

222

pentru x [0, 1]. Rezolvare Considerm 10 noduri echidistante xi pe [0, 1], de pas h = 0.1. Din (14), pentru fiecare i = 0, 1, , 9, avem :

yi+1 = yi + ( )14131211 kk2k2k61 +++

zi+1 = zi + ( )24232221 kk2k2k61 +++ unde:

)z,y,x(hfk iiijj1 = ; j = 1, 2

)2kz,

2ky,

2hx(hfk

21

i

11

iijj2 +++= ; j = 1, 2

)2kz,

2ky,

2hx(hfk

22

i

12

iijj3 +++= ; j = 1, 2

)kz,ky,hx(hfk 23i13iij

j4 +++= ; j = 1, 2

Obinem urmtoarele valori:

x1 = 1.1 y1 = 1,1000000000; z1 = 1,2099979386 x2 = 1.2 y2 = 1,2000000000; z2 = 1,4399959714 x3 = 1.3 y3 = 1,3000000000; z3 = 1,6899940413 x4 = 1.4 y4 = 1,4000000000; z4 = 1,9599921095 x5 = 1.5 y5 = 1,5000000000; z5 = 2,2499901493 x6 = 1.6 y6 = 1,6000000000; z6 = 2,5599881417 x7 = 1.7 y7 = 1,7000000000; z7 = 2,8899860729 x8 = 1.8 y8 = 1,8000000000; z8 = 3,2399839328 x9 = 1.9 y9 = 1,9000000000; z9 = 3,6099817135 x10 = 2 y10 = 2,0000000000; z10 = 3,9999794092

soluia exact fiind , deci

==

2x)x(zx)x(y

pentru x1 = 1,1 y(x1) = 1,1; z(x1) = 1,21 x2 = 1,2 y(x2) = 1,2; z(x2) = 1,44 x3 = 1,3 y(x3) = 1,3; z(x3) = 1,69 x4 = 1,4 y(x4) = 1,4; z(x4) = 1,96 x5 = 1,5 y(x5) = 1,5; z(x5) = 2,25 x6 = 1,6 y(x6) = 1,6; z(x6) = 2,56 x7 = 1,7 y(x7) = 1,7; z(x7) = 2,89 x8 = 1,8 y(x8) = 1,8; z(x8) = 3,24 x9 = 1,9 y(x9) = 1,9; z(x9) = 3,61 x10 = 2,0 y(x10) = 2,0; z(x10) = 4,00

2. REZOLVAREA NUMERICA A ECUATIILOR DIFERENTIALE DE ORDIN SUPERIOR Fie ecuaia diferenial de ordinul scris sub form explicit n

(4) ))(),...,('),(,()( )1()( xyxyxyxfxy nn =

cu condiiile asociate

(5)

=

==

10)1(

10

00

)(..................

)(')(

nn txy

txytxy

Introducem notaiile

(6) n1,2,..,j ),()( )1( == xyxz jj Remarc: 1,..,2,1 ),()( 1

' == + njxzxz jj Tinnd cont de (4), (5) i (6) avem:

(7)

==

==

))(),...,(,()(

)()(......................)()(

)()(

1'

'1

3'2

2'1

xzxzxfxzxzxz

xzxzxzxz

nn

nn

respectiv:

(8)

=

==

10

102

001

)(..................)()(

nn txz

txztxz

In concluzie (7) i (8) reprezint un sistem avnd n ecuaii difereniale ordinare cu n condiii iniiale. Putem astfel s l rezolvm cu tehnica prezentat mai sus. 7. Folosind o formul Runge-Kutta de ordinul 4 s se rezolve problema Cauchy:

==

+=

1)0('y1)0(y

'yyx''y

pentru x [0, 1].

Rezolvare Transformm ecuaia diferenial de ordinul doi n sistemul de ordinul nti:

=+=====

1)0(z ;zyx)z,y,x(f'z1)0(y ;z)z,y,x(f'y

2

1

Considerm 10 noduri echidistante xi pe [0,1] de pas h = 0.1. Din (14), pentru fiecare i = 0, 1, , 9 avem:

( )14131211i1i kk2k2k61yy ++++=+ ( 24232221i1i kk2k2k61zz ++++=+ ) - variabile auxiliare

unde )z,y,x(hfk iiij

j1 = ; j = 1,2

)2

kz,2ky,

2hx(hfk

21

i

11

iijj2 +++= ; j = 1,2

)2

kz,2

ky,2hx(hfk

22

i

12

iijj3 +++= ; j = 1,2

)kz,ky,hx(hfk 23i13iij

j4 +++= ; j = 1,2

Obinem urmtoarele valori : x1 = 0.1 y1 = 0.90968333333; x2 = 0.2 y2 = 0.83758567267; x3 = 0.3 y3 = 0.78223901323; x4 = 0.4 y4 = 0.74247458668; x5 = 0.5 y5 = 0.71738325383; x6 = 0.6 y6 = 0.70628213912; x7 = 0.7 y7 = 0.70868659523; x8 = 0.8 y8 = 0.72428672347; x9 = 0.9 y9 = 0.75292779297; x10 = 1 y10 = 0.79459400090;

3. METODE NUMERICE CU PAS VARIABIL PENTRU REZOLVAREA NUMERICA A ECUATIILOR DIFERENTIALE

Sunt tehnici utilizate pentru controlul erorii unei metode numerice pentru rezolvarea numeric a unei

probleme de tip Cauchy:

=+=

00

00

)(],[ )),(,()('

yxyhxxxxyxfxy

fcnd apel la o manier eficient de alegere a mulimii de puncte n care se va aproxima soluia exact . )(xyIdeea de baz este utilizarea unor metode numerice de ordine diferite n scopul caracterizrii erorii de consisten alegerii pasului de integrare astfel nct eroarea global s fie inferioar unei anumite precizii impuse. Considerm n continuare dou metode numerice care conduc la aproximaiile urmtoare ale soluiei exacte

n : )(xy hxx ii +=+1 (a) 0 ),,,(1 >+=+ ihyxhyy iiii

0 ),,~,(~~ 1 >+=+ ihyxhyy iiii (b)

Presupunem c iii yyxy ~)( = i c (a) i (b) sunt obinute cu acelai pas . h

Atunci:

)()),(,()()( ),,()()(

11

111

hhhxyxhxyxyhyxhyxyyxy

iiiii

iiiiii

++

+++=

=

Tinnd cont c i c )()(1

ni hOh =+ )()(~ 11 ++ = ni hOh rezult

( 111 ~1)( +++ iii yyhh ) (c) de unde

(d) ni khh + )(1k constant ce nu depinde de . h

In consecin pentru un nou pas de integrare : qh

)()( 11 hqqh in

i ++ Impunnd + )(1 qhi obinem:

n

ii yyhq

/1

11~

++

Un exemplu de metod numeric de acest tip este metoda Runge-Kutta-Fehlberg. Ea const n utilizarea unei metode Runge-Kutta de ordin 5:

654311 552

509

5643028561

128256656

135126~ kkkkkyy ii ++++=+

n scopul estimrii erorii de consisten pentru o metod Runge-Kutta de ordin 4:

54311 51

41042197

25651408

21625~ kkkkyy ii +++=+ unde

+++=

+++=

+++=

+++=

++==

543216

43215

3214

213

12

1

4011

41041859

256535442

278,

2

4104845

51336808

216439,

21977296

21977200

21971932,

1312

329

323,

83

41,

4

),(

kkkkkyhxfhk

kkkkyhxfhk

kkkyhxfhk

kkyhxfhk

kyhxfhk

yxfhk

ii

ii

ii

ii

ii

ii

4. ECUATII DIFERENTIAL-ALGEBRICE

Sub form implict acestea se scriu:

0))('),(,( =xyxyxF unde funcia necunoscut este o funcie scalar sau vectorial. )(xyy =

Vom aborda n continuare cazul ecuaiilor diferenial-algebrice care sub form explicit devin:

)()()()(')( xfxyxAxyxM += Exemplu:

Dndu-se ecuaia diferenial-algebric , aceasta se poate explicita sub forma

, unde , ,

=++=

+=

2424'

467'

wvuvuv

xvuu

)()()()(')( xfxyxAxyxM +=

=

000010001

M

=111924067

A

=

240

4xf

Ecuaiile diferenial-algebrice se clasific n funce de diveri parametrii:

-numrul de condiii de tip algebric: din ecuaia diferenial-algebric; an

-numrul de condiii de tip diferenial: din ecuaia diferenial-algebric; dn

-indexul diferenial: - numrul de operaii de difereniere necesare convertirii ecuaei diferenial-algebrice ntr-un sistem de ecuaii difereniale.

di

Exist o mare diferen ntre ecuaiile diferenial-algebrice i ecuaiile difereniale n privina condiiilor la iniiale ataate pentru asigurarea unicitii soluiei: astfel, pentru o ecuaie diferenial tim precis cte condiii iniiale sunt necesare, pe cnd pentru o ecuaie diferenial-algebric situaia devine neclar. De exemplu, pentru ecuaia (trivial), privind-o doar ca ecuaie neliniar nu necesit nici o condiie iniial asociat. 0=ye

Pentru rezolvarea numeric a ecuaiei diferenial-algebric

0))('),(,( =xyxyxF presupunnd cunoscute valorile aproximative ale soluiei exacte )(xyy = n punctele , pentru determinarea valorii aproximative n punctul , utiliznd o metod numeric de rezolvarea a ecuaiilor

difereniale, vom nlocui printr-o relaie de forma .

knii xxx ...,, ,11+ix

)(' 1+ixy =

+k

ii xyc0

1 )(

Astfel, ecuaia iniial, pentru devine: 1+= ixx

=

+++ =k

iiii ycyxF0

111 0),,(

adic o ecuaie neliniar n raport cu . Rezolvarea acesteia se va efectua prin intermediul metodei Newton sau metodei aproximaiilor succesive, sau a oricrei metode specifice acestui tip de ecuaii.

)( 11 ++ ii xyy

5. METODA DIFERENTELOR FINITE PENTRU PROBLEMA STURM-LIOUVIILLE

Considerm problema bilocal (Sturm-Liouville) dat de:

(P)

==

++=

)()(

],[ ),()()()(')()(''

byay

baxxrxyxqxyxpxy

Teorem: Dac: (i). sunt continue pe ; )( ),( ),( xrxqxp ],[ ba (ii). ],[)( ,0)( baxxq >atunci problema (P) admitt o soluie unic.

Fie o diviziune echidistant de pas h a lui [a, b], x0 = a, xi = a + ih, 0 i n cu xn = b.

Teorem a) Fie y C3 [a, b], x (a, b), h > 0 astfel nct x h, x + h (a, b). Atunci (x h, x + h) astfel nct: )(

h2)hx(y)hx(y + = )('''y

6h)x('y

2

+ b) Fie y C4 [a, b], x (a, b), h > 0 astfel nct

x h, x + h (a, b). Atunci (x h, x + h) astfel nct: )(

2h)hx(y)x(y2)hx(y ++ = )(y

12h)x(''y IV

2

+ Demonstraie. Se folosete formula dezvoltrii n serie Taylor i Teorema de medie

Pentru i = 1n,1 facem aproximrile:

h2)x(y)y(x )x('y 1i1ii +

= (*)

21ii1i

i h)x(y)x(y2)x(y)x(''y + += (**)

Dac n ecuaia (P) considerm x = xi, 1 i n 1 i folosim (*) i (**) avem

211 )()(2)(

hxyxyxy iii + + =

hxyxyxp iii 2

)()()( 11 + )()()( iii xrxyxq ++

iar condiiile (25) devin

==

)x(y)x(y

n

0

Obinem n final schema:

(S) 211 2

hyyy iii + + -

hyyp iii 2

11 + - 11 , = niryq iii == nyy ,0

unde )( ),( ),( ),( iiiiiiii xyyxrrxqqxpp ==== i)( , care reprezint un sistem liniar de n-1 ecuaii i n-1 necunoscute y1, y2, ..., yn-1, matricea sistemului fiind tridiagonal dominant diagonal.

Exemplu: Fie problema

==

++=

2)1(1)0(

]1,0[,44 3

yy

xeyyy xIII

S se aproximeze valorile soluiei n punctele , hixi = 31 i , 4/1=h , cu ajutorul metodei

diferenelor finite. Soluie: Schema (S) ne conduce la sistemul urmtor:

==

==

+

+

+ +

21

3,1 ,214221

4

0

312212

yy

ieyhh

yh

yhh

ixiii

Folosind o metod numeric, pentru 4/1=h rezult:

===

5910,01754,06468,0

3

2

1

yyy

BIBLIOGRAFIE i Webografie: 1. Berbente C., Mitran S., Zancu, S., Metode Numerice, Editura Tehnic, Bucureti, 1997 2. R. Burden, J. Faires, Numerical Analysis, PWS-Kent, 2005. 3. C. Carasso, Analyse Numrique, Lidec, Canada, 1970. 4. P.Ciarlet, J.Lions, Finite Difference Methods, North-Holland, Amsterdam, 1989. 5. B. Demidovitch. I. Maron, lments de Calcul Numrique, Mir, Moscow, 1973. 6. D. Ebnc, Metode de Calcul Numeric, Editura Sitech, Craiova, 1994. 7. Kinkaid, D., Cheney, W., Numerical Analysis Mathematics of Scientific Computing,

American Mathematical Society, 2009. 8. R. Militaru Mthodes Numriques. Thorie et Applications Ed. Sitech, 2008 9. J.P. Nougier, Mthodes de Calcul Numrique, Hermes Sciences Publication, Paris, 2001. 10. M. Popa, R. Militaru, Metode Numerice n pseudocodaplicaii, Ed. Sitech, Craiova, 2010. 11. M. Popa, R. Militaru, Analiz numeric note de curs, Editura Sitech, Craiova, 2009. http://www.ac-nancy-metz.fr/enseign/physique/divers/MethodNum/method-num.htm http://www.ac-nancy-metz.fr/enseign/physique/divers/MethodNum/Schw- doc/EULER03.pdf

http://ta.twi.tudelft.nl/nw/users/vuik/wi211/disasters.html

Exemplul 1RezolvareExemplul 2Exemplul 3Rezolvare.2. Metode de tip Runge-KuttaCazuri particulare2.2. Metod de tip Runge-Kutta implicit4.2 Metoda Adams-Bashforth-Moulton de ordinul patru