Carmen Georgescu ,Tudor Petrovici

Catedra de Hidraulica si Masini Hidraulice, Universitatea Politehnica din Bucuresti

Laborator 12. Rezolvarea numerica a ecuatiilor si sistemelor de ecuatii

diferentiale ordinare cu conditii initiale (probleme Cauchy)

1.Ecuatii diferentiale ordinare de ordinul 1

O ecuatie diferentiala de ordinul I, este o relatie de tip F(x,y,y)=0, unde y=y(x) si y(x)=dy/dx. Daca se adauga si conditia initiala y(x0)=y0 , atunci solutia y(x) este perfect determinata (problema Cauchy). Daca x si y apartin lui Rn, atunci avem un sistem de ecuatii diferentiale de ordin I. In cazul general o problema Cauchy de ordinul n, este o relatie de tip F(x,y,y,y, . . . ,y(n))=0, unde y=y(x) , y(x)=dy/dx, . . . , y(n)=dny/dxn la care se adauga conditiile initiale y(x0)=y0 , . . . , y

(n)(x0)=yn Abordare numerica in OCTAVE

Aplicatia #1 Sa se rezolve problema Cauchy: y' + 2xy =0 unde y=y(x) y(0)=3. Sa se scrie solutia folosind 20 de puncte intermediare intre 0 si 1 Solutie teoretica:

y'=- 2xy dy/dx=- 2xy dy/y=-2xdx ln(y)=-x2+C y(x)=Cexp(-x2) C=3 y(x)=3exp(-x2)

Solutie numerica: Pas 1: Se scrie functia care descrie ecuatia diferentiala >>function dy=aplicatie1(y,x) >dy=-2*x*y; >endfunction Pas 2 : >>x=linspace(0,1,20);[y,istate,msg]= lsode ("aplicatie1",3,x); >>y1=3*exp(-x.^2);plot(x,y,x,y1),legend(Solutie teoretica,Solutie lsode),grid

De retinut : Functia Octave: lsode Sintaxa: [y,istate,msg]= lsode ("f ", y0, x);

Aplicatia #2 Sa se rezolve problema Cauchy: y-2y=-x2 unde y=y(x) y(0)=0.25 Sa se scrie solutia folosind 22 de puncte intermediare intre 0 si 1 Solutie teoretica: y(x)=0.25(2*x^2+2*x+1) Solutie numerica: Pas 1: Se scrie functia care descrie ecuatia diferentiala >>function dy=aplicatie2(y,x) >dy=2*y-x^2; >endfunction Pas 2 : >>x=linspace(0,1,22);[y,istate,msg]= lsode ("aplicatie2", 0.25, x); >>y1=0.25*(2*x.^2+2*x+1);plot(x,y,x,y1) >> legend(Solutie teoretica,Solutie numerica),grid

2.Sisteme de ecuatii diferentiale ordinare de ordinul 1

Un sistem de ecuatii diferentiale de ordinul I, este este un sistem de forma dx1/dt=f1(t,x1,x2, ... ,xn) dx2/dt=f2(t,x1,x2, ... ,xn) .......................... dxn/dt=fn(t,x1,x2, ... ,xn) Daca se aduga conditiile initiale x1(t0)=x10 , x2(t0)=x20 , ..., xn(t0)=xn0 se defineste o problema Cauchy Aplicatia #3 Atractorul lui Lorentz Sa se rezolve sistemul de ecuatii diferentiale x=a(y-x) y=x(b-z)-y z=xy-cz intre t=0 si t=50 cu conditiile initiale x(0)=y(0)=z(0)=1 (a=10, b=28, c=8/3 ) si x=x(t), y=y(t), z=z(t)

Sa se foloseasca functia lsode cu 5000 de pasi de integrare. Sa se reprezinte grafic solutia folosind functia plot3 Pas_1: Se editeaza in linia de comanda functia numita spre exemplu func, care descrie sistemul de ecuatii diferentiale si notam: x=[x(1) x(2) x(3)], xdot=[xdot(1) xdot(2) xdot(3)] >>function xdot=func(x,t) >a=10;b=28;c=8/3; >xdot(1)=a*(x(2)-x(1)); >xdot(2)=x(1).*(b-x(3))-x(2); >xdot(3)=x(1).*x(2)-c*x(3); >endfunction Pas_2: Program >>x0=[1;1;1];t=linspace(0,50,5000);x=lsode("func",x0,t); >> plot3 (t,x(:,1),t,x(:,2),t,x(:,3)) >>plot3 (x(:,1),x(:,2),x(:,3))

Aplicatia #4. Problema Vanator-Prada Sa se rezolve problema Cauchy (sistem de ecuatii diferentiale cu conditii initiale) dy1/dt=y1(r1-a1y2) dy2/dt=y2(a2y1- r2) pe intervalul [0,5] cu 10 pasi intermediari ( r1, r2 sunt ratele de crestere ale speciilor y1, y2 iar a1, a2 sunt ratele de mortalitate) Program: Se considera ca a1= a2=1 si r1= r2=0.5 >>function xdot = vp(y, x) >xdot(1)= y(1)*(0.5-y(2)*1);

>xdot(2)=y(2)*(-0.5+y(1)*1); >endfunction >>x0=[0.1;0.1];t=0:500; >>x=lsode("vp",x0,t) ; >>plot(t,x(:,1),t,x(:,2)) >> plot(x(:,1),x(:,2)) % Traiectoria

3.Eecuatii diferentiale ordinare de ordinul superior

Orice ecuatie diferentiala de ordin mai mare decat unu, se transforma mai intai intr-un sistem de ecuatii diferentiale de ordinul unu (numit sistemul ecuatiilor normale) si se rezolva ca in paragraful 2

Exemplu: Aplicatia #5 Ecuatia Van der Pol Sa se rezolve ecuatia diferentiala y'' + (y2-1)y+y=0 cu y=y(x), pe intervalul [0;20] cu conditiile initiale y(0)=0.25 si y'(0)=0 . Sa se reprezinte grafic solutia numerica Pas1 : Se transforma ecuatia diferentiala de ordinul doi intr-un sistem de ecuatii diferentiale de ordinul unu, notand y2=y1 , rezulta sistemul (sistemul de ecuatii diferentiale de ordin 1, echivalent cu ecuatia diferentiala de ordin superior, se numeste sistemul ecuatiilor normale): y1=y1(1-y2

2)-y2 y2=y1 Pas2: Se editeaza in linia de comanda functia numita ecdif, care descrie sistemul de ecuatii normale >>function yprim = ecdif(y,x) >yprim(1)=y(1).*(1-y(2).^2)-y(2) ; >yprim(2)=y(1) ; >endfunction Pas3: Program >>x=0:0.1:20;y=lsode("ecdif",[0.25;0] ,x) >> plot(x,y(:,1))

Probleme propuse

1. Se se rezolve numeric ecuatia diferentiala cu conditie initiala y=3*x^3-2, y(1)=1 pe intervalul [1 ,3] .Sa se compare solutia numerica obtinuta cu solutia teoretica. Sa se reprezinte pe acelasi grafic solutia numerica si solutia teoretica 2. Se se rezolve numeric ecuatia diferentiala cu conditie initiala : v+2v+3 v-4 = 0 , v(0) = 12, v(0)=0 pe intervalul [0,10] . Sa se reprezinte grafic v=v(t) pe intervalul [0,10] 3. Sa se rezolve numeric sistemul omogen de ecuatii diferentiale, cu conditii initiale si sa se compare cu solutia teoretica x=x+y ; x(0)=0 y=-2*x+4*y ; y(0)=-1 HELP : Solutia teoretica x(t)=exp(2*t)-exp(3*t), y(t)= exp(2*t)- 2*exp(3*t) 5. Sa se rezolve numeric sistemul omogen de ecuatii diferentiale, cu conditii initiale si sa se compare cu solutia teoretica x=x+y+t-t^2-2 ; x(0)=0 y=-2*x+4*y+2*t^2-4*t-7 ; y(0)=2 HELP : Solutia teoretica x(t)=t^2, y(t)= t+2