INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE CU DERIVATE PARTIALE

8/7/2019 INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE CU DERIVATE PARTIALE

1/8

Carmen-Sanda Georgescu, Tudor Petrovici, Radu PopaMetode numerice. Fisa nr. 13:

INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE CUDERIVATE PARTIALE

8 . INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE CU DERIVATE PARTIALEDE ORDINUL DOI CU CONDITII LA LIM ITA

Fie o ecuatie diferentiala cu derivate partial e de ordinul doi in care functia u depinde de doua variabileindependente x si y este de forma:

0),,,,),,(,,(2

2

2

2

2

=

yx

u

y

u

x

u

y

u

x

uyxuyxF (13.1)

Ecuatia diferentiala cu derivate partiale de ordinul doi (13.1) se scrie sub forma

0),),,(,,(),(),(),(2

22

2

2

=

+

+

+

y

u

x

uyxuyxf

y

uyxC

yx

uyxB

x

uyxA (13.2)

Daca termenul ),),,(,,(

y

u

x

uyxuyxf

este liniar atunci ecuatia (13.2) se numeste liniar a ; daca

termenul ),),,(,,(y

u

x

uyxuyxf

este neliniar atunci ecuatia (13.2) se numeste cvasiliniara. Ecuatia

(13.2) caracterizeaza fenomene fizice care se desfasoara intr-un dom eniu2

RD si este insotita deconditii initiale si/ sau la limita pe frontiera dupa caz.Ecuatiile diferentiale cu derivate partiale de ordin ul doi se clasifica in functie de semnul discriminant ului

ACB 42 =

Daca 0 ecuatia este de tiphiperbolic (semnul discriminantulu i trebuie sa fie acelasi in toate punctele domeniulu i ; in caz contrar esteposibil ca ecuatia sa-si schimbe tipul in diferite subdomenii ale dom eniului D). Rezolvarea numerica aecuatiilor diferentiale cu derivate partiale impune transform area acestora in ecuatii cu diferente finite(ecuatii algebrice).

8.1. OPERATORI CU DI FERENTE FINI TE PARTIALE IN COORDONATE CARTEZIENE

Fie2RD domeniul in care are loc fenomenul fizic. Transformarea ecuatiei diferentiale cu derivate

partiale de ordinul doi (13.2) in ecuatii algebrice cu diferente finite in coordinate carteziene, impunediscretizarea domeniul D intr-o retea cu subdomen ii dreptunghiulare (in particular patrat ice). Daca notamcu hdistantele constante dintre nodurile ret elei in directia x si cu k distantele constante in directia y dintrenodurile retelei (Fig.1)

y

i-1,j+1 i,j+1 i+1,j+1

i-1,j i,j i+1,j x

i-1,j-1 i,j-1 i+1,j-1

Fig.1. Discretizarea domen iului2

RD in coordonate carteziene. Indicele i indica coloana iar indicele jindica linia

Operatori cu diferente finite partiale centrate

k

uu

y

u

h

uu

x

u jijiji

jiji

ji2

)(,2

)(1,1,

,

,1,1

,

++

(13.3)


2/8

hk

uuuu

yx

u

k

uuu

y

u

h

uuu

x

u

jijijiji

ji

jijiji

ji

jijiji

ji

4)(

2)(,

2)(

1,11,11,11,1

,

2

2

1,,1,

,2

2

2

,1,,1

,2

2

++++

++

+

+

+

(13.4)

8.2. ECUATII DIFERENTIALE CU DERIVATE PARTIALE DE TIP ELIPTICECUATIA LUI LAPLACE

O mare varietate de pr obleme bidimen sionale din fizica si tehnica sunt descrise de ecuatia lui Laplace. Informa canonica ecuatia lui La place se scrie sub forma

datauy

u

x

u==

+

|,02

2

2

2

(13.5)

unde este frontiera domeniului inchis D ( o astfel de problema se numeste problema Dirichlet).Conditii le la limita pentru ecuatia (13.5) au forme specifice functie de fenomenul fizic modelat mate maticprin acest tip de ecuatie.Exemplul 1. Se considera o placa subtire complet izolata si T(x,y) temperatura placii intr -un punct M(x,y) aldomeniului plan (dreptunghi). In ipoteza ca temperatura este independenta de timp sa se determinedistributia temperatu rii in interiorul placii daca pe frontiera domeniului temperat ura are valorile : T=100oCpe latura din stanga si T=0oC pe celelalte trei laturi. Se va considera (Fig.2) o discretizare a domeniului insubdomenii patrate cu 6 noduri interioare (notate cu 1,2,3,4,5,6) si 10 noduri pe frontiera echistante (inaceste noduri se cunoaste tem peratura T conform condit ii lor la limita dat e ; T7=0, T8=0 , T9=0, T10=0,T11=0 ,T12=0, T13=0, T14=0 ,T15=100, T16=100).Rezolvare : Daca se introduc aproximatiile (13.4) pentru derivatele partiale de ordin ul doi in ecuatia (13.5)se obtine ecuatia algebrica in diferente finite (Fig.1)

4

,1,11,1,

,

jijijiji

ji

uuuuu

++ +++ (13.6)

T7=0 T8=0 T9=0

T16=100 1 2 3

T15=100 4 5 6

T10=0

T11=0

T14=0 T13=0 T12=0Fig. 2. Discretizarea domeniul ui de integrare

Aplicand relatia (13.6) in fiecare din cele 6 noduri interioare si notand t emperaturile din aceste noduri prinT1, T2, T3, T4, T5, T6 se obtine un sistem de sase ecuatii cu sase necunoscute

=+

=+

=+

=+

=+

=

0404

1004

04

04

1004

653

6542

541

632

5321

421

TTTTTTT

TTT

TTT

TTTT

TTT

(13.7)

Rezolvarea sistemului (13.7) se poate face prin mai multe metode : Pivotare Gauss. Daca se rezolva sistemul prin pivotare Gauss se obtin solutiile:T1=38.0952 , T2=14.2857, T3=4.7619, T4=38.0952, T5=14.2857, T6=4.7619 Metoda iterativa. Se observa ca sunt indeplinite conditii le de rezolvare a sistemului prin metodeleiterative Jacobi si Gauss-Seidel (matricea sistemului este diagonal dominant a adica


3/8

niaa ii

n

ijj

ij ,...,2,1|,|||1

==

). Aceste metode sunt preferate atun ci cand discretizarea domeniului de

integrare contine un numar mar e de noduri interioare Metoda iterativa LiebmannLiebmann a propus o metoda iter ative cu accelerarea convergentei procesului de iteratie. Astfel el

introduce un coefficient w cu valori cuprinse intre 1 si 2 in ecuatia cu diferente finite. Coeficientul poate fischimbat la fiecare pas de iteratie. Astfel procedeul lui Liebman n poate fi descris prin relatia

[ ])(4

1 )(1,

)(

1,

)(

,1

)(

,1

)(

,

)(

,

)1(

,

k

ji

k

ji

k

ji

k

ji

k

ji

k

ji

k

ji uuuuuwuu +++ +++= (13.8)

Exemplul 2. Sa se scrie un program Octave pentru probleme de tip (13.5) folosind procedeul lui Liebmannsi sa se rezolve problema : Se considera o placa subtire complet izolata si T(x,y) temperatura placii intr-unpunct M(x,y) al domeniului plan de forma pat rata ( x[0 ;4], y[0 ;4]). In ipoteza ca tem peratura esteindependenta de timp sa se determine distribu tia temperaturii in interior ul placii daca pe frontieradomeniului tem peratura are valorile : T=20 oC pe latura din stanga , T=180 oC pe latura din dreapt a, T=80 oCpe latura de sus si T=0oC pe latura de jos. Se va considera (Fig.2) o discretizare a domeniului insubdomenii patrate cu 64 nodur i interioare . Sa se reprezinte grafic functia u(x,y)Algoritm:

1. Informatia initiala : Citeste functiile f1(x), f2, f3, f4 si marimile a=lungimea intervalului pe Ox, b=lungimea intervalului pe Oy, h=lungimea pasului, tol=toleranta admisa, max=numarul maxim de iteratii

2. Calculeaza numarul nodurilor pe Ox si respectiv pe Oy, n fix(a/h)+1;mfix(b/h)+1; si valoarea destart a procesului de iteratie media (a*(feval(f1,0)+feval(f2,0)) + b*(feval(f3,0)+ feval(f4,0)))/(2*a+2*b);

3. Defineste parametrul w 4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2));4. Pentru j=2:m-1 si i=2:n-1 calculeaza u(i,j) dupa formula (13.8)

Functie - Program Octave

octave#1>function u = liebmann(f1, f2, f3, f4,a,b,h,tol,max)n = fix(a/h)+1;m = fix(b/h)+1;media = (a*(feval(f1,0)+feval(f2,0)) + b*(feval(f3,0)+feval(f4,0)))/(2*a+2*b);u = media*ones(n,m);for j=1:m,

u(1,j) = feval(f3,h*(j-1)); u(n,j) = feval(f4,h*(j-1));endforfor i=1:n,

u(i,1) = feval(f1,h*(i-1)); u(i,m) = feval(f2,h*(i-1));endforu(1,1) = (u(1,2) + u(2,1))/2;u(1,m) = (u(1,m-1) + u(2,m))/2;u(n,1) = (u(n-1, ) + u(n,2))/2;u(n,m) = (u(n-1,m) + u(n,m-1))/2;1err = 1;cnt = 0;w = 4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2));while ((err>tol)&(cnt


4/8

10.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 90.00000

Fig.3. Reprezentarea grafica a repartitiei tem peraturii in placa plana de la Exemplul 2

8.3. ECUATII DIFERENTIALE CU DERIVATE PARTIALE DE TIP PA RABOLI CECUATIA CALDURII

Consideram o bara subtire complet izolata (Fig.4) cu o distributie initiala a temperaturii la timpul t= 0 data(vom considera variabila x adimen sionala). La x=0 si respectiv x=1 si t>0 cu

se dau conditii le la limita

]1,0[),()0,( = xxfxu

],0( Tt

)(),1(),(),0(10tgtutgtu == . In forma canonica ecuatia caldurii se scrie sub forma

02

2

=

t

u

x

uc (13.9)

cu conditii le initiale si la limite

=

=

=

],0(),(),1(

],0(),(),0(

]1,0[),()0,(

1

0

Tttgtu

Tttgtu

xxfxu

(c este difuzibilitatea termica a barei)

Daca se considera o discretizare pentru timp un pas t , iar pentru x un pas x si se aplica ecuatiei (13.9)

operatorii in diferente finitet

uu

t

u nini

+ ,1,si (13.4) pentru tn se obtin ecuatiile algebrice in diferente

finite sub forma explicita

2

,1,,1,1,

)(

2

x

uuu

t

uu ninininini

+=

++(13.10)


5/8

(i este indicele pe directia axei Ox, n este indicele pe axa Ot).Daca definim constanta de discretizare a

domeniului2

)( x

t

= atunci ecuatiile (13.10) se scriu sub forma explicita

nininini uuuu ,1,,11, )21( ++ ++= (13.11)cu conditii le la limite si initiale uo,n+1=g o(t n+ 1), u 1,n+1=g 1(t n+1 ), u i,o=f(x i)

Daca se considera o discretizare pentru timp un pas t , iar pentru x un pas x si se aplica ecuatiei (13.9)

operatorii in diferente finitet

uu

t

u nini

+ ,1,si (13.4) pentru tn+ 1 se obtin ecuatiile algebrice in

diferente finite sub forma implicita

2

1,11,1,1,1,

)(

2

x

uuu

t

uu ninininini

+=

+++++(13.12)

(i este indicele pe directia axei Ox, n este indicele pe axa Ot).Daca definim constanta de discretizare a

domeniului2

)( x

t

= atunci ecuatiile (13.10) se scriu sub forma implicita

1,11,1,1,)21( ++++ ++= nininini uuuu (13.13)

cu conditii le la limite si initiale u o,n+1=g o(t n+ 1), u 1,n+1=g 1(t n+ 1), u i,o=f(x i)t

T

u(0,t)=g0(t) u(1,t)=g1(t )

0 u(x,0)=f(x) xFig.4

Exemplul 4. Sa se scrie un program Octave care sa rezolve problema (13.9) si sa se rezolve problema

(13.9) daca , ,a=1, f(x)=x*sin(x), g]1,0[x ]1,0[t 0=0 oC, g 1=0oC . Se considera o discretizare aintervalului pe x de pas h=0.2 si discretizare a intervalului pe t de pas k =0.2

octave#1>function U = caldura(f,g1,g2,a,b,c,n,m)h = a/(n-1);k = b/(m-1);r = c^2*k/h^2;s = 1 - 2*r;U = zeros(n,m);for j = 1:m,U(1,j) = feval("g1",k*(j-1)); U(n,j) = feval("g2",k*(j-1));

endforfor i = 2:(n-1),U(i,1) = feval("f",h*(i-1));

endfor

for j = 2:mfor i = 2:(n-1)U(i,j) = s*U(i,j-1) + r*(U(i-1,j-1) + U(i+1,j-1));

endforendforendfunction

Program Octave :octave#1> function z = f(x),z = x.*sin(x);endfunctionoctave#2>function z = g1(x),z = 0;endfunctionoctave#3>function z = g2(x),z = 0;endfunctionoctave#4>u=caldura("f","g1","g2",1,1,1,5,5)


6/8

u =0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+006.1851e-02 5.2589e-01 -1.2239e+00 -4.2135e+01 9.6302e+022.3971e-01 6.1433e-01 -1.2676e+01 1.6702e+02 -2.1227e+035.1123e-01 -2.6198e+00 2.0796e+01 -1.9627e+02 2 .0420e+030.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00

8.4. ECUATII DIFERENTIALE CU DERIVATE PARTIALE DE TIP HIPERBOLICECUATIA COARDEI VIBRAN TEO coarda elastica perfect flexibila de lungime L este intinsa intre doua pun cte capete x=0 si x=L intr-unplan orizontal. Coarda are o sageata initiala verticala u(x,0)=f(x) si este perturbata din p ozitia initiala dereapus la timpul t=0 pastrand capetele fixe. Sa se determine sageata u(x,t) la un timp t>0 (pentr u oscilatiimici ale coardei). Ecuatia care descrie acest fenomen se scrie sub forma

02

2

2

22 =

t

u

x

ua (13.14)

(

va = este viteza de propagare a un delor in lungul coardei, cu v reprezantand tensiunea la capetele

coardei si mas ape unitatea de lungime)

Conditii le initiale si la limita atasate ecuatiei (13.14) : u(x,0)=f(x), )()0,( xgxt

u

=

, u(0,t)=0,u(L,t)=0.

Daca se considera o discretizare pentru timp un pas t , iar pentru x un pas x si se aplica ecuatiei

(13.14) operatorii in diferente finite2

1,,1,

2

2

)(

2

t

uuu

t

u jijiji

+

+si operatorul (13.4) se obtin

formulele de recureanta in diferente finite pentru determinarea deplasarii in punctul xi la timpul tj+1

(13.15)2,1,11,12

1,

2

,)1(2 + ++= jijijijiji uuuuu

undex

ta

= . Solutia in diferente finite pentru ecua tia (13.13) este stabila pentru 1 .

Exemplul 4. Sa se scrie un program Octave care sa rezolve problema (13.14) si sa se rezolve problema(13.13) daca L=1, t(0 ;0.5], a=2, f(x)=sin(pi*x)+sin(2*pi* x), g(x)=0.

Algoritm:

5. Informatia initiala : Citeste functiile f(x), g(x) si marimile L, tfinal, a, numarul nodurilor pe directia x intrex=0 si x=L, numarul intervalelor de timp dintre t=0 si t=tfinal

6. Calculeaza pasii h pe Ox si respectiv k pe Ot, h L/(n-1); ktfinal/(m-1). Calculeaza lambda a*k/h;7. Pentru i=2:n-1, calculeaza formulele de pornire pentru u(i,1) si u(i,2) folosind conditiile initiale8. Pentru j=3:m si i=2:n-1 calculeaza u(i,j) dupa formula (13.15)

Functie - Program Octave

octave#1> function u = coarda(f,g,L,tfinal,a,n,m)h = L/(n-1);k = tfinal/(m-1);lambda = a*k/h;lambda2 = lambda^2;u = zeros(n,m);

for i=2:(n-1),u(i,1) = feval(f,h*(i-1));u(i,2) = (1 lambda2)*feval(f,h*(i-1)) + k*feval(g,h*(i-1)) ...

+ (lambda2/2)*(feval(f,h*(i)) + feval(f,h*(i-2)));endforfor j=3:m,

for i=2:(n-1)u(i,j) = (2-2*lambda2)*u(i,j-1) + lambda2*u(i-1,j-1) + u(i+1,j-1)) - u(i,j-2)

endforendforendfunction


7/8

Program Octaveoctave#1>function z = f(x),z = sin(pi*x) + sin(2*pi*x);endfunctionoctave#2>function z = g(x),z = 0;endfunctionoctave#3>u=coarda(f, g,1,0.5,2,11,11)Rezultate:

x=0 x=0.1 x=0.2 x=0.3 x=0.4 x=0.5 x=0.6 x=0.7 x=0.8 x=0.9 x=1t=0.0 0 0.89680 1.53884 1.76007 1.53884 1.00000 0.36327 -0.14204 -0.36327 -0.27877 0t=0.05 0 0.76942 1.32844 1.53884 1.38004 0.95106 0.42898 -0.00000 -0.21040 -0.18164 0

t=0.1 0 0.43164 0.76942 0.94840 0.95106 0.80902 0.58779 0.36062 0.18164 0.06836 0t=0.15 0 0.00000 0.05160 0.18164 0.37738 0.58779 0.74065 0.76942 0.63938 0.36327 0t=0.2 0 -0.38004 -0.58779 -0.51942 -0.18164 0.30902 0.76942 1.01942 0.95106 0.57102 0t=0.25 0 -0.58779 -0.95106 -0.95106 -0.58779 -0.00000 0.58779 0.95106 0.95106 0.58779 0t=0.3 0 -0.57102 -0.95106 -1.01942 -0.76942 -0.30902 0.18164 0.51942 0.58779 0.38004 0t=0.3.5 0 -0.36327 -0.63938 -0.76942 -0.74065 -0.58779 -0.37738 -0.18164 -0.05160 -0.00000 0t=0.4 0 -0.06836 -0.18164 -0.36062 -0.58779 -0.80902 -0.95106 -0.94840 -0.76942 -0.43164 0t=0.45 0 0.18164 0.21040 0.00000 -0.42898 -0.95106 -1.38004 -1.53884 -1.32844 -0.76942 0t=0.5 0 0.27877 0.36327 0.14204 -0.36327 -1.00000 -1.53884 -1.76007 -1.53884 -0.89680 0

Fig.5 Repartitia sagetii coardei vibrante din Exemplul 4


8/8

INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE CU DERIVATE PARTIALE

Documents

Transcript of INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE CU DERIVATE PARTIALE