INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE DE ORDINUL I, CU CONDITIE INITIALE (Probleme Cauchy)

8/7/2019 INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE DE ORDINUL I, CU CONDITIE INITIALE (Probleme Cauchy)

1/5

Carmen-Sanda Georgescu, Tudor Petrovici, Radu PopaMetode numerice. Fisa de laborator nr. 10:

INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE DEORDINUL I, CU CONDITIE INITIALE (Probleme Cauchy)

7. INTEGRAREA NUMERIC A ECUAIILOR DIFERENIALE ORDINARE7 . 1 . M E TO D E D I R E CTE D E I N TE G R A R E A E CU A TI E I D I F E R E N TI A LE

D E O R D I N U L I N TA I , CU CO N D I TI E I N I TI A LA

Fie ecuaia diferenial ordinar de ordinul nti

( yxfy )x

y,

d

d== . (10.1)

Se caut soluia pe intervalul [x( )xy 0, xf], care s satisfac ecuaia (10.1), cu condiia iniial:( ) 00 yxy = . (10.2)

Ecuatia diferentiala (10.1) impreuna cu conditia initiala (10.2) formeaza o problema Cauchy.

Fiex[x0, xf] unde xf=xn . Se mparte intervalul [x0, xn] n n subintervale egale, distanate cu pasul

n

abh

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

xk x0 x1 ... xnyk y0 y1 ... yn

7.1.1. METODE DE TIP RUNGE-KUTTA

RK1 = METODA RUNGE-KUTTA DE ORDINUL INTAI ( METODA EULER)Solutia ecuaei difereniale de ordinul nti (10.1) cu conditia initiala (10.2) se obtine prin schema iterativa

( kkkk yxfhyy ,1 +=+ ) k=1,2,3, . . . si y(xo)=dat (10.3)Exemplul 1: Sa se rezolve problema Cauchy y=(x-y)/2, y(0)=1 pe intervalul [0,1] cu pasul h=0.2 folosindmetoda lui Euler. Notam K1=h*f(xk,yk)k xk y0=1, yik+1= yk+K1 K1=h*f(xk,yk)0 0 y0=1 K1 = -0.100001 0.2 y1=0.90000 K1 = -0.070002 0.4 y2=0.83000 K1 = -0.0430003 0.6 y3=0.78700 K1 =-0.0187004 0.8 y4=0.76830 K1 = 0.00317005 1 y5=0.77147Eroarea comisa este de ordinul lui h2 adica 0.04Program Octave pentru rezolvarea unei ecuatii diferentiale de ordinul I prin metoda lui Euler pe intervalulx[a,b] cu pas h=(b-a)/m si y(a)=y0

octave#1> function [X,Y] = euleredo(f,a,b,y0,m)h = (b - a)/m;

X = zeros(1,m+1);X = zeros(1,m+1);X(1) = a;Y(1) = y0;for j=1:m,Y(j+1) = Y(j) + h*feval(f,X(j),Y(j))X(j+1) = a + h*jEndforendfunction

octave#1>function z = fedo(x,y) ,z = (x-y)/2;endfunction


2/5

octave#1>a=0,b=1,y0=1,m=100, [X,Y]=euleredo("fedo",a,b,y0,m)Rezultat: y(1)=0.81731; Eroarea comisa este de ordinul lui 0.0001

RK2 = METODA RUNGE-KUTTA DE ORDINUL AL DOILEASolutia ecuaei difereniale de ordinul nti (10.1) cu conditia initiala (10.2) se obtine prin schema iterativa

+++=

+ kkkkk y

h

y

h

xfhyy 2,

2

1unde

( )kkk yxfy,= , k=1,2,3, . . . si y(x

o)=dat (10.4)

Exemplul 2: Sa se rezolve problema Cauchy y=(x-y)/2, y(0)=1 pe intervalul [0,1] cu pasul h=0.2 folosindmetoda Runge-Kutta de ordinul II. Notam K1=h*f(xi,yi) si K2=h*f(xk+h/2,yk+K1/2). Atunci solutia seobtine cu schema iterativa y0=1, yk+1= yk+K2k Xk y0=1, yk+1= yk+K2 K1=h*f(xi,yi), K2=h*f(xk+h/2,yk+K1/2)0 0 y0=1 K1 = -0.10000, K2 = -0.0850001 0.2 y1=0.91500 K1 = -0.071500, K2 = -0.0579252 0.4 y2=0.85708 K1 = -0.045708, K2 = -0.0334223 0.6 y3= 0.82365 K1 = -0.022365, K2 = -0.0112474 0.8 y4= 0.81241 K1 = -0.0012406, K2 = 0.00882145 1 y5=0.82123Program Octave pentru rezolvarea unei ecuatii diferentiale de ordinul I prin metodaRunge-Kutta II, peintervalul x[a,b] cu pas h=(b-a)/m si y(a)=y0

octave#1>function z = fedo(x,y) ,z = (x-y)/2;endfunctionoctave#1>a=0,b=1,y0=1,m=100, [X,Y]=rk2("fedo",a,b,y0,m)Rezultat: y(1)=0.81960

RK3 = METODA RUNGE-KUTTA DE ORDINUL AL TREILEASolutia ecuaei difereniale de ordinul nti (10.1) cu conditia initiala (10.2) se obtine prin schema iterativa

4

3

4

311

KKyy kk ++=+ (10.5)

unde

( )kk yxfhK ,1 =

++=

3,

3 12

Ky

hxfhK kk (10.6)

++=

3

2,

3

2 23

Ky

hxfhK kk , k=1,2,3, . . . si y(xo)=dat

Exemplul3: Sa se rezolve problema Cauchy y=(x-y)/2, y(0)=1 pe intervalul [0,1] cu pasul h=0.2 folosindmetoda Runge-Kutta de ordinul III. Notam K1=h*f(xi,yi) si K2=h*f(xk+h/2,yk+K1/2) siK3=h*f(xk+h/2,yk+K2/2) Atunci solutia se obtine cu schema iterativa y0=1, yk+1= yk+ K1/4+3*K3/4k Xk y0=1, yk+1= yk+ K1/4 + 3*K3/4 K1=h*f(xi,yi), K2=h*f(xk+h/2,yk+K1/2),K3=h*f(xk+h/2,yk+K2/2)0 0 y0=1 K1 = -0.10000, K2 = -0.085000, K3 = -0.0857501 0.2 y1=0.91069 K1 = -0.071069, K2 = -0.057515, K3 = -0.0581932 0.4 y2=0.84928 K1 = -0.044928, K2 = -0.032681, K3 = -0.0332933 0.6 y3= 0.81307 K1 = -0.021307, K2 = -0.010242, K3 = -0.0107954 0.8 y4= 0.79965 K1 = 3.4972e-05, K2 = 0.010033, K3 = 0.00953335 1 y5= 0.80681

octave#1> function [X,Y] = rk2(f,a,b,ya,m)

h = (b - a)/m; X = zeros(1,m+1); Y = zeros(1,m+1); X(1) = a; Y(1) = ya;for j=1:m,xj = X(j);yj = Y(j);k1 = h*feval(f,xj,yj)k2 = h*feval(f,xj+h/2,yj+k1/2)Y(j+1) = yj + k2X(j+1) = a + h*jendforendfunction


3/5

Program Octave pentru rezolvarea unei ecuatii diferentiale de ordinul I prin metoda Runge-Kutta III, peintervalul x[a,b] cu pas h=(b-a)/m si y(a)=y0

octave#1> function [X,Y] = rk3(f,a,b,ya,m)h = (b - a)/m;X = zeros(1,m+1);Y = zeros(1,m+1);X(1) = a;Y(1) = ya;for j=1:m,xj = X(j);yj = Y(j);

k1 = h*feval(f,xj,yj);k2 = h*feval(f,xj+h/2,yj+k1/2)k3 = h*feval(f,xj+h/2,yj+k2/2)Y(j+1) = yj + k1/4 +3*k3/4X(j+1) = a + h*jendforendfunction


RK4 = METODA RUNGE-KUTTA DE ORDINUL PATRU (cel mai frecvent utilizat n literatura despecialitate)

Solutia ecuaei difereniale de ordinul nti (10.1) cu conditia initiala (10.2) se obtine prin schema iterativa

( )43211 226

1KKKKyy kk ++++=+ (10.7)

unde

( )kk yxfhK ,1 =

++=

2,

2 12

Ky

hxfhK kk

++=

2,

2 23

Ky

hxfhK kk (10.8)

( 34 , KyhxfhK kk ++= ) , k=1,2,3, . . . si y(xo)=datExemplul 4: Sa se rezolve problema Cauchy y=(x-y)/2, y(0)=1 pe intervalul [0,1] cu pasul h=0.2 folosindmetoda Runge-Kutta de ordinul IV. Notam K1=h*f(xi,yi) si K2=h*f(xk+h/2,yk+K1/2),K3=h*f(xk+h/2,yk+K2/2) si K4=h*f(xk+h,yk+K3). Atunci solutia se obtine cu schema iterativa y0=1, yk+1=yk+(1/6)*(K1+2*K2+2*K3+K4)k Xk y0=1, yk+1=

(1/6)*(K1+2*K2+2*K3+K4)

K1=h*f(xi,yi), K2=h*f(xk+h/2,yk+K1/2),K3=h*f(xk+h/2,yk+K2/2)

0 0 y0=1 K1 = -0.10, K2 = -0.0850, K3 = -0.085750, K4 = -0.0714251 0.2 y1=0.91451 K1 = -0.07145, K2 = -0.057879, K3 = -0.058557,K4 = -0.045592 0.4 y2=0.85619 K1 = -0.04562, K2 = -0.03333, K3 = -0.03395, K4 = -0.0222243 0.6 y3= 0.82246 K1 = -0.02224,K2 = -0.01113,K3 = -0.01168,K4 = -0.00107664 0.8 y4= 0.81096 K1 = -0.001096,K2 = 0.008958,K3 = 0.008456,K4 = 0.0180585 1 y5= 0.81959Program Octave pentru rezolvarea unei ecuatii diferentiale de ordinul I prin metoda Runge-Kutta III, peintervalul x[a,b] cu pas h=(b-a)/m si y(a)=y0

octave#1> function [X,Y] = rk4(f,a,b,ya,m)h = (b - a)/m;X = zeros(1,m+1);Y = zeros(1,m+1);X(1) = a;Y(1) = ya;for j=1:m,xj = X(j);yj = Y(j);k1 = h*feval(f,xj,yj);k2 = h*feval(f,xj+h/2,yj+k1/2)k3 = h*feval(f,xj+h/2,yj+k2/2)k4 = h*feval(f,xj+h,yj+k3)

Y(j+1) = yj + (k1 + 2*k2 + 2*k3 + k4)/6X(j+1) = a + h*jendforendfunction


4/5


Exemplul 5. Sa se rezolve problema Cauchy pe intervalul x[0,1] cu pas h=0.1y=y-2*x/y , y(0)=1folosind metodele lui Euler si Runge-Kutta de ordin 4 (RK4)si sa se compare grafic rezultatele obtinute cusolutia obtinuta daca se foloseste functia Octavelsode si solutia teoretica y=sqrt(2*x+1)X =

0.00000 0.10000 0.20000 0.30000 0.40000 0.50000 0.60000 0.70000 0.80000 0.90000 1.00000Yteoretic=1.00000 1.09545 1.18322 1.26491 1.34164 1.41421 1.48324 1.54919 1.61245 1.67332 1.73205

YEuler =1.00000 1.10000 1.19182 1.27744 1.35821 1.43513 1.50897 1.58034 1.64978 1.71778 1.78477

YRK4=1.0 1.09545 1.18322 1.26491 1.34164 1.41422 1.48324 1.54920 1.61246 1.67332 1.73206

Ylsode =1.0 1.09545 1.18322 1.26491 1.34164 1.41421 1.48324 1.54919 1.61245 1.67332 1.73205

Erori relative ale solutiei la capatul intervalului de integrare ( in x=1) :erEuler=(Yteoretic(1)-Yeuler(1))/ Yteoretic(1)*100%=-3.0438 %erRK4=(Yteoretic(1)-YRK4(1))/ Yteoretic(1)*100%= -3.2087e-04 %erLSODE=(Ylsode(1)-Ylsode(1))/ Yteoretic(1)*100%=2.3101e-06 %


5/5

Functii Octave pentru rezolvarea ecuatiilor diferentiale cu conditii initiale

octave#1> [x,istate,msg]= lsode ("f", x0, t,{t_cr});

Rezolva o ecuatie/sistem de ecuatii, diferentiale ordinare cu conditii initiale. Intrari :f - este numele functiei care descrie ecuatia/sistemul

x0 este vectorul conditiilor initialet este intervalul in care se rezolva ecuatiat_cr este optionalIesiri : x vectorul solutie ; istate scalar de validare a integrarii ; msg mesaj text

octave#1> [x,istate,msg]= lsode (["f" "j"], x0, t,{t_cr});

Rezolva o ecuatie/sistem de ecuatii, diferentiale ordinare cu conditii initiale, folosindJacobianul functiei f. In acest caz functia lsode trebuie sa fie acompaniata de functialsode_options (opt, val) (help lsode_options (opt, val) )

Alte functii pentru rezolvarea ecuatiilor/sistemelor diferentiale ordinare cu conditii initiale (Problema Cauchy)octave#1> [tout,xout]=ode23(FUN,tspan,x0,{alte_argumente_optionale})octave#1> [tout,xout]=ode45(FUN,tspan,x0,{alte_argumente_optionale})octave#1> [tout,xout]=ode78(FUN,tspan,x0,{alte_argumente_optionale})Rezolva o ecuatie/sistem de ecuatii, diferentiale ordinare cu conditii initiale

Intrari: FUN - este un string care contine numele functiei care descrie ecuatia diferentiala tspan- esteintervalul de integrarex0- vectorul conditiilor initiale{alte_argumente_optionale} pair, ode_fcn_format,tol,trace,count,hmax... (help ode..)Iesiri : xout vectorul solutie ; tout vectorul variabilei independente

APLICAII

S se integreze numeric ecuaiile difereniale ordinare, cu condiie iniial, de mai jos, plecnd de lax0pn la valoarea finalxn, utiliznd metodele Runge-Kutta. Comparai rezultatele obinute cu RK1,..., RK4.In fiecare caz in parte sa se determine solutia utilizand functiile Octave: lsode, ode23, ode45, ode78

Problema 1.y

xy = , cu i( ) 42 =y 0=nx

Problema 2. , cu iyy = ( ) 10 =y 2=nx

Problema 3.x

yy

e

2

= , cu i( ) 10 =y 1=nx

S se integreze numeric ecuaiile difereniale ordinare, cu condiie iniial, de mai jos, plecnd de lavaloarea iniial, pn la valoarea final, utiliznd metoda Runge-Kutta de ordinul 4.Sa se determine de asemenea solutia utilizand functiile Octave: lsode si ode23

Problema 4.y

xyy

2= , cu i( ) 10 =y 1=nx

INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE DE ORDINUL I, CU CONDITIE INITIALE (Probleme Cauchy)

Documents

Transcript of INTEGRAREA NUMERICA A ECUATIILOR DIFERENTIALE DE ORDINUL I, CU CONDITIE INITIALE (Probleme Cauchy)