Cursul 14
Rezolvarea numerică a
ecuațiilor diferențiale
Introducere
Ecuaţiile diferenţiale sau cu derivate parţiale constituie
modelele matematice pentru majoritatea problemelor
inginereşti: studiul tensiunilor din: 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
impune o soluţie aproximativă.
Introducere
O ecuaţie diferenţială este o ecuaţie în care necunoscuta
este o funcţie ş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ă.
Introducere
Expresia generală a unei ecuaţii diferenţiale, sub formă
implicită este:
Unde y=y(x)
Forma explicită se poate scrie:
Există două tipuri de condiţii la limită :
1. condiţii Cauchy: se cunosc într-un punct x0 atât valoarea
funcţiei necunoscute cât şi valorile derivatelor, până la
ordinul cel mai mare ce figurează în ecuaţie;
2. se cunosc valorile funcţiei necunoscută în puncte diferite.
Integrarea ecuaţiilor diferenţiale cu condiţii iniţiale
Metode cu paşi separaţi
Se dau:
intervalul închis I=[x0, x0+a] R,
funcţia continuă f:IxRR,(x,y)f(x,y)
ecuaţia diferenţială P:y’ = f(x,y),
Problema diferenţială de ordinul 1 constă în determinarea funcţiei derivabile
y:IR, xy(x)
cu proprietatea că pentru xI avem y’(x)f(x,y(x))
Integrarea ecuaţiilor diferenţiale cu condiţii iniţiale
Pentru un sistem de ecuaţii diferenţiale de ordinul 1
se cunosc funcţiile continue:
fj:IxRnRn,(x,y1,…,yn) (z1,…,zn), j=1:n.
şi ecuaţiile diferenţiale
y1’=f1(x,y1,…,yn),
y2’=f2(x,y1,…,yn),
yn’=fn(x,y1,…,yn).
şi interesează determinarea funcţiilor derivabile
yj:IRn, xyj(x),
astfel încât yj’(x)fj(x,y1(x),…,yn(x)), j=1:n.
Integrarea unui sistem de ecuaţii diferenţiale de ordinul p
Integrarea ecuaţiilor diferenţiale cu condiţii iniţiale
yj(p)=fj(x,y1,…,yn,y1’,…,yn’,y1
(p-1),…,yn(p-1)), j=1:n
y(p)=f(x,y,y’,…,y(p-1)),
cu fj:Ix(Rn)pR, sau f:Ix(Rn)pRn.
presupune determinarea funcţiilor derivabile
yj:IRn, xyj(x)(sau y:xy(x))
Substituţiile:
y=v1,
y’=v2,
y(p-1)=vp.
reduc sistemul la ordinul 1:
Integrarea ecuaţiilor diferenţiale cu condiţii iniţiale
v1’=v2,
v2’=v3,
vp’=f(x,v1,v2,…,vp).
Problema diferenţială cu condiţii iniţiale (problema Cauchy) :
P1:y’=f(x,y),
P2:y(x0)=, cu R, dat (condiţia iniţială).
Presupunem că funcţiile f satisfac o condiţie Lipschitz
xI, u,vRn,L>0, astfel încât
|f(x,u)-f(x,v)|<L.|u-v|.
Condiţia Lipschitz asigură existenţa şi unicitatea
soluţiei.
Integrarea ecuaţiilor diferenţiale cu condiţii iniţiale
Problema diferenţială cu condiţii iniţiale este echivalentă cu determinarea unei funcţii y continue pe I, cu proprietatea că:
echivalenţă cunoscută şi sub numele de metoda constructivă
Picard.
Metoda lui Euler
Se împarte intervalul I=[x0,x0+a] în N intervale echidistante
de lungime prin punctele de abscise
şi se aproximează soluţia y(xi)Yi , în care
dt)t(y,tf)x(y
x
0x
N
ah
N
aixihxx 00i
.1N:0i,y,xhfYY
,Y
iii1i
0
Metode cu paşi separaţi: Metoda lui Euler
y = Euler(x0, a, n, f, y0)
% Intrări:
%x0=capătul stâng al intervalului de integrare
% a = lungimea intervalului
% n = numărul de puncte
% y0 = condiţia iniţială
% f = funcţia de integrat y'=f(x,y)
%Ieşiri:
% y = vectorul aproximaţiilor soluţiei
y(0) = y0 h = a / n
for i = 1 : n
x = x0 + i*h
y(i) = y(i-1) +h*f(x, y(i-1))
end
Metode cu paşi separaţi: Metoda lui Euler
Se consideră problema Cauchy:
Să se determine o soluţie aproximativă a acestei probleme folosind
metoda Euler cu pasul h = 0,2.
Rezolvare
Folosind formulele din slidul 9 pentru x0 = 0; x1 = 0,2; x2 = 0,4; x3 =
0,6; x4 = 0,8; x5 = 1 şi y0 = 1, se obţine:
y1 = y0 + hf(x0, y0); f(x0, y0) = y0 – 2x0/y0= 1.
Deci y1 = 1 + 0,2 · 1 = 1,2.
y2 = y1+0,2f(x1, y1) = y1+0,2 (y1 – 2x1/y1)= 1,2 + 0,2·0,8667 =
1,3733.
Samd.
Metode cu paşi separaţi: Metoda lui Euler
Soluția exactă este
Rezultatele comparative ale metodei axproximative cu cea exactă:
xi yi (Euler) yi (exact)
0 1 1
0,2 1,2 1,1832
0,4 1,3733 1,3416
0,6 1,5294 1,4832
0,8 1,6786 1,6124
1 1,8237 1,7320
12 xy
Metode cu paşi separaţi: Metoda lui Euler
Consistenţă, stabilitate, convergenţă
O metodă cu paşi separaţi determină aproximaţia soluţiei în pasul următor Yi+1, folosind numai
informaţia din pasul curent i.
Eroarea într-un pas (în xi) se defineşte ca:
ei=y(xi)-Yi,
iar eroarea globală
,1,0i,Y,xhfYY
,Y
iihi1i
h0
i
N:0i
N emaxE
Metode cu paşi separaţi
O metodă cu paşi separaţi este consistentă dacă
Stabilitatea unei metode cu paşi separaţi impune
ca variaţia condiţiilor iniţiale să nu producă variaţii
mari în rezultate.
Fie problema diferenţială cu condiţii iniţiale
y’=f(x,y), y(x0)=,
problema diferenţială perturbată
z’=f(x,z)+(x), z(x0)=y0+0
şi metodele cu paşi separaţi corespunzătoare:
y,xfy,xflim h0h
Metode cu paşi separaţi
yi+1=yi+h.f(xi,yi),
zi+1=zi+h.[fh(xi,yi)+iN].
Metoda cu paşi separaţi este stabilă dacă
astfel încât:
y(0)=y0,
y’(x)=-A.y(x), A>0
y(x)=y0e-Ax (soluţia exactă)
Metoda Euler furnizează aproximaţia
Yn=(1-h.A)n.Y0,
Deoarece este necesar ca:
0K,K,0lim21iN
N
iN
Ni0
2001iiNi0
maxKzyKzymax
0tylimt
Metode cu paşi separaţi: metode Runge-Kutta
care ne conduce la:
adică metoda Euler este condiţional stabilă, aceasta
furnizează rezultate corecte numai dacă se alege
pasul suficient de mic. Condiţia de mai sus defineşte
A-stabilitatea metodei.
Metode de tip Runge – Kutta
O metodă de tip Runge-Kutta este o metodă cu paşi separaţi în care funcţia fh(x,y) se determină astfel:
se împarte intervalul [xi,xi+1] în q subintervale
cu abscisele:
0Ylim nn
A
2h01Ah1
Metode Runge-Kutta
xij=xi0+uj.h. 0uj1, , u0=0, uq=1
Numărul subintervalelor q defineşte rangul
metodei.
se calculează aproximaţiile soluţiei în punctele
intermediare de forma:
Punctele intermediare xij şi constantele Kjl se
obţin din condiţia ca în dezvoltarea Taylor a lui yij
după puterile lui h, să coincidă cât mai mulţi
termeni cu cei din dezvoltarea Taylor a soluţiei
exacte.
1j
0l
ililjliij
i0i
.q:1j,y,xfKhyy
,yy
Metode Runge-Kutta
Metoda este de ordinul p, dacă în cele două
dezvoltări termenii coincid până la hp inclusiv.
Metoda Runge-Kutta de ordin 1 şi rang 1 este de
forma:
Dezvoltarea în serie Taylor a soluţiei exacte este
Din identificare se obţine K10=1,care conduce la
metoda lui Euler. Metoda nu este utilizată în mod
practic, deoarece eroarea într-un pas, de ordinul h2, este importantă.
.y,xfKhyy
,yy
0i0i10i1i
i0i
0i0ii
2
iii1i y,xhfy!2
hy
!1
hyhxyxy
Metode Runge-Kutta
In metodele cu paşi separaţi de rang 2 şi ordin 2 se ia
xij=xi0+ujh, j=0:2, u0=0, u1[0,1], u2=1
Pentru j=1 avem yi1=yi+hK10f(xi0,yi0)
Dezvoltarea Taylor a soluţiei este
y(xi1)=y(xi0+u1h)=yi+u1hyi’+…=yi+u1hf0+…
Din identificare rezultă K10=u1.
Pentru j=2 avem:
yi2=yi+hK20f(xi0,yi0)+hK21f(xi1,yi1)
.2,1j,y,xfKhyy
,yy1j
0l
ililjliij
i0i
Metode Runge-Kutta
în care f(xi1,yi1) se înlocuieşte cu dezvoltarea
în serie Taylor în vecinătatea punctului (xi0,yi0)
y
fyy
x
fxxy,xfy,xf
0
0i1i
0
0i1i0i0i1i1i
y
fhfK
x
fhufy,xf
0
010
0
101i1i
y
ffKKh
x
fKhuhfKKyy
0
01021
20
21
2
102120i2i
i
2
iii2i y2
hyhyhxyxy
y
ff
x
f
2
hhfy
0
0
0
2
0i
Metode Runge-Kutta
K20 + K21=1,
yi0=yi,
yi1=yi+hu1f(xi0,yi0),
Cazuri particulare
metoda tangentei ameliorate,
2
1Ku 211
1
21
1
20u2
1K,
u2
11K
1i1i
1
0i0i
1
i2i y,xfu2
hy,xf
u2
11hyy
2
1u1
2
hxhuxx i10i1i
Metode Runge-Kutta
yi+1=yi + hf(xi1,yi1).
- metoda Euler-Cauchy, în care u1=1
xi1=xi+h
yi1=yi+hf(xi,yi)
metoda Heun
iii1i y,xf
2
hyy
1i1iiii1i y,xfy,xf
2
hyy
3
2u1
h3
2xx i1i
1i1ii1i y,xhf3
2yy
1i1iiii1i y,xf
4
h3y,xf
4
hyy
Metode Runge-Kutta
Eroarea într-un pas în metodele Runge-Kutta de
ordin 2 este
Metoda Runge-Kutta folosită are ordin 4:
K1=h.f(xi,yi),
y
fyu
2
3yu
2
31
6
h11
3
2
Ky,
2
hxfhK 1
ii2
2
Ky,
2
hxfhK 2
ii3
3ii4 Ky,hxfhK
6
KK2K2Kyy
4321
i1i
Metode Runge-Kutta
Pentru ca metoda Runge-Kutta
yi+1=yi+h.(xi,yi,h)
să fie convergentă către soluţia y(t) a problemei
diferenţiale, atunci când h0 este suficient ca
să verifice o condiţie Lipschitz globală şi
(x,y,0)=f(x,y).
Verificarea A-stabilităţii pentru metoda Runge-
Kutta de ordinul 4 conduce la
K1=h.f(tn,yn)=-hAyn,
2
hA1hAy
2
KyhA
2
Ky,
2
hthfK n
1
n
1
nn2
4
Ah
2
hA1hAy
2
Ky,
2
hthfK
22
n
2
nn3
Metode Runge-Kutta
yn+1=p4(hA).yn,
yn=[p4(hA)]n
Pasul de discretizare nu poate fi ales arbitrar
4
Ah
2
AhhA1hAyKy,hthfK
3322
n3nn4
,hA24
1hA
6
1hA
2
1hA1y
6
KK2K2Kyy
432
n
4321
n1n
0ylim0tylim nnt
A
78.2h01hAp4
Metode cu paşi legaţi (Metode multipas)
Metodele cu paşi separaţi
sunt simple
necesită puţine informaţii iniţiale
prezintă dezavantajul lipsei de precizie.
Metodele cu paşi legaţi (sau metodele multipas)
folosesc mai multe informaţii iniţiale,
sunt mai precise.
Problema diferenţială cu condiţii iniţiale
în care f(x,y) satisface condiţia Lipschitz globală
|f(x,u)-f(x,v)|L|u-v|
.y,xfy
,xy 0
Metode cu paşi legaţi (Metode multipas)
este echivalentă cu determinarea funcţiei
Presupunem cunoscute valorile aproximative
yky(xk), fk=f(xk,yk),
yk-1y(xk-1), fk-1=f(xk-1,yk-1),
yk-ry(xk-r), fk-r=f(xk-r,yk-r).
dtty,tfxyx
x0
,dtty,tfxy
,dtty,tfxy
k
0
1k
0
x
xk
x
x1k
dtty,tfxyxy1k
k
x
xk1k
Metode cu paşi legaţi (Metode multipas)
Folosim cea de-a 3-a formulă de interpolare
Newton-Gregory :
Integralele calculate cu
metoda seriei generatoare au valorile
dxxxxx!1r
y,f
fj
u1dxxy,xf
rkk
x
x
kk
1r
k
jx
x
r
0j
jx
x
1k
k
1k
k
1k
k
.duru1uu!1r
y,fh
duj
u1hf
1
0
kk
1r2r
1
0
jr
0j
k
j
duk
u1c
1
0
k
k
Metode cu paşi legaţi (Metode multipas)
duc la relaţia
definesc o metodă explicită, cunoscută sub
numele de metoda Adams-Bashforth.
In acelaşi mod se obţine metoda implicită Adams-
Moulton
Dezvoltarea diferenţelor regresive conduce la
relaţii de recurenţă liniare de forma:
,720
251,
8
3,
12
5,
2
1,1
k
r
0j
j
jk1k fchyy
1k
r
0j
j
jk1k fdhyy
Metode cu paşi legaţi (Metode multipas)
O metodă explicită are 0=0 în timp ce metodele
implicite se caracterizează prin 00.
Pentru determinarea coeficienţilor j, j vom
impune ca relaţia (28) să reprezinte soluţia exactă
pentru polinoame de grad cât mai ridicat.
Astfel pentru
y(x)=1, y’(x)=0,
y(x)=x, y’(x)=1,
y(x)=xp, y’(x)=pxp-1.
se aleg punctele xk-r,…,xk echidistante cu h=1 şi
originea în xk-r.
jk
r
0j
jjk
r
1j
jk fhyy
Metode cu paşi legaţi (Metode multipas)
Determinarea celor 2r+1 coeficienţi impune tot
atâtea relaţii, aşadar formula va fi exactă pentru
polinoame de grad 2r.
Metode explicite şi implicite
Pentru r=2, metoda Adams-Bashforth este
1r
1j
j
rjrr
0j
j
r
1j
j
pr
0j
1p
j
r
1j
p
j rjrpjr
2
fffhyf
2
1fhyy 1kk
kkkkk1k
Metode cu paşi legaţi (Metode multipas)
pentru r=3
pentru r=4
1kkk1k ff3
2
hyy
k
2
1k yh12
5e
2k1kkk1k f5f16f23
12
hyy
k
43
1k yh8
3e
3k2k1kkk1k f9f37f59f55
24
hyy
k
54
1k yh720
251e
Metode cu paşi legaţi (Metode multipas)
pentru r=5
Metodele implicite Adams-Moulton sunt :
r=2
r=3
4k3k2k1kkk1k f251f1274f2616f2774f1901
720
hyy
k
65
1k yh288
95e
1kk1kk1k ff8f
12
hyy
k
43
1k yh24
1e
2k1kk1kk1k ff5f19f9
24
hyy
k
54
1k yh720
19e
Metode cu paşi legaţi (Metode multipas)
r=4
metoda explicită Milne
metoda implicită Simpson
3k2k1kk1kk1k f19f106f264f646f251
720
hyy
k
65
1k yh160
3e
2k1kk3k1k f2ff2
3
h4yy
k
54
1k yh45
14e
1kk1k1k1k ff4f2
3
hyy
k
54
1k yh90
1e
Metode cu paşi legaţi (Metode multipas)
Metode predictor – corector.
Metodele implicite asigură aproximaţii mai bune
ale soluţiilor decât metodele explicite. De aceea
ele se folosesc pentru a creşte precizia
rezultatelor obţinute prin metode explicite.
Aplicarea unei metode implicite conduce la o
ecuaţie neliniară de forma
yk+1=hC1f(xk+1,yk+1)+C2.
Rezolvarea ecuaţiei neliniare se face de obicei
prin metoda aproximaţiilor succesive.
Metoda explicită calculează o predicţie a soluţiei,
iar metoda implicită - o corecţie, făcând o singură
iteraţie, cu valoarea de pornire, predicţia oferită de
metoda explicită.
Metode cu paşi legaţi (Metode multipas)
Metoda explicită Adams-Bashforth de ordinul 3
furnizează predicţia
formula implicită Adams-Moulton de ordinul 2
corectează această valoare
Convergenţa şi stabilitatea metodelor multipas
Dacă yk reprezintă o estimare obţinută printr-o
metodă multipas a soluţiei exacte y(xk), atunci se
spune că metoda multipas este de ordinul p dacă
2k1kkk
p
1k f5f16f2312
hyy
1kk
p
1kk
c
1k yy8y512
hyy
Metode cu paşi legaţi (Metode multipas)
|y(xk)-yk|C.hp+1 .
O metodă multipas generală aproximează soluţia
problemei diferenţiale cu condiţii iniţiale cu soluţia
ecuaţiei recurente
Metodele explicite se iau cu 0=0, în timp ce
pentru metodele implicite 00.
Astfel metodei explicite Adams-Bashforth
şi metodei implicite Adams-Moulton
r
0j
r
0j
jkjjkj fhy
r
1j
jkj1kk fhyy
r
0j
jkj1kk fhyy
Metode cu paşi legaţi (Metode multipas)
0=1, 1=-1, j=0, j=2:r,
Metoda explicită Nystrom
şi metoda implicită Simpson-Milne
folosesc 0=1, 1=0, 2=-1, j=0, j=3:r. Se definesc polinoamele de gradul k
r
1j
jkj2kk fhyy
r
0j
jkj2kk fhyy
r
0j
jk
jzz
r
0j
jk
jzz
Metode cu paşi legaţi (Metode multipas)
Ecuaţia recurentă poate fi exprimată ca
(E)yk-h.(E)fk=0, unde Eyk=yk+1
Algoritmul multipas definit prin relaţia recurentă este stabil dacă
- rădăcinile ecuaţiei (z)=0 sunt de modul subunitar i.e.
(zi)=0 |zi|1;
- rădăcinile de modul 1 sunt simple
Stabilitatea este puternică dacă singura rădăcină de modul 1 este zi=1.
Metoda Adams este puternic stabilă deoarece
(z)=zk-zk-1=0,
are singura rădăcină de modul 1 zi=1.
Metode cu paşi legaţi (Metode multipas)
Metoda Nystrom este slab stabilă căci
(z)=zk-zk-2=0.
are două rădăcini de modul 1 zi=1.
Impunem ca ecuaţia recurentă să aibă ordinul p.
Se dezvoltă această ecuaţie în serie Taylor în vecinătatea unui punct t
în care
C0=0+1+…+r=(1)
C1=1+22+…+rr-(0+1+…+r)=’(1)-(1)
r
1p
2
1p
1r
p
2
p
1p r2!1p
1r2
!p
1C
tyhCtyEhtyEii
0i
i
Metode cu paşi legaţi (Metode multipas)
Dacă metoda este de ordin p atunci
O metodă cu paşi legaţi este consistentă dacă
este de ordin cel puţin 1.
Condiţia de consistenţă impune deci ca
O metodă cu paşi legaţi stabilă şi consistentă este
convergentă.
.0C
,0CCC
1p
p21
.11
,01