C14-Rezolvarea numerica a ecuatiilor diferentiale.pdf

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))


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


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:


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.


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


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.


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


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


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


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


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


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


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


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


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


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


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


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 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ă.


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


|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


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


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.


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


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

C14-Rezolvarea numerica a ecuatiilor diferentiale.pdf

Documents

Transcript of C14-Rezolvarea numerica a ecuatiilor diferentiale.pdf