Rezolvarea ecuatiilor diferentiale

Aceasta lucrare descrie modul de rezolvare numerica in MATLAB a ecuatiilor diferentiale. MATLAB-ul are urmatoarele functii pentru rezolvarea ecuatiilor diferentiale:

Functia    Tipul problemei de rezolvat Metoda de integrare

ode45 Ecuatie diferentiala cu o scara de timp Runge-Kutta

ode23 Ecuatie diferentiala cu o scara de timp Runge-Kutta

ode113 Ecuatie diferentiala cu doua scari de timp Adams

ode15s Ecuatie diferentiala cu doua scari de timp FDN

ode23s Ecuatie diferentiala cu doua scari de timp Rosenbrock

ode23t Ecuatie diferentiala cu doua scari de timp Metoda trapezelor

ode23tb Ecuatie diferentiala cu doua scari de timp TR-BDF2

unde FDN sunt formule de diferentiere numerica iar TR-BDF2 aplica mai intai o metoda Runge-Kutta iar in faza a doua o metoda de diferentiere inapoi de ordinul doi.

Sintaxa acestor functii este urmatoarea:

[T,Y] = solver(Functie,tspan,y0)

[T,Y] = solver(Functie,tspan,y0,optiuni)

unde solver este una din functiile MATLAB ode45, ode23, ode113, ode15s, ode23s, ode23t, sau ode23tb.

Argumentul Functie este o functie MATLAB ce descrie membrul drept al ecuatiei diferentiale: f(t,y). Toate aceste functii pot rezolva ecuatii de forma: y'=f(t,y)sau sisteme de ecuatii diferentiale de forma: M(t,y)y'=f(t,y), unde M(t,y) este o matrice. Functia ode23s poate rezolva numai sisteme cu matrici M constante, iar functiile ode15s si ode23t pot rezolva sisteme pentru care matricea M este singulara.

Parametrul tspan este un vector ce specifica intervalul de integrare, [t0 tf]. Pentru a obtine solutiile in anumite momente de timp se utilizeaza forma:

tspan = [t0,t1,,tf]

Parametrul y0 este un vector ce contine conditiile initiale.

Pentru utilizarea unor argumente optionale de integrare (cum ar fi precizia de integrare) se utilizeaza functia odeset.

Exemplul 1: Un exemplu de sistem cu o singura scara de timp este urmatorul, in care este descrisa miscarea unui rigid fara actiunea fortelor externe.

y₁'=y₂y₃ , y₁(0)=0

y₂'=-y₁y₃ , y₂(0)=1

y₃'=-0.51y₁y₂ , y₃(0)=1

Pentru a simula acest sistem se creeaza o functie sistem ce contine ecuatiile:

function dy = sistem(t,y)

dy = zeros(3,1);    % aloca spatiu pentru un vector coloana

dy(1) = y(2) * y(3);

dy(2) = -y(1) * y(3);

dy(3) = -0.51 * y(1) * y(2);

Precizia cu care sunt efectuate calculele poate fi modificata utilizand functia odeset. Integrarea se face pe intervalul [0 12] cu conditiile initiale [0 1 1].

optiuni = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);

[T,Y] = ode45(@sistem,[0 12],[0 1 1],optiuni);

plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')

In cazul unei ecuatii diferentiale de ordin superior, aceasta trebuie transformata intr-un sistem de ecuatii diferentiale de ordinul I de forma urmatoare:

y'=f(t,y)

Orice ecuatie de forma y⁽ⁿ⁾=f(t,y,y', . ,y^(n-1)) poate fi transformata in urmatorul sistem de ecuatii diferentiale de ordinul I, cu ajutorul urmatoarelor substitutii: y₁=y, y₂=y'_, y_n=y^(n-1). Se obtine urmatorul sistem de ecuatii diferentiale de ordinul I:

y₁'=y₂

y₂'=y₃

:

:

y_n'= f(t,y₁,y₂, . ,y_n)

Pentru rezolvarea ecuatiilor diferentiale cu conditii pe frontiera se utilizeaza functia bvp4c. Aceasta functie rezolva o ecuatie de forma y'=f(t,y)pe intervalul [a, b] care respecta conditiile pe frontiera ci(y(a),y(b))=0. Functia bvp4c are sintaxa urmatoare:

sol = bvp4c(Functie,ConditiiInitiale,solinit)

Argumentul Functie este o functie MATLAB ce descrie membrul drept al ecuatiilor diferentiale f(x,y).

Argumentul ConditiiInitiale este o functie MATLAB ce descrie conditiile pe frontiera. Ea poate avea una din formele urmatoare:

res = ConditiiInitiale(ya,yb)

res = ConditiiInitiale(ya,yb,p1,p2,)

res = ConditiiInitiale(ya,yb,parameters)

res = ConditiiInitiale(ya,yb,parameters,p1,p2,)

unde ya si yb sunt vectori coloana corespunzatori lui y(a) si y(b). Argumentul parameters este vectorul parametrilor necunoscuti iar p₁, p₂, . sunt parametrii cunoscuti. Iesirea res este un vector coloana.

Argumentul solinit este o structura cu urmatoarele campuri:

- x nodurile retelei de puncte in care se calculeaza solutiile. Conditiile la frontiera sunt impuse astfel:

a = solinit.x(1) si b = solinit.x(end).

- y aproximarile initiale ale solutiei astfel ca solinit.y(:,i) este o aproximare a solutiei in nodul solinit.x(i).

Se poate forma solinit cu ajutorul functiei bvpinit ce determina o aproximare initiala pentru functia bvp4c, asa cum este prezentat in exemplul urmator.

Exemplul 2: Sa se rezolve ecuatia y"+|y|=0, unde , cu conditiile pe frontiera:

y(0)=0

y(4)=-2

Solutie:

Mai intai ecuatia diferentiala de ordinul doi trebuie scrisa sub forma unui sistem de doua ecuatii diferentiale de ordinul I:

y₁'=y₂

y₂'=-|y₁|

unde s-au facut notatiile: y₁=y si y₂=y'.

Functia MATLAB ce descrie sistemul de ecuatii diferentiale este urmatoarea:

function dydx = sistem(x,y)

dydx = [ y(2)

-abs(y(1))];

Functia MATLAB ce descrie conditiilor initiale pe frontiera este urmatoarea:

function res = ci(ya,yb)

res = [ ya(1)

yb(1) + 2];

Apoximarea initiala a solutiei in cinci puncte uniform distribuite in intervalul [0,4] si cu aproximarile initiale y₁(x)=1 si y₂(x)=0 este data de functia:

solinit = bvpinit(linspace(0,4,5),[1 0]);

Problema este rezolvata cu comanda:

sol = bvp4c(@sistem,@ci,solinit);

Solutia numerica este evaluata in 100 de puncte uniform distribuite si este plotata cu comenzile:

y = deval(sol,linspace(0,4));

x= linspace(0,4);

plot(x,y(1,:));

O alta solutie a acestei probleme poate fi obtinuta schimband conditiile initiale:

solinit = bvpinit(linspace(0,4,5),[-1 0]);

TEME:

Sa se determine solutiile ecuatiei diferentiale urmatoare:

y"+5y'y²+y=0

pe intervalul [0, 10] in conditiile initiale: y'(0)=0, y(0)=0.5. Sa se reprezinte grafic si evolutia in planul fazelor.

Sa se determine solutiile ecuatiei diferentiale urmatoare:

pe intervalul [0, 20] in conditiile initiale: y'(0)=0, y(0)=2. Parametrul are valoarea 1.

Aceeasi problema pentru =1000.

Sa se determine solutiile ecuatiei diferentiale urmatoare:

y"'+y"+y'sin(y')+y+1=0

pe intervalul [0, 10] in conditiile initiale: y"(0)=0.2, y'(0)=0, y(0)=0.3. Sa se imparta ecranul in patru cadrane utilizand functia subplot, si sa se vizualizeze simultan dependentele (y,y'), (y,y"), (y',y") precum si o reprezentare tridimensionala a matricei y folosind functia mesh.

Aceeasi problema de la punctul 4 pentru ecuatia:

y"'+7y"+5y'+3y+1=0

in conditiile initiale y"(0)=0.4, y'(0)=0, y(0)=0.2, pe intervalul [0, 5].

Sa se determine solutiile ecuatiei diferentiale urmatoare:

y"-3y'+y=0

pe intervalul [0,10], cu conditiile pe frontiera y(0)=0; y(10)=-5.