INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE ORDINARE CU CONDITII LA LIMITA (Probleme...

8/7/2019 INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE ORDINARE CU CONDITII LA LIMITA (Probleme Sturm-Liouville s

1/6

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

INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE ORDINARECU CONDITII LA LIMITA (Probleme Sturm-Liouville sau bilocale)

7 . 4 . INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE ORDINARE DEORDINUL II, CU CONDITII LA LIMITA

Ecuatia diferentiala de ordinul doi, definita pe intervalul x[a, b]

[ ] 0)()()( =++

yxrxq

dx

dyxp

dx

d

(12.1)

cu conditiile la limita (12.2)

=+

=+

221

121

)(')(

)(')(

bybbyb

ayaaya

unde p(x)>0,r(x)>0,q(x) functie de pondere, a1, a2, b1, b2, 1, 2 constante date si parametru real

formeaza o problema Sturm-Liouville. Problemele Sturm-Liouville se intalnesc in multe aplicatii cum suntspre exemplu functiile Bessel sau functiile Legendre .Sistemul diferential (12.1) si (12.2) se mai numeste si problema bilocala.

7 . 4 . 1 . METODA TIRULUIMETODA TIRULUI reduce problema bilocala (12.1) si (12.2) la o problema Cauchy, rezolvabila cu una dinmetodele prezentate in Fisa_10 sau Fisa_11.Pas1: Se aleg doua valori y0 si y0 care verifica prima din conditiile la limita (12.2) si se rezolva problemaCauchy formata din ecuatia diferentiala (12.1) cu conditiile initiale y(a)= y0, y(a)= y0Pas2: Cu solutia aproximativa obtinuta verificam egalitatea a doua din conditiile la limita (12.2)Daca egalitatea este verificata cu o aproximatie convenabila , se opreste algoritmul si solutia gasita estesolutia cautata, daca nu se continua algoritmul cu pasii 1 si 2 pana cand egalitatea este verificata cu oaproximatie convenabil aleasa.Exemplul 1. Sa se afle solutia aproximativa pe intervalul x[0,1], a problemei bilocale neliniare y'' = 0.5*y-2*(y)^2/yy(0) = 1 , y(1) = 1.2

Rezolvare:Pas1: Se transforma ecuatia diferentiala de ordinul doi intr-un sistem de doua ecuatii diferentiale deordinul unu. Se notam cu z=y sistemul normal asociat ecuatiei diferentiale de ordinul doi (ProblemaCauchy) se scrie sub formay=z , y(0)=1.25z=0.5*y-2*z^2/y, z(0)=se alege conform Pas1 din algoritmul prezentat mai susPas2: Se rezolva problema Cauchy pe intervalul [0,4], folosind functia Octave lsodeoctave#1>function zdot=ftir(y,x)octave#2> zdot(1)=y(2); zdot(2)= 0.5*y(1)-(2*y(2). 2)./y(1); endfunctionoctave#3>x=0:0.1:1 ;# Se alege y(0)=1 si y(0)=0 si se obtineoctave#4> [Y,istate,msg]=lsode("ftir",[1; 0],x) # se obtine y(1)=1.227283506205422octave#5># Daca se alege y(0)=1 si y(0)=-0.05octave#6>[Y,istate,msg]=lsode("ftir",[1; -0.05],x) # se obtine y(1)=1.18360810910837Ce valoare trebuie sa luam pentru y(0) la urmatoarea incercare ?

7 . 4 . 2 . METODA DIFERENTELOR FINITE CENTRATEFie ecuatia diferentiala de ordinul doi (12.1) adusa la forma

( ) ( ) ( )xRyxQx

yxP

x

y=++

d

d

d

d2

2

(12.3)

cu conditiile la limita (12.2)

Se cere sa afle solutia aproximativa a problemei bilocale pentru [ ]bax , Se discretizeaza intervalul [a=x0 ;b= xn] ntr-un numr de intervale egale, distanate cu pasul

( ) nxxh n 0= , cele (n+1) valori discrete fiind hkxxk 0 += , cu nk ,0= .


2/6

Daca pentru derivatele y(x) si y(x) facem aproximarile obtinute prin diferente finite centrate:

ii

kkk

k

kk

k yxycuh

xyxyxyxy

h

xyxyxy =

+=

= ++ )(,

)()(2)()(",

2

)()()('

2

1111(12.4)

ecuatia (12.3) construita cu diferene finite centrate se scrie :

( ) ( ) ( kkkkk

k

kkk xRyxQ

h

yyxP

h

yyy=+

+

+ ++

2

2 11

2

11 ), cu 1,1 = nk ,yk=y(xk) (12.5)

Ecuatiile in diferente finite (12.5) se scriu sub forma:

( )[ ] ( ) ( )[ ] ( )kkkkkkk xRhxPhyxQhyxPhy

2

1

2

1 22242 =++++ + , 1,1 = nk (12.6)In transcriere matricila, sistemul liniar cu necunoscutele y1, y2, , yn-1 (12.6) se scrie sub forma A*y=funde A este o matrice tribanda de forma

+

++

++

=

)(24)(2....000

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

0...0)(2)(24)(20

0...0)(2)(24)(2

1

2

1

33

2

3

22

2

2

nnxQhxhP

xhPxQhxhP

xhPxQhxhP

A

si y=[ y1, y2, , yn-1]T, f=[2h2R(x1)-y0(2-hP(x1), 2h2R(x2), . . . , 2h2R(xn-2), 2h2R(xn-1)-yn(2+hP(xn-1)]T

Ecuatiile difereniale ordinare neliniare se construiesc, in mod asemanator, cu diferente finite centrate si se

rezolva pentru nodurile interioare domeniului de calcul ( 1,1 = nk ).

Exemplul 2. Sa se afle solutia aproximativa a problemei bilocale prin metoda diferentelor finite peintervalul x[0,1] cu pas h=0.2(1+x2)y-xy-3y=6x-3 cu conditiile la limita y(0)=3, y(1)=2

Rezolvare : Pentru P(x)=-x/(1+x2), Q(x)=-3/(1+x2), R(x)=(6x-3)/( 1+x2) sistemul (12.6) se scrie subforma

4,3,2,1),1/()3x*(6*2.0*2

)1/()(-x*2.02)1/()(-3*2.0*24)1/()(-x*0.22

2

k

2

2

k1

222

k1

=+=

+++++++ +

kx

xyxyxy

k

kkkkkk

(12.7)Daca se considera tabelul care contine diviziunea intervalului [0 ;1] de pas h=0.2si cunoscandu-se valorilela limita y0=y(x0)=3 si y5=y(x5)=2, se cer valorile intermediare y1, y2, y3, y4xk x0=0.00000 x1=0.20000 x2=0.40000 x3=0.60000 x4=0.80000 x5=1.00000yk y0=y(x0)=3 y1=? y2=? y3=? y4=? y5=y(x5)=2

Din (12.7) se obtine sistemul A*Y=B unde A =[ 2.1154 -0.98077 0.00000 0.00000; -1.0345 2.1034 -0.96552 0.00000

0.00000 -1.0441 2.0882 -0.95588; 0.00000 0.00000 -1.0488 2.0732]B=[3.126923 0.020690 -0.017647 1.858537]T, Y=[ y1 y2 y3 y4]

T

octave:33> A\B # solutia esteSolutie : y1=2.4351 y2=2.0639 y3=1.8660 y4=1.8404Observatie : Metoda diferentelor finite se poate aplica si daca problema bilocala este complet aca cea din(12.2)Exemplul 3. Sa se afle solutia aproximativa a problemei bilocale prin metoda diferentelor finite peintervalul x[0,1] cu pas h=0.2(1+x2)y-xy-3y=6x-3 cu conditiile la limita y(0)-y(1)=3, y(1)=2Rezolvare : Se aplica metoda diferentelor finite pentru pasul h=0.2. Necunoscutele sunt y1, y2, y3, y4obtinute pentru xn=0.2*n, n=1,2,3,4.contine diviziunea intervalului [0 ;1] de pas h=0.2si cunoscandu-sevalorile la limita y0=y(x0)=3 si y5=y(x5)=2, se cer valorile intermediare y1, y2, y3, y4 dar si y0Conform relatiei (12.5) , prima conditie la limita y(0)-y(1)=3 se scrie sub forma

14.0

11

1 =

yy

y (12.8)


3/6

Daca facem in ecuatia (1+x2)y-xy-3y=6x-3 , x=0 se obtine y-3y=-3 si conform (12.5) aceasta relatie se

scrie sub forma 3304.0

20

101 =+ yyyy

(12.9)

Se elimina y-1 dintre (12.9) si (12.10) si se obtine prima ecuatie a sistemului-2.52*y0 +2*y1=-0.52. In continuare pentru n=1,2,3,4 din (12.8) se obtin urmatoarele 4 ecuatii alesistemului. In final pentru necunoscutele y0, y1, y2, y3, y4 se obtine sistemul A*Y=B unde A=[-2.52 2 0 0 0 ;1.06 -2.20 1.02 0 0 ;0 1.2 -2.44 1.12 0 ;0 0 1.42 -2.84 1.2 ; 0 0 0 1.72 -3.4],B=[-0.52 -0.072 -0.024 0.024 -3.048]T, Y=[ y0 y1 y2 y3 y4]

T. Din rezolvarea sistemului (spre exempluprin metoda Gauss) se obtine solutia : y0=1.0132, y1=1.0167, y2=1.0693, y3=1.2188, y4=1.513Observatie : Pentru imbunatatirea solutiei aproximative obtinute se reia algoritmul prezentat pentru unpas mai fin pe intervalul [0 ;1], spre exemplu h=0.1

7 . 4 . 3 . METODA CU FUNTII SPLINE

Fie problema bilocala ( ) ( ) ( )xRyxQx

yxP

x

y=++

d

d

d

d2

2

, pentru x [a,b]

(12.10)

=+

=+

221

121

)(')(

)(')(

bybbyb

ayaaya

Consideram o diviziune oarecare pentru intervalul [a;b], bxxxxa nn =


4/6

{ }

{ } 2])1()1([21)1(

])1()1([6

1)1()1(

1

)()((

*215*2

04032

*3

15*3

04

2

03021

21

*15*043

=+++

+++++

=+

=++

xKxKxKK

xKxKxKxKK

KK

xxxKxxKK

(12.15)

Daca in sistemul (12.15) se face x=xi, i=0,1,2 se obtine sistemul de 5 ecuatii cu 5 necunoscute A*Y=BUnde A=[0 0 1 0 0 ;0 0 1 0.5 0 ;0 0 1 1 0.5 ;1 1 0 0 0 ;1 2 2 2/3 2/3], Y=[ K1 K2 K3 K4 K5]

T ,B=[0 0.5 1 1 2]T. Din rezolvarea sistemului (spre exemplu prin metoda Gauss) se obtinK1=0.66667, K2=0.33333, K3= 0.00000, K4=1.00000, K5=0.00000 si solutia se scrie sub formay(x)= 0.66667+0.33333*x+(1/6)*x 3

7 . 5 . INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE ORDINARE DE ORDINMA MARE DECAT II, CU CONDITII LA LIMITA

7 . 5 . 1 . METODA DIFERENTELOR FINTE CENTRATEDaca pentru derivatele y(x) ,y(4)(x), . . . facem aproximarile obtinute prin diferente finite centrate:

ii

kkkkk

k

kkkk

k

yxycuh

xyxyxyxyxyxy

h

xyxyxyxy

xy

=++

=

+

=

++

++

)(...,,)()(4)(6)(4)(

)(

,2

)()(2)(2)(

)('''

4

2112)4(

3

2112

(12.14)

atunci metoda diferentelor finite centrate poate fi aplicata pentru o ecuatie diferentiala cu conditii la limitade orice ordin n , numar natural (pentru aproximarea derivatelor y si y se folosesc formulele (12.4) )Exemplul 5. Sa se afle solutia aproximativa a problemei bilocale prin metoda diferentelor finite centrate peintervalul x[0,1] cu pas h=1/4.y(4)+6*y=100 cu conditiile la limita y(0)=y(0)=0, y(1)=y(1)=0Rezolvare : Folosind formulele (12.14) si tabelul

xk x0=0 x1=1/4 x2=1/2 x3=3/4 x4=1yk y0= 0 y1= ? y2= ? y3= ? y4=0

ecuatia data se scrie in diferente centrate sub forma

3,2,1,100)()()(4)(6)(4)(

4

2112 ==+++ ++ kxy

h

xyxyxyxyxyk

kkkkk(12.15)

cu conditiile la limita : y(x0)= 0, y(x-1)= y(x1), y(x4)= 0, y(x5)=- y(x3) (12.16)Daca notam y(xk)=yk, atunci pentru k=1,2,3 sistemul (12.15) impreuna cu conditiile la limita (12.16) sescrie

43421

43241

43214

4

100]1)

4

11(6[4

4

1004)

4

11(64

4

1004]1)

4

11(6[

=++

=++

=+++

yyy

yyy

yyy

(12.17)

Se obtini solutiile : y1=0.34438, y2=0.63592, y3=0.51557.


5/6

APLICATII

Problema 1. Program Octave pentru rezolvarea problemei bilocaley = p(x)y`(x)+q(x)y(x)+r(x), y(a) = alpha, y(b) = beta

Algoritm:1. Informatia initiala : Citeste functiile p,q,r; citeste extremitatile intervalului de integrare a si b ; citeste

conditiile bilocale alpha ,beta si numarul n de subintervale de integrare2. Initializeaza vectorul X; vectorul solutie Y si vectorii Va (supradiagonal), Vd (subdiagonal), Vc

(diagonal) , Vb (termeni liberi);defineste pasul de integrare h(b - a)/n3. Pentru j=1:n-1, calculeaza vectorul nodurilor de integrareVx(j) a + h*j , vectorul termeni liberi Vb(j)

-h^2*feval("r",Vx(j)) si apoi completeaza Vb(j) cu conditiile bilocale (la capete)4. Pentru j=1:n-1, calculeaza vectorul diagonal:Vd(j) 2 + h^2*feval("q",Vx(j))5. Pentru j=1:n-2, calculeaza vectorii subdiagonal si supradiagonal:

Va(j) -1 - h/2*feval("p",Vx(j+1)) , Vc(j) = -1 + h/2*feval("p",Vx(j));6. Rezolva sistemul tridiagonal A*Y=B si afiseaza vectorul nodurilor X si vectorul solutie Y, calculate in

fiecare element al lui X

octave#1> function [X,Y] = edobilocal(p,q,r,a,b,alpha,beta,n)X = zeros(1,n+1);Y = zeros(1,n-1);Va = zeros(1,n-2);Vb = zeros(1,n-1);Vc = zeros(1,n-2);Vd = zeros(1,n-1); h = (b - a)/n;# Se calculeaza matricea tridiagonala si vectorul termini liberifor j=1:n-1,

Vx(j) = a + h*j;endforfor j=1:n-1,Vb(j) = -h^2*feval("r",Vx(j));

endforVb(1) = Vb(1) + (1 + h/2*feval("p",Vx(1)))*alpha;Vb(n-1) = Vb(n-1) + (1 - h/2*feval("p",Vx(n-1)))*beta;for j=1:n-1,Vd(j) = 2 + h^2*feval("q",Vx(j));

endforfor j=1:n-2,Va(j) = -1 - h/2*feval("p",Vx(j+1));

endforfor j=1:n-2,Vc(j) = -1 + h/2*feval("p",Vx(j));

endfor

# Se rezolva sistemul tridiagonaln = length(Vb);for k = 2:n,

m = Va(k-1)/Vd(k-1);Vd(k) = Vd(k) - m*Vc(k-1);Vb(k) = Vb(k) - m*Vb(k-1);

endforY(n) = Vb(n)/Vd(n);for k = (n-1):-1:1,Y(k) = (Vb(k) - Vc(k)*Y(k+1))/Vd(k);

endforX = [a,Vx,b];Y= [alpha,Y,beta];endfunction

Problema 2. S se determine solutia problemei bilocale de mai jos pe intervalul x [0 ;1], cu metodadiferentelor finite centrate, considerand pasul 2,0=h :

( )( )

=

=

=+

11

00

024

y

y

xyyy

.

Comparati solutia obtinuta cu solutia numerica obtinuta folosind in mediul de programare Octave functiaodebilocal pentru p(x)=-4, q(x)=2*x, r(x)=0, a=0, b=1, n=4,alpha=0, beta=1


6/6

Problema 3.S se determine solutia problemei bilocale de mai jos pe intervalul x [0 ;3], cu metodadiferentelor finite centrate, considerand pasul 75,0=h :

( )( )

=

=

=

13

10

1

y

y

yy

.

Comparati solutia obtinuta cu solutia numerica obtinuta folosind in mediul de programare Octave functiaodebilocal pentru p(x)=0, q(x)=1, r(x)=1, a=0, b=3, n=3,alpha=-1, beta=-1

Problema 4.Fie ecuatia diferentiala unidimensionala a radiatiei (ecuatie neliniara) pe intervalul x[0 ;1] :

( )( )

=

=

=

01

00

0e

y

y

yy

, a carei forma in diferente finite centrate este : 0e2

2

11 =+ + kykkk

h

yyy.

Se cere sa se determine solutia aproximativa considerand 3 intervale egale pe domeniul x [0, 1]. Sistemulde ecuatii neliniare rezultat se va rezolva iterativ cu metoda Newton-Raphson, pornind de la aproximatia

initiala .0)0(

2)0(

1== yy

Sa se compare solutia obtinuta cu solutia numerica obtinuta prin metoda tirului (pentru problemele Cauchyobtinute folositi functia Octave lsode )

Problema 5.Fie ecuatia y=-2*y*y cu conditia bilocala y(0)+y(0)=0, y(1)=0.5. Consideram diviziuneapentru x[0;1] de pas h=0.2 sa se scrie sistemul neliniar in yk, k=0,1,2,3,4,5 atasat problemei bilocalefolosind metoda diferentelor finite centrate.Indicatie: Pentru derivatele y(x) si y(x) se folosesc aproximarile (12.4). Sistemul atasat se scrie

5,4,3,2,1,0,04.0

)()()(2

04.0

2 1111 ==

++ ++ k

xyxyxy

yyy kkk

kkkcu conditiile la limita

y5=0.5 iar pentru prima ecuatie a sistemului se elimina y-1 din relatiile

04.0

04.0

204.0

2

11

0

11

0

211

=

+

=

++

yyy

yyy

yyy

Problema 6.Sa se rezolve problema bilocala de ordin trei, pentru x [0, 1], prin metoda diferentelorfinite centrate si considerand o diviziune a intervalului x[0;1] de pas h=0.25y=4(y+x) cu conditiile la limita y(0)=1, y(1)=2

INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE ORDINARE CU CONDITII LA LIMITA (Probleme...

Documents

Transcript of INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE ORDINARE CU CONDITII LA LIMITA (Probleme...