INTEGRAREA ECUATIILOR DIFERENTIALE ORDINARE CU CONDITII INITIALE

INTEGRAREA ECUATIILOR DIFERENTIALE ORDINARE CU CONDITII INITIALE

O mare parte dintre sistemele dinamice pot fi modelate cu ajutorul ecuatiilor / sistemelor de ecuatii diferentiale, ordinare, cu conditii initiale. Prin rezolvarea (numerica) a acestora se realizeaza de fapt o simulare numerica, utila atat in etapa de analiza cat si in cea de sinteza (proiectare) a sistemelor dinamice.

Lucrarea de fata isi propune studiul, prin simulare numerica, a unui sistem dinamic de ordinul I, descris de o ecuatie diferentiala ordinara, de acelasi ordin, pe baza metodei numerice Runge-Kutta de ordinul IV. Partea a-II-a a lucrarii este dedicata studiului unui sistem dinamic descris de o ecuatie diferentiala de ordinul al-II-lea, transformata intr-un sistem de doua ecuatii diferentiale de ordinul I. Integrarea se realizeaza cu functia ode23 (MATLAB) bazata pe metoda numerica Runge-Kutta de ordin II-III, cu pas autoadaptabil.

1. SIMULAREA COMPORTARII UNUI TRADUCTOR DE TEMPERATURA

Dinamica incalzirii si racirii traductoarelor de temperatura, in ipoteza tratarii acestora drept sisteme cu parametri concentrati, este descrisa de catre ecuatia diferentiala ordinara:

unde s-au utilizat urmatoarele notatii: t este variabila independenta (timpul), T(t) este temperatura asociata semnalului generat de catre traductor, a este o constanta (de timp), iar T (t) reprezinta temperatura mediului unde este introdus traductorul. Se considera momentul initial si conditia initiala . Notand y(t):=T(t) drept iesirea sistemului dinamic si u(t):= T (t) drept intrarea sistemului dinamic, se obtine urmatoarea descriere:

Notand membrul drept al ecuatiei anterioare cu:

rezulta problema clasica a lui Cauchy, anume rezolvarea unei ecuatii diferentiale ordinare, de ordinul intai, cu conditie initiala cunoscuta.

Cunoscand constanta de timp a sistemului dinamic, a, rezulta urmatorii parametri pentru integrarea ecuatiei:

pas de integrare: ;

interval de simulare: .

APLICATIE:

Sa se realizeze, intr-un limbaj de programare (de exemplu MATLAB), un program pentru simularea comportarii unui traductor de temperatura, rezultatele reprezentandu-se grafic. Se va utiliza metoda Runge-Kutta de ordinul IV cu pas de integrare constant.

In realizarea programului se vor avea in vedere cele descrise in continuare.

Functia f se introduce printr-o functie externa utilizator, numita FCS, cu parametrii formali t si y, cu semnificatia prezentata anterior. Parametrii sistemului, anume constanta de timp a si valoarea intrarii u, se vor specifica drept variabile globale in programul principal. Functia este descrisa in limbaj pseudocod, la punctul 13 care urmeaza.

Programul, descris partial in limbajul pseudocod, va parcurge etapele prezentate in continuare.

Specificatii asupra tipului variabilelor (limbaj MATLAB):

variabile globale a si u: global a u;

Initializare preliminara, numar maxim de puncte pentru integrare:

atribuie NpMax 200

Introducere parametru model, anume constanta de timp, in variabila a. Se va testa conditia: si in caz ca aceasta nu este indeplinita, se va relua introducerea acestui parametru al modelului.

Introducere intrare in sistem (in variabila u) si conditie initiala (in variabila y0).

Introducere parametri de simulare: pas de integrare (in variabila Hint), pas de afisare sau tiparire (in variabila Htip), durata de simulare sau interval de observare (in variabila Hobs).

Calcul: numar pasi de simulare (variabila Nstep), numar pasi la care se extrag datele pentru reprezentarea grafica (variabila Ntip) si numar puncte pentru reprezentarea grafica (variabila Np). Se va verifica urmatoarea conditie: Np <= NpMax. Calculele care se efectueaza sunt urmatoarele:

atribuie Nstep fix ( Hobs/Hint )

atribuie Ntip fix ( Htip/Hint )

atribuie Np fix ( Nstep/Ntip ) + 1

* daca (Np > NpMax), atunci se reia de la punctul 5.

Initializari integrare:

atribuie t0 0

atribuie h2 Hint/2

atribuie h6 Hint/6

* initializare vectori linie: ys si timp, cu Np zerouri

* initializare vector linie: us cu Np valori egale cu u

atribuie timp(1) t0

atribuie ys(1) y0

atribuie t t0

atribuie y y0

atribuie kp 1

atribuie h Hint

scrie t, y

Calcul integrare, aplicand metoda Runge-Kutta de ordin IV:

pentru Kstep = 1,Nstep executa

atribuie t t + Hint

ê atribuie r1 FCS ( t, y )

ê atribuie r2 FCS ( t+h2, y+h2*r1 )

ê atribuie r3 FCS ( t+h2, y+h2*r2 )

ê atribuie r4 FCS ( t+h, y+h*r3 )

ê atribuie y y + h6 * ( r1 + 2*r2 + 2*r3 + r4 )

ê daca ( rem ( Kstep, Ntip ) = 0 ) atunci

ê atribuie kp kp + 1

ê ê atribuie timp(kp) t

ê ê atribuie ys(kp) y

ê ê scrie t, y

ê ð

ð

Reprezentare grafica: se reprezinta, in acelasi grafic, perechile de vectori (timp, us) [linie continua] si (timp, ys) [linie intrerupta].

Reluare program, daca se doreste, de la punctul 5, cu alti parametri de integrare (simulare).

Reluare program, daca se doreste, de la punctul 4, cu alte conditii initiale.

Reluare program, daca se doreste, de la punctul 3, cu alt parametru a al sistemului dinamic.

functia r = FCS ( t, y ) este

* variabile globale a si u: global a u;

ê atribuie aa 1/a

ê atribuie r -aa * y + aa * u

ð

Nota: pentru limbajul MATLAB, se va edita aceasta functie intr-un fisier cu numele FCS.m

Observatii:

In descrierea anterioara sunt folosite urmatoarele funcii MATLAB: global (declarare variabile globale), fix (calcul parte intreaga a unui numar), rem (restul impartirii a doua numere).

Programul se va rula pentru urmatoarele date de intrare, cu referinta la punctele 3, 4 si 5 din descrierea anterioara:

a = 1 (min.);

a = 2 (min.);

u = 90 (grd. Celsius);

y0 = -20 grd. Celsius);

u = 90 (grd. Celsius);

y0 = 0 (grd. Celsius);

u = 90 (grd. Celsius);

y0 = 150 (grd. Celsius);

Hint = 0.05 (min.);

Htip = 0.1 (min.);

Hobs = 5 (min.);

Hint = 0.1 (min.);

Htip = 0.2 (min.);

Hobs = 10 (min.);

Reprezentarile grafice obtinute se vor nota, analiza si interpreta.

2. FUNCTIA ode23 (MATLAB)

Scop Rezolvarea sistemelor de ecuatii diferentiale ordinare, de ordinul I, cu conditii initiale, prin metoda Runge-Kutta cu pas autoadaptabil, de ordinul II-III.

Apel

[tout,yout] = ode23('FCU',[t0,tfinal],y0)

sau:

[tout,yout] = ode23('FCU',[t0,tfinal],y0,tol)

sau:

[tout,yout] = ode23x('FCU',[t0,tfinal],y0,tol,trace)

Parametrii formali :

> INTRARE:

FCU numele fisierului .m ce contine functia utilizator pentru evaluarea functiilor ce constituie membrul drept al sistemului de ecuatii diferentiale ; functia FCU trebuie sa fie furnizata de catre utilizator in fisierul FCU.m si trebuie sa fie de forma :

function yprime = FCU(t,y)

parametri intrare

t variabila independenta, scalar;

y variabila dependenta, vector coloana;

parametri iesire

yprime derivata calculata, vector coloana.

FCU trebuie sa evalueze yprime(i)=dy(i)/dt, i=1,length(y) fiind date : t si y.

t0 valoarea initiala a lui t.

tfinal valoarea finala a lui t.

y0 valoarea (conditia) initiala a lui y, vector linie.

tol acuratetea dorita (valoare de test pentru adaptarea pasului), implicit = 1.e-3.

trace controleaza afisarea rezultatelor obtinute ( t, pas autoadaptat h, y ):

(fara afisare, valoare implicita);

cu afisare)

> IESIRE:

tout punctele in intervalul [t0,tfinal] in care se calculeaza solutia, vector linie. In general,aceste puncte nu sunt echidistante.

yout solutia calculata, length(tout) linii si cu atatea coloane cate componente are vectorul y.

OBSERVATIE:

Daca parametrul tol nu se specifica, atunci nu se specifica nici parametrul trace.

3. SIMULAREA COMPORTARII UNUI SISTEM MECANIC

Se considera un corp de masa m, suspendat de un resort avand coeficientul de elasticitate k, si care se misca cu frecare (coeficientul de frecare r) intr-un mediu dat, sub actiunea unei forte externe u(t). Fie y(t) deplasarea la momentul t a centrului de greutate al sistemului si viteza de deplasare a acestuia. Sistemul se considera in urmatoarele ipostaze:

(A) miscare libera, dupa ce actiunea fortei u(t) a incetat, sistemul avand o stare initiala nenula:

(B) actiunea unei forte externe constante u(t) = u, pornind din starea initiala nula

Modelul matematic al sistemului este reprezentat prin ecuatia diferentiala ordinara de ordinul al-II-lea:

Drept variabile de stare ale sistemului se considera deplasarea si viteza de deplasare . Rezulta urmatorul sistem de doua ecuatii diferentiale de ordinul intai:

unde s-au utilizat urmatoarele notatii:

Pentru alegerea conditiilor de simulare, ecuatia omogena asociata ecuatiei diferentiale se scrie sub forma:

unde se numeste factor de amortizare iar reprezinta pulsatia naturala. Rezulta relatiile de calcul:

Pe de alta parte, perioada de variatie a solutiei ecuatiei omogene este:

.

Asadar, rezulta urmatorii parametri pentru integrarea ecuatiei diferentiale:

pas de integrare: ;

interval de simulare: .

APLICATIE:

Sa se realizeze, in limbajul de programare MATLAB, un program pentru simularea evolutiei sistemului mecanic descris anterior, rezultatele reprezentandu-se grafic. Se va utiliza functia MATLAB ode23.

Programul, descris partial in limbajul pseudocod, va parcurge etapele prezentate in continuare.

Specificatii asupra tipului variabilelor (limbaj MATLAB):

* variabile globale A, B si u:    global A B u;

Initializare preliminara, ordin sistem de ecuatii diferentiale de ordinul intai:

atribuie n 2

Introducere 1 - coeficienti fizici model: coeficient de elasticitate (in variabila k) si masa corpului (in variabila m).

Introducere 2 - coeficienti fizici model: coeficient de frecare (in variabila r).

Calcul model de stare:

matrice A, vector B (8.37);

parametri model: (8.39, 8.40).

Introducere conditii initiale pentru simulare:

valoare intrare in sistem (in variabila u);

conditii initiale pentru variabilele de stare (in vectorul linie x):

pentru i = 1,n executa

scrie i

ê citeste x(1,i)

ð

Introducere parametri simulare

* scriere valoare T_n

* citire durata simulare (in variabila Hobs).

Initializari integrare:

atribuie t0 0

atribuie Tstart t0

atribuie Tstop Tstart + Hobs

scrie

Integrare (limbaj MATLAB):

[ Tout, Xout ] = ode23 ( FCU Tstart, Tstop], x, 1.e-3, 1

Afisare informatii referitoare la integrare, anume numar elemente (puncte) solutie, utilizand functia MATLAB length:

atribuie Np length (Tout)

scrie '=> numar elemente solutie = Np

Reprezentare grafica.

Se reprezinta grafic evolutiile celor doua stari ale sistemului dinamic, stocate in vectorii x1 (coloana 1 a matricei de iesire Xout), grafic cu linie continua, si respectiv x2 (coloana 2 a matricei de iesire Xout), grafic cu linie intrerupta, avand vectorul Tout furnizat pentru axa absciselor.

Reluare reprezentare grafica, daca se doreste (punctul 11).

Reluare program, daca se doreste, de la punctul 7, cu alt interval de observare Hobs.

Reluare program, daca se doreste, de la punctul 6, cu alte conditii initiale.

Reluare program, daca se doreste, de la punctul 4, cu alt coeficient de frecare r.

Reluare program, daca se doreste, de la punctul 3, cu alti coeficienti fizici ai modelului.

functia xp = FCU ( t, x ) este

* variabile globale A, B si u: global A B u;

ê atribuie xp A * x + B * u

ð

Nota pentru limbajul MATLAB, se va edita aceasta functie intr-un fisier cu numele FCU.m

OBSERVATIE:

Programul se ruleaza pentru urmatoarele date de intrare:

3. k = 16. (N/m)

m . (kg)

4.1 r = 7. (N/m/s)

4.2 r = 3. (N/m/s)

4.3 r = 0. (N/m/s)

6.1 u = 0. (N)

x1(0) = 2. (m)

x2(0) = -0.64 (m/s)

6.2 u = 5. (N)

x1(0) = 0. (m)

x2(0) = 0. (m/s)

7. Hobs = 10. (s)