REZOLVAREA NUMERICĂ A ECUAŢIILOR DIFERENŢIALE Introducere Acest tip de probleme provine din cadrul vast al analizei funcţionale. Ecuaţiile diferenţiale sau cu derivate
parţiale constituie modelele matematice pentru majoritatea problemelor inginereşti: studiul eforturilor la care sunt supuse elementele de rezistenţă: bare, grinzi, plăci subţiri, groase, conducte; studiul problemelor de câmp electric în dielectrici, câmp magnetic, câmp termic, propagarea undelor, curgerea fluidelor etc.
Odată stabilit fenomenul fizico-tehnic şi ecuaţiile diferenţiale care îl guvernează, ca formă, coeficienţi, condiţii la limită (pe frontieră) rămâne de rezolvat ultima problemă: rezolvarea acestui model matematic. Din diverse motive: neomogenităţile fizice, frontiere cu geometrie dificilă, număr de necunoscute, etc., rezolvarea o vom face căutând o soluţie aproximativă cu ajutorul unui cod numeric, folosind calculatorul.
O ecuaţie diferenţială este o ecuaţie în care necunoscuta este o funcţîe şi în care intervine funcţia
necunoscută, derivatele ei de diverse ordine şi variabile independente de care depind aceste funcţii. In cazul în care funcţia necunoscută depinde de o singură variabilă independentă, ecuaţia se numeşte
ecuaţie diferenţială ordinară, iar în situaţia în care funcţia necunoscută depinde de mai multe variabile independente, ecuaţia se numeşte cu derivate parţiale.
Ordinul unei ecuaţii diferenţiale este cel mai înalt ordin a derivatei funcţiei necunoscute ce figurează în
ecuaţia respectivă. Expresia generală a unei ecuaţii diferenţiale, sub formă implicită este:
0,...,,,, 2
2
=⎟⎟⎠
⎞⎜⎜⎝
⎛n
n
dxyd
dxyd
dxdyyxF
unde . )(xyy = Forma explicită se poate scrie:
⎟⎟⎠
⎞⎜⎜⎝
⎛= −
−
1
1
2
2
,...,,,, n
n
n
n
dxyd
dxyd
dxdyyxf
dxyd
Privitor la condiţiile la limită există două tipuri: (i). Condiţii Cauchy: se cunosc într-un punct atât valoarea funcţiei necunoscute cât şi valorile derivatelor, până la ordinul cel mai mare ce figurează în ecuaţie;
0x
(ii). Condiţii la limită: se cunosc valorile funcţiei necunoscută în puncte diferite. Rezolvarea numerică a unei probleme asociate unei ecuaţii diferenţiale poate fi privită sub două aspecte:
(a). determinarea unei funcţii )(~~ xyy = , aparţinând unei anumite clase de funcţii (în general polinoame, dată fiind importanţa lor teoretică fundamentală), şi care aproximează „suficient de bine” soluţia exactă )(xyy = ; (b). determinarea valorilor aproximative ale soluţiei exacte )(xyy = , într-o mulţime de puncte date
. ,...1,0 , =ixi
Vom expune în continuare principalele metode numerice al căror algoritm are un cost de calcul redus şi se
pretează la implementarea pe calculator pentru rezolvarea numerică a ecuaţiilor diferenţiale. Pentru ecuaţii diferenţiale ordinare acestea se pot clasifica în două mari tipuri: - metode unipas (de tip Euler, Runge-Kutta) în care determinarea valorii aproximative a soluţiei în fiecare
punct se va obţine direct, pe baza informaţiilor din punctul precedent;
- metode multipas, în care valoarea soluţiei exacte în fiecare punct este aproximată folosind informaţiile din mai multe puncte anterioare.
Evident este vorba de soluţii aproximative pe care nu avem cum să le comparăm cu o soluţie exactă, deoarece practic aceasta este imposibil de găsit.
De aceea în practică trebuie să procedăm cu atenţie pentru alegerea algoritmilor cei mai potriviţi pentru problema concretă de rezolvat.
1. METODE DE TIP EULER
1.1. Metoda Euler
Se consideră ecuaţia diferenţială:
y′ = ƒ(x, y) (1) cu condiţia iniţială:
y(x0) = y0, (2) unde funcţia ƒ este definită într-un domeniu D din planul xOy. Perechea (1),(2) constituie o problemă Cauchy. Presupunem asigurate existenţa şi uniciitatea soluţiei. Se defineşte un câmp de direcţii în D dacă în fiecare punct M(x, y) ∈ D se ia direcţia α = arctg ƒ(x, y) (α fiind unghiul format de direcţie cu sensul pozitiv al axei Ox).
Acest câmp de direcţii are următoarea interpretare: graficul soluţiei ecuaţiei (1) cu condiţia (2) trece prin punctul M(x0, y0) şi este tangent în orice punct al său direcţiilor câmpului.
Metoda lui Euler propune aproximarea soluţiei printr-o linie poligonală în care fiecare segment este coliniar cu direcţia câmpului definită de extremitatea sa stângă. Astfel se consideră nodurile echidistante xi = x0 + ih, n0i ,= . În punctul M0(x0, y0) se calculează direcţia câmpului definită de M0 şi se scrie ecuaţia dreptei determinate de M0 şi de această direcţie:
y = y0 + ƒ(x0, y0)(x – x0) (3) Funcţia (3) se propune ca aproximantă a soluţie problemei (1)+(2) pe [x0,x1]. Valoarea aproximativă a
soluţiei în x1 este dată de: y1 = y0 + ƒ(x0, y0)(x1 – x0) = y0 + ƒoh Repetăm procedeul şi presupunând că în xi s-a calculat valoarea aproximativă yi, atunci pe intervalul [xi,
xi+1] se aproximează soluţia cu: y = yi + ƒ(xi, yi)(x – xi) = yi + ƒi(x – xi),
iar în punctul xi+1 se obţine valoarea aproximativă: yi+1 = yi + hƒi.
Aproximarea este justificată şi de următoarea 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ξ−=−−
Demonstraţie. Din formula lui Taylor avem:
y(x + h) = y(x) + )(''y2
h)x('hy 1
2
ξ+ , ∈ξ1 (x, x + h)
y(x – h) = y(x) – )(''y2
h)x('hy 2
2
ξ+ , ∈ξ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ε+
−−
Considerând cunoscută aproximarea yi a soluţiei problemei (1)+(2) în xi, procedeul de aproximare Euler, poate fi acum rezumat astfel:
(4) ,...1,0),(
1
1 =⎪⎩
⎪⎨
⎧
+=+=
=
+
+ ihfyyhxxyxff
iii
ii
iii
Observaţii: (i).Neglijarea termenilor de ordin superior în (4) face ca metoda să fie comodă în calcul, dar puţin precisă, erorile cumulându-se la fiecare pas. (ii).Metoda se poate aplica şi dacă nodurile xi nu sunt echidistante, având la fiecare iteraţie alt pas h în acest caz. In general o metodă unipas poate fi scrisă sub forma:
);,(1 hyxhyy iiii Γ+=+
Pentru soluţie a problemei Cauchy (1)+(2) (şi pentru orice )(xyy = }1,...,2,1,0{ −∈ ni ) considerăm:
( ) ));(,()()(1)( 1 hxyxxyxyh
h iiiii Γ−−= +τ
Cantitatea de mai sus se numeşte eroare de consistenţă a metodei în şi reprezintă o măsură a calităţii metodei de aproximare.
ix
Pentru metoda Euler, ţinând cont de algoritmul acesteia, de expresia erorii de consistenţă şi de faptul că
obţinem 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ă determinăm o soluţie 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, obţinem:
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,2·0,8667 =
= 1,3733. Obţinem în final următorul tabel, ultima coloană reprezentând valorile exacte ale soluţiei 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 câteva variante îmbunătăţite ale algoritmului Euler. 1.2. Metoda Euler modificată Considerăm pe [xi, xi+1] ca direcţie a segmentului MiMi+1 direcţia definită de punctul de la mijlocul
segmentului (nu de extremitatea stângă ca în formula iniţială) se obţine metoda Euler modificată. Dacă xi, yi sunt valori calculate, procesul iterativ este următorul:
⎪⎪
⎩
⎪⎪
⎨
⎧
+=⎟⎟⎠
⎞⎜⎜⎝
⎛=+=
+==+=
++
++++
+
;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 aceleiaşi probleme de la exemplul 1 pentru a putea face
o comparaţie rapidă privind eficienţa 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 +
2h 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
Continuând în acelaşi mod obţinem 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 Euler–Cauchy Introducând y′ ca o nouă variabilă, se poate înlocui ecuaţia diferenţială (1) cu următorul sistem de ecuaţii:
⎪⎩
⎪⎨⎧
=
=
).y,x(f'y
'ydxdy
Integrând prima ecuaţie a sistemului de la xi la xi+1, unde xi = x0 + ih se obţine:
⎪⎪⎩
⎪⎪⎨
⎧
=
+= ∫+
+
)y,x(fy
;dx'yyy
ii'i
1x
xi1i
i
i
Folosind metoda dreptunghiului pentru aproximarea integralei din prima ecuaţie obţinem: yi+1 = yi + hy′i = yi + hƒ(xi, yi). (5) Dacă aplicăm metoda trapezului pentru aproximarea aceleiaşi integrale se obţine:
yi+1 = yi + )(2
'1
'++ ii yyh = yi + ))y,x(f)y,x(f(
2h
1i1iii +++ =
= yi + )(2 1++ ii ffh . (6)
Ecuaţia (6) poate fi rezolvată iterativ, relativ la necunoscuta yi+1. Dacă se alege ca aproximaţie iniţială dată de (5), obţinem pentru yi+1 valoarea:
)0(1+iy
)),(),((2
)0(11
)1(1 +++ ++= iiiiii yxfyxfhyy . (7)
Deci pentru construirea liniei poligonale care aproximează soluţia problemei (1) + (2) este dat algoritmul Euler–Cauchy:
,...1,0
).(2
);,(
;
;
)0(1
)1(11
)0(11
)0(1
)0(1
1
=
⎪⎪⎪
⎩
⎪⎪⎪
⎨
⎧
++==
=
+=
+=
+++
+++
+
+
i
ffhyyy
yxff
hfyy
hxx
iiiii
iii
iii
ii
Pentru această metodă )(hiτ este de ordinul . )( 2hO
Observaţie: In scopul ameliorării aproximaţiilor se pot utiliza metodele prezentate de o manieră iterativă. Astfel, metode euler-Cauchy ne conduce la următorul algoritm de calcul:
;)0(1 iii hfyy +=+
)(2
)1(1
)(1
−++ ++= k
iiik
i ffhyy , k = 1, 2, …
cu ( ))1k(1i1i
)1k(1i y,xff −
++−
+ = Se demonstrează că dacă funcţia este lipsichitziană în y, de constantă L şi dacă h este suficient de mic cu
hL < 2 atunci , pentru k . 1)(1 ++ → i
ki yy ∞→
În practică, se continuă iteraţiile după indicele k până când un criteriu de eroare impus de exemplu:
)k(1i
)1k(1i
)k(1i
yyy
+
−++ − < ε, ε impus, 0, este satisfăcut. ≠+
)k(1iy
Exemplul 3 Să se găsească soluţia 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,xf
205,0yy +++ =
= 1 + 205,0 (1 + 1,1025) = 1,0525625
Deci fiind stabile primele 4 zecimale rezultă: y1 = 1,0525. Similar cu x2 = 0,1 obţinem:
)y,x(hfyy 111)0(
2 += =1,1077 )1(
2y = 1,11036; = 1,11042. )2(2y
Deci y2 = 1,1104. Soluţia exactă fiind y = 2ex – x – 1, obţinem în x = 0,1 valoarea: y(0,1) = 1,1103 eroarea fiind de 0,0001.
1.4. Metoda Euler–Heun
Vom încheia trecerea în revistă a variantelor Euler pentru rezolvarea ecuaţiilor de tipul (1) prezentând o
ultimă variantă de ordinul doi a metodei lui Euler şi anume aceea propusă de Heun, [26]. Presupunând calculată valoarea yi la pasul xi, se propune pentru calcularea soluţiei la pasul xi+1 expresia:
41hyy ii +=+ ⎥
⎦
⎤⎢⎣
⎡⎟⎠⎞
⎜⎝⎛ +++ iiiiii yxhfyhxfyxf ,(
32,
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 scăzut
al erorii de consistenţă conduce la o aplicabilitate limitată. In scopul obţinerii unor metode de ordin ridicat trebuie renunţat fie la proprietatea de a fi unipas şi pastrată liniaritatea, fie viceversa. Metodele de tip Runge-Kutta sunt neliniare şi conservă caracteristicile metodelor unipas, având un ordin ridicat. Ele au trei proprietăţi principale:
1. Sunt metode directe, adică pentru determinarea aproximării soluţiei la pasul i+1 avem nevoie de
informaţiile existente în punctul precedent xi, yi. 2. Sunt identice cu seriile Taylor până la termenii hn, unde h este pasul curent iar n este diferit pentru
metode diferite din această familie şi defineşte ordinul metodei. 3. în procesul de calcul nu necesită decât evaluarea funcţiei 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 întâi iar metodele Euler-Cauchy şi Euler-Heun sunt metode R-K de ordinul 2.
2.1. Construirea formulelor Runge-Kutta Ne ocupăm în special de rezolvarea ecuaţiilor diferenţiale ordinare cu condiţii iniţiale. Adică pentru
ecuaţii diferenţiale ordinare de ordinul întâi 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 soluţiei 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 – urmând a fi determinaţi (i = m,2 ; j = 1i,1 − )
Observaţie: Pentru a fi consisente, metodele Runge-Kutta trenuie să satisfacă condiţia: . ∑=
=m
iic
1
1
Notăm ε(h) = y(x+h) – z(h) (eroarea de aproximare). Vom determina parametrii c1, ci, ai, bij din condiţiile:
⎢⎢⎣
⎡
≠ε
=ε==ε=ε+ f funcţie anumită opentru ,0)0(
şi f funcţie 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 mărime 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ă uşor că ε(0) = 0 pentru orice funcţie f. ε’(h) = y’(x+h) – c1f(x,y) Obţinem ε’(0) = y’(x) – c1f(x,y) = (1 – c1)f(x,y). Se observă că pentru c1 = 1, ε’(0) = 0 pentru orice funcţie 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ă funcţie f. De exemplu, pentru f(x,y) = y se obţine ε’’(0) = y. Deci p = 1.
Astfel, pentru c1 = 1, se obţine 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ă uşor că ε(0) = 0 pentru orice funcţie 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
+++ ]
Notăm cu: expresia inclusă în parantezele drepte. /hf )kby,hax( 1212 ++
Obţinem: ε’(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 ′′Obţinem: ε’’(0) = y’’(x) – 2c2[a2fx (x,y) + b21f(x,y) · fy(x,y)] =
=(fx + f·fy) – 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 ++′′′
Obţinem: ε’’’(0) = y’’’(x)–3c2( fxx+a2b21f·fxy+a2b21f·fyx+ f2·fyy) = 22a 2
21b= (fxx + f·fxy + f·fyx + fx·fy + f·f2
y + f2·fyy) – 3c2( fxx + a2b21f·fxy+ a2b21f·fyx + f2·fyy) = (1 – 3c2 ) fxx + (1 – 3c2a2b21) f·fxy+
22a 2
21b 22a
+ (1 – 3c2a2b21)f·fyx + (1 – 3c2 ) f2 · fyy + fx·fy + f ·221b 2
yf Dacă:
⎪⎩
⎪⎨
⎧
=−=−=−−
0bc210ac210cc1
212
22
21
sau ⎪⎩
⎪⎨
⎧
==
=+
2212
21
c21ba
1cc (10)
atunci ε(0) = ε’(0) = ε’’(0) = 0 pentru orice funcţie f. Sistemul (10) este compatibil nedeterminat. De asemenea, ε’’’(0) nu este identic nulă pentru orice funcţie f. De exemplu, pentru funcţia f(x,y) = y, se
obţine ε’’’(0) = y. Deci p= 2. Pentru m = p = 2 se pot obţine oricâte formule dorim, alegând parametrii c1, c2, a2, b21 astfel încât să
verifice (10). Din (9) rezultă că eroarea în aceste formule este de ordinul lui h2.
Astfel, o soluţie a sistemului (10) este: c1 = c2 = 21 , a2 = b21= 1, şi pentru această soluţie se obţine formula
Runge-Kutta simplă:
⎪⎪
⎩
⎪⎪
⎨
⎧
++==
++=+
)ky,hx(hfk)y,x(hfk
)kk(21)x(y)hx(y
12
1
21
(10’)
Observaţie. Dacă f nu depinde de y obţinem:
[ ])hx(f)x(f2h)x(y)hx(y ++=−+
Pe de altă parte:
∫ ∫+ +
==−+hx
x
hx
x
dt)t(fdt)t('y)x(y)hx(y
Din cele două relaţii de mai sus rezultă formula trapezului (vezi § 5.1.).
[ ]∫+
++=hx
x
)hx(f)x(f2hdt)t(f
O altă formulă Runge-Kutta de ordinul 2 se obţine pentru
c1 = 41 , c2 =
43 , 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 soluţiei unei probleme Cauchy (8) pe un interval [a, b] de noduri
echidistante de pas h procedăm astfel : - presupunem cunoscute valorile xi, yi (yi ≅ y(xi)) – de exemplu, pentru i = 0 cunoaştem din (8) x0, y0; - calculăm, 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 soluţia problemei date în xi = 0.1·i, 1 ≤ i 3. ≤
Rezolvare f(x,y) = xy; x0 = 0; y0 = 1; h = 0.1 Se observă că soluţia exactă este y(x) = . Deci 2/x2
e
y(0,1) 1,005013 ≅y(0,2) 1,020201 ≅y(0,3) 1,046027 ≅Vom întocmi următorul 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) Efectuând calcule asemănătoare 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 soluţii. O soluţie a acestui sistem este:
c1 = 61 ; c2 =
32 ; c3 =
61 ; a2 =
21 , a3 = 1
b21 = 21 , b31 = -1, b32 = 2
obţinându-se următoarea 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 următoarea:
⎪⎪⎪⎪
⎩
⎪⎪⎪⎪
⎨
⎧
⎟⎠⎞
⎜⎝⎛ ++=
⎟⎠⎞
⎜⎝⎛ ++=
=
++=+
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 funcţie f. Într-adevăr pentru f(x,y) = y se obţine εIV(0) = y.
Observaţii. (i). Din (9) rezultă că eroarea în aceste formule este de ordinul lui h3. (ii).Dacă funcţia f nu depinde de y, atunci
∫+hx
x
dt)t(f = ∫+hx
x
dt)t('y =y(x+h)–y(x) = ( )⎥⎦
⎤⎢⎣
⎡++⎟
⎠⎞
⎜⎝⎛ ++ hxf
2hxf4)x(f
6h
adică formula de integrare numerică Simpson. 4) m = 4. Raţionamente şi calcule asemănătoare celor anterioare conduc la construirea unei mulţimi 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 Observaţie. Aceste formule sunt de ordinul patru. .
Pentru calculul aproximativ al soluţiei unei probleme Cauchy (8), pe un interval [a, b], cu pasul h procedăm astfel:
x0 = a, xi = x0 + ih, i = N,0 cu xN = b , deci h = N
ab −
cunoscând pentru xi, valorile yi (yi y(xi)) – de exemplu pentru i = 0, cunoaştem x0, y0. ≈Calculăm: xi+1 = xi + h
yi+1 = yi + 61 (k1 +2k2 + 2k3 + k4)
k1 = hf(xi, yi)
k2 = hf ⎟⎠⎞
⎜⎝⎛ ++
2ky,
2hx 1
ii
k3 = hf ⎟⎠⎞
⎜⎝⎛ ++
2ky,
2hx 2
ii
k4 = hf(xi + h, yi +k3) Exemplul 5
Se consideră problema Cauchy: . ⎩⎨⎧
==
1)0(yxy'y
Să aproximăm soluţia acestei ecuaţii în punctele xi = ih, i = 9,1 , h = 0.1, folosind o metodă Runge-Kutta de ordin 4. xi yi valoarea soluţiei 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/x2
e
2.2. Metodă de tip Runge-Kutta implicit Metodele de tip Runge-Kutta, expuse anterior, sunt explicite. Pentru a ameliora proprietăţile de
stabilitate ale acestor metode se consideră cele de tip implicit.
Prin proprietăţi de stabilitate ne referim la restricţiile impuse asupra pasului de integrare în situaţia
utilizării 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, funcţiile kj numai sunt definite explicit, ci printr-un set de m
ecuaţii implicite, în general neliniare. În practică, se foloseşte cazul metodelor Runge-Kutta implicite cu
= 2
m
m
Astfel:
kj = f(x + haj, y + bj1hk1 + bj2hk2), j = 2,1
Considerând k1, k2 dezvoltabile sub forma
kj = Aj + hBj + h2Cj + h3Dj + O(h4), j = 2,1
pe baza dezvoltării în serie Taylor în raport cu (x, y) a lui kj, j = 2,1 se obţine prin identificarea puterilor lui
h:
⎪⎪⎩
⎪⎪⎨
⎧
==
±=
;41bb
;63
21a
2211
1
41ab;
41ab
63
21a
221112
2
−=−=
= m (16)
Observaţie. Prin simetrie se observă faptul că valorile date prin alternarea semnelor, în relaţiile
(16), conduc la aceeaşi 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 existenţei unei metode Runge-Kutta de tip semi-explicit de ordinul
patru, descrisă de relaţiile 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 ≤<<<≤ ...10 . Forma generală a unei metode multipas este dată printr-o relaţie de recurenţă ce exprimă valoarea în
raport cu valorile lui şi în şi în punctele precedente: 1+iy
)(xy )(' xy 1+ix
(*
110111
101211
N ,...,,1
),(...),(),( ...
∈−=
++++ )++++=
−+−+−++
−+−−−+
mmmi
yxfbyxfbyxfbhyayayay
mimiiimiim
miimimi
(20) Observatie: Dacă numim metoda explicită; altfel o numim implicită. 0=mb Se numeşte eroare de consistenţă a metodei în , cantitatea: ix
( )
))(,(..))(,())(,(
)(...)()(1)(
110111
10111
mimiiimiim
miimii
xyxfbxyxfbxyxfb
xyaxyaxyh
h
−+−+−++
−+−−+
−−−−
−−−−=τ
(21)
Printre metodele explicite multipas, cele mai utilizate sunt cele de tip Adams-Bashforth: (i). de ordinul trei, cu trei paşi:
,...3,2
)),(5),(16),(23(12 22111
=
+−+= −−−−+
i
yxfyxfyxfhyy iiiiiiii
(22) (ii). de ordinul patru, cu patru paşi:
,...4,3
)),(9),(37),(59),(55(24 3322111
=
−+−+= −−−−−−+
i
yxfyxfyxfyxfhyy iiiiiiiiii
(23)
Cele mai cunoscute metode implicite, multipas, sunt cele de tip Adams-Moulton: (i). de ordinul trei, cu doi paşi:
,...2,1
)),(),(8),(5(12 11111
=
−++= −−+++
i
yxfyxfyxfhyy iiiiiiii
(24) (ii). de ordinul patru, cu trei paşi:
,...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 combinaţie între o metodă numerică explicită şi o
metodă numerică implicită. Metoda explicită permite predicţia unei valori aproximative şi cea implicită corectează (odată sau de mai multe ori) această predicţie.
Metodele predictor-corector, oferă o precizie superioară faţă de metodele prezentate în paragrafele anterioare, fără a implica sporiri ale numărului de operaţii aritmetice.
În ceea ce priveşte 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 kiii
ki
Criteriul de oprire este:
ε≤− −++
)1(1
)(1
ki
ki yy sau ε≤
−
+
−++
)(1
)1(1
)(1
ki
ki
ki
y
yy cu şi 0)(
1 ≠+k
iy ε precizia de calcul impusă.
4.1 Metoda lui Milne
Algoritmul de calcul este următorul:
( )
( )⎪⎪⎪
⎩
⎪⎪⎪
⎨
⎧
+++=
+−+=
+=
+−−+
−−−+
+
)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 ameliorării rezultatelor putem aplica formulele de calcul de manieră iterativă:
( ) ,...2,1 ,43
)1(111
)(1 =+++= −
+−−+ kfffhyy ik
iiik
i
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 utilizăm o metodă de ordinul patru. iiii yyyy ,, 1,23 −−−
Corectăm 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 Numérique, Lidec, Canada, 1970. 3. B. Demidovitch. I. Maron, Éléments de Calcul Numérique, Mir, Moscow, 1973. 4. D. Ebâncă, Metode de Calcul Numeric, Editura Sitech, Craiova, 1994. 5. R. Militaru Méthodes Numériques. Théorie et Applications Ed. Sitech, 2008 6. J.P. Nougier, Méthodes de Calcul Numérique, Hermes Sciences Publication, Paris, 2001. 7. M. Popa, R. Militaru, Metode Numerice –algoritmi si aplica�ii, 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ă ecuaţie diferenţială, anterior prezentate, se pot extinde şi în cazul sistemelor de ecuaţii diferenţiale. Considerăm situaţia metodelor numerice de tip Runge-Kutta. Metodele Runge-Kutta, prezentate pentru ecuaţii diferenţiale, pot fi aplicate cu uşurinţă şi la sisteme de ecuaţii diferenţiale. Fie sistemul de ecuaţii diferenţiale:
(1) ⎪⎩
⎪⎨
⎧
=′
=′
)y...,y,x(fy................................),y,...,y,x(fy
n1nn
n111
şi urmărim să determinăm soluţia care satisface condiţiile iniţiale: yi(x0) = yi,0 , i = n,1 (2)
Presupunând că dispunem de soluţia problemei (1) + (2) la pasul i: y1,i , ..., yn,i metoda Runge-Kutta de ordinul 4 calculează soluţia în pasul i+1 cu formulele:
y1,i+1 = y1,i + Δy1,i ...........................
(3) yn,i+1 = yn,i + Δyn,i Corecţiile: Δy1,i , Δy2,i , ... , Δyn,i care intervin în formulele precedente se calculează cu ajutorul relaţiilor:
Δy1,i = 14
13
12
11 k
61k
31k
31k
61
+++ ;
....................................................
Δyn,i = n4
n3
n2
n1 k
61k
31k
31k
61
+++ ;
iar coeficienţii au următoarea formă: n4
n1
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,n
13i,1ij
j4 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 Considerăm 10 noduri echidistante xi pe [0, 1], de pas h = 0.1. Din (14), pentru fiecare i = 0, 1, …, 9, avem :
yi+1 = yi + ( )14
13
12
11 kk2k2k
61
+++
zi+1 = zi + ( )24
23
22
21 kk2k2k
61
+++
unde: )z,y,x(hfk iiij
j1 = ; 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 23i
13iij
j4 +++= ; j = 1, 2
Obţinem următoarele 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 ⇒
soluţia exactă fiind , deci ⎩⎨⎧
=
=2x)x(z
x)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 ecuaţia diferenţială de ordinul scrisă sub formă explicită n
(4) ))(),...,('),(,()( )1()( xyxyxyxfxy nn −=
cu condiţiile asociate
(5)
⎪⎪⎩
⎪⎪⎨
⎧
=
==
−−
10)1(
10
00
)(..................
)(')(
nn txy
txytxy
Introducem notaţiile
(6) n1,2,..,j ),()( )1( == − xyxz jj
Remarcă: 1,..,2,1 ),()( 1
' −== + njxzxz jj
Tinând cont de (4), (5) şi (6) avem:
(7)
⎪⎪⎪
⎩
⎪⎪⎪
⎨
⎧
=
=
=
=
−
))(),...,(,()(
)()(......................)()(
)()(
1'
'1
3'2
2'1
xzxzxfxz
xzxz
xzxz
xzxz
nn
nn
respectiv:
(8)
⎪⎪⎩
⎪⎪⎨
⎧
=
==
−10
102
001
)(..................)()(
nn txz
txztxz
In concluzie (7) şi (8) reprezintă un sistem având n ecuaţii diferenţiale ordinare cu n condiţii iniţiale. Putem astfel să îl rezolvăm 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 Transformăm ecuaţia diferenţială de ordinul doi în sistemul de ordinul întâi:
⎩⎨⎧
−=−+=====
1)0(z ;zyx)z,y,x(f'z1)0(y ;z)z,y,x(f'y
2
1
Considerăm 10 noduri echidistante xi pe [0,1] de pas h = 0.1. Din (14), pentru fiecare i = 0, 1, …, 9 avem:
( )14
13
12
11i1i kk2k2k
61yy ++++=+
( 24
23
22
21i1i kk2k2k
61zz ++++=+ ) - 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 23i
13iij
j4 +++= ; j = 1,2
Obţinem următoarele 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
făcând apel la o manieră eficientă de alegere a mulţimii de puncte în care se va aproxima soluţia exactă . )(xyIdeea de bază este utilizarea unor metode numerice de ordine diferite în scopul caracterizării erorii de consistenţă şî alegerii pasului de integrare astfel încât eroarea globală să fie inferioară unei anumite precizii ε impuse. Considerăm în continuare două metode numerice care conduc la aproximaţiile următoare ale soluţiei 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 obţinute cu acelaşi pas . h
Atunci:
)()),(,()()( ),,()()(
11
111
hhhxyxhxyxyhyxhyxyyxy
iiiii
iiiiii
++
+++
=Φ−−≅≅Φ−−=−
τ
Tinând cont că şi că )()(1
ni hOh =+τ )()(~ 1
1+
+ = ni hOhτ rezultă
( 111~1)( +++ −≅ iii yy
hhτ ) (c)
de unde (d) n
i khh ≅+ )(1τk constantă ce nu depinde de . h
In consecinţă pentru un nou pas de integrare : qh
)()( 11 hqqh in
i ++ ≅ ττ
Impunând ετ ≤+ )(1 qhi obţinem: 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 estimării 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 funcţia necunoscută este o funcţie scalară sau vectorială. )(xyy =
Vom aborda în continuare cazul ecuaţiilor diferenţial-algebrice care sub formă explicită devin:
)()()()(')( xfxyxAxyxM +=
Exemplu:
Dându-se ecuaţia diferenţial-algebrică , aceasta se poate explicita sub forma
, unde , ,
⎪⎩
⎪⎨
⎧
=++−=
+−=
2424'
467'
wvuvuv
xvuu
)()()()(')( xfxyxAxyxM +=⎟⎟⎟
⎠
⎞
⎜⎜⎜
⎝
⎛=
000010001
M⎟⎟⎟
⎠
⎞
⎜⎜⎜
⎝
⎛−−
=111924067
A⎟⎟⎟
⎠
⎞
⎜⎜⎜
⎝
⎛
−=
240
4xf
Ecuaţiile diferenţial-algebrice se clasifică în funcţîe de diverşi parametrii:
-numărul de condiţii de tip algebric: din ecuaţia diferenţial-algebrică; an
-numărul de condiţii de tip diferenţial: din ecuaţia diferenţial-algebrică; dn
-indexul diferenţial: - numărul de operaţii de diferenţiere necesare convertirii ecuaţîei diferenţial-algebrice într-un sistem de ecuaţii diferenţiale.
di
Există o mare diferenţă între ecuaţiile diferenţial-algebrice şi ecuaţiile diferenţiale în privinţa condiţiilor la iniţiale ataşate pentru asigurarea unicităţii soluţiei: astfel, pentru o ecuaţie diferenţială ştim precis câte condiţii iniţiale sunt necesare, pe când pentru o ecuaţie diferenţial-algebrică situaţia devine neclară. De exemplu, pentru ecuaţia (trivială), privind-o doar ca ecuaţie neliniară nu necesită nici o condiţie iniţială asociată. 0=ye
Pentru rezolvarea numerică a ecuaţiei diferenţial-algebrică
0))('),(,( =xyxyxF
presupunând cunoscute valorile aproximative ale soluţiei exacte )(xyy = în punctele , pentru
determinarea valorii aproximative în punctul , utilizând o metodă numerică de rezolvarea a ecuaţiilor
diferenţiale, vom înlocui printr-o relaţie de forma .
knii xxx −− ...,, ,1
1+ix
)(' 1+ixy ∑=
−+
k
ii xyc0
1 )(α
α
Astfel, ecuaţia iniţială, pentru devine: 1+= ixx
∑=
−+++ =k
iiii ycyxF0
111 0),,(α
α
adică o ecuaţie neliniară în raport cu . Rezolvarea acesteia se va efectua prin intermediul metodei Newton sau metodei aproximaţiilor succesive, sau a oricărei metode specifice acestui tip de ecuaţii.
)( 11 ++ ≅ ii xyy
5. METODA DIFERENTELOR FINITE PENTRU PROBLEMA STURM-LIOUVIILLE
Considerăm 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 soluţie 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 încât x – h, x + h ∈ (a, b). Atunci (x – h, x + h) astfel încât: ∈ξ∃)(
h2)hx(y)hx(y −−+ = )('''y
6h)x('y
2
ξ+
b) Fie y ∈ C4 [a, b], x ∈ (a, b), h > 0 astfel încât x – h, x + h ∈ (a, b). Atunci (x – h, x + h) astfel încât: ∈ξ∃)(
2h)hx(y)x(y2)hx(y −+−+ = )(y
12h)x(''y IV
2
ξ+
Demonstraţie. Se foloseşte formula dezvoltării în serie Taylor şi Teorema de medie
Pentru i = 1n,1 − facem aproximările:
h2)x(y)y(x )x('y 1i1i
i−+ −
= (*)
21ii1i
i h)x(y)x(y2)x(y)x(''y −+ +−
= (**)
Dacă în ecuaţia (P) considerăm x = xi, 1 ≤ i ≤ n – 1 şi folosim (*) şi (**) avem
211 )()(2)(
hxyxyxy iii −+ +−
=h
xyxyxp iii 2
)()()( 11 −+ − )()()( iii xrxyxq ++
iar condiţiile (25) devin ⎩⎨⎧
β=α=
)x(y)x(y
n
0
Obţinem în final schema:
(S) 211 2
hyyy iii −+ +−
-h
yyp iii 2
11 −+ −- 11 , −≤≤= niryq iii
βα == nyy ,0
unde )( ),( ),( ),( iiiiiiii xyyxrrxqqxpp ==== i)(∀ , care reprezintă un sistem liniar de n-1 ecuaţii ş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 soluţiei în punctele , hixi ⋅= 31 ≤≤ i , 4/1=h , cu ajutorul metodei
diferenţelor finite. Soluţie: Schema (S) ne conduce la sistemul următor:
⎪⎪⎪
⎩
⎪⎪⎪
⎨
⎧
−==
==⎟⎠⎞
⎜⎝⎛ −+⎟
⎠⎞
⎜⎝⎛ −−+⎟
⎠⎞
⎜⎝⎛ + +−
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ă, Bucureşti, 1997 2. R. Burden, J. Faires, Numerical Analysis, PWS-Kent, 2005. 3. C. Carasso, Analyse Numérique, Lidec, Canada, 1970. 4. P.Ciarlet, J.Lions, Finite Difference Methods, North-Holland, Amsterdam, 1989. 5. B. Demidovitch. I. Maron, Éléments de Calcul Numérique, Mir, Moscow, 1973. 6. D. Ebâncă, 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 Méthodes Numériques. Théorie et Applications Ed. Sitech, 2008 9. J.P. Nougier, Méthodes de Calcul Numérique, Hermes Sciences Publication, Paris, 2001. 10. M. Popa, R. Militaru, Metode Numerice în pseudocod–aplica�ii, 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
Top Related