Calculul Functiilor de Matrice

Politehnica

Calculul Functiilor de Matrice

Byadmin ianuarie 1, 2024

Cuprins

Capitolul 1. INTRODUCERE

Capitolul 2. APLICAȚII DE INTRARE – IEȘIRE PENTRU SISTEME DINAMICE

2.1. Aplicații de intrare-iesire pentru sisteme dinamice liniare continue

2.2. Sisteme izomorfe.

2.3. Sisteme discrete

2.4. Aplicații în Matlab

Capitolul 3. SISTEME LINIARE. CALCULUL FUNCȚIILOR DE MATRICE

3.1. Funcții de matrice

3.2. Calculul exponențialei matriceale

3.3. Proprietãțile matricei exponențiale

3.4. Calculul matricei exponențiale folosind transformãrile Laplace

3.5. Importanța polinomului caracteristic

3.6. Teoria Dunford-Taylor pentru matricea f(A)

Capitolul 4. METODE DE CALCUL PENTRU DETERMINAREA COMPORTÃRII DINAMICE A SISTEMELOR LINIARE

4.1. Metode de calcul a matricei de tranziție

Concluzii

Bibliografie

Capitolul 1

INTRODUCERE

Capitolul 2

APLICAȚII DE INTRARE – IEȘIRE PENTRU SISTEME DINAMICE

2.1 Aplicații de intrare-iesire pentru sisteme dinamice liniare continue

Fie sistemul liniar descris de

ẋ(t)=A(t)x(t)+B(t)u(t) (2.1)

y(t)=C(t)x(t)+D(t)u(t) (2.2)

unde x(t0)=x€n, t0=momentul inițial.

Fie ɸ(t,t0) matricea fundamentalã pentru A(t). ɸ(t,t0) este matricea care are coloanele x1(t), x2(t), … ,xn(t), unde xi(t) este soluția problemei Cauchy ẋi(t)=A(t)xi(t), xi(t0)=ei=(0…0 0…0)’.

Matricea fundamentalã are urmãtoarele proprietãți:

i) ɸ(t,t0)=A(t)ɸ(t,t0);

ii) ɸ(t,t0)=I;

iii) x(t)=ɸ(t,t0)x0 este soluția problemei Cauchy

ẋ(t)=A(t)x(t), x(t0)=x0;

iv) ɸ(t,t0) ɸ(t1,t0)= ɸ(t,t0);

v) ɸ(t,t0)-1= ɸ(t0,t).

Soluția sistemului de ecuații diferențiale (2.1) este datã de formula variației constantelor:

x(t)=ɸ(t,t0)x0+ɸ(t,s)b(s)ds.

Folosind formula aceasta pentru f(t)=B(t)u(t) obținem formule pentru sistemul intrãrii u și stãrii inițiale x(t0)=x€n:

x(t)= ɸ(t,t0)x0+ (2.3)

Utilizând formula y(t)=C(t)x(t)+D(t)u(t), se obține aplicația intrare-ieșire a sistemului ca fiind datã de

y(t)=C(t) ɸ (t,t0)x0+ ɸ(t,s)B(s)u(s)ds+D(t)u(t) (2.4)

Pentru sistemele staționare, matricea fundamentalã este obținutã cu ajutorul matricei exponențiale a matricii A:

ɸ(t,t0)=eA(t-t0) ∆=I+

și in acest caz relația y(t)=C(t) ɸ (t,t0)x0+ɸ(t,s)B(s)u(s)ds+D(t)u(t) devine:

y(t)= CeA(t-t0)x0+. (2.5)

Dacã u(t)≡0 ieșirea este

y(t) |liber=C(t) ɸ (t,t0)x0

și se numește rãspunsul la intrare nulã al sistemului.

Dacã x0=0 se obține rãspunsul la starea nulã:

y(t)|perturbat= ɸ(t,s)B(s)u(s)ds+D(t)u(t).

Dacã se aplicã un impuls de tip Dirac ui(t)=δ(t-t0)ei, întrucât () pentru orice ieșirea corespunzãtoare este coloana i a matricei:

R(t,t0)=

matrice care se numește rãspunsul la impuls al sistemului. Cu ajutorul acestei matrice rãspunsul la starea nulã se scrie:

y(t)=

Pentru un sistem staționar rãspunsul la impuls este:

R(t)= ,

iar aplicația intrare-ieșire are expresia

y(t)=

Exemplul 2.1

Se considerã circuitul electric:

Fig.2.1

Intensitatea curentului i(t) și tensiunea v(t) verificã ecuația diferențialã:

L

Considerând drept stare, dar și ieșire x(t)=i(t)=y(t) și drept intrare u(t)=v(t) obținem reprezentarea în spațiul stãrilor a acestui circuit:

Aplicația de intrare-ieșire pentru sistemul cu o singurã intrare și o singura ieșire (SISO) este descris în acest exemplu, conform formulei (2.5)

Exemplul 2.2

Aplicația de intrare-ieșire a sistemului

este, pentru t0=0

și continuând calculele rezultatul va fi urmãtorul

y(t)=e-2tx10+2e-3tx20+2.2 Sisteme izomorfe

În continuare voi arãta efectul pe care îl are asupra unui sistem staționar o schimbare de bazã în spațiul stãrilor.

Se va considera o transformare de forma x=T, unde T este o matrice n x n nesingularã (inversabilã). Prin înlocuirea lui x în (2.1) și (2.2) obținem un sistem izomorf, și anume:

=T-1 AT +T-1Bu

ŷ=CT +Du

Sistemele =(A,B,C,D) și=() sunt izomorfe dacã și numai dacã existã o matrice nesingularã T astfel încât sã fie îndeplinite relațiile:

Â=T-1AT, =T-1B, =CT, =D. (2.6)

Deoarece eÂt=T-1eAtT, aplicația intrare-ieșire a sistemului (calculatã conform (2.5)) este:

=

=.

Aplicațiile intrare-ieșire pentru douã sisteme izomorfe sunt egale (același lucru este valabil și pentru rãspunsurile la impuls). Relația intrare-ieșire este independentã de alegerea stãrilor (interne).

Se va considera sistemul staționar (2.1), (2.2) și schimbarea de coordonate x(t)=T(t), (unde T: ia ca valori matrice nesingulare, este diferențiabilã și în plus, T, T-1, sunt mãrginite pe fiecare interval mãrginit de numere reale (T se mai numește transformare Liapunov).

Înlocuind și efectuând calculele se va obține:

ẋ(t)=Ṫ(t)(t)+T(t)(t)=A(t)T(t)(t)+B(t)u(t),

de unde, mai departe
(t)=[T(t)-1A(t)T(t)-T(t)-1Ṫ(t)](t)+T(t)-1B(t)u(t)

(t)=C(t)T(t)(t)+D(t).

Rezultã cã se poate considera sistemele nestaționare și =() ca fiind izomorfe dacã existã o transformare Liapunov T(t) astfel încât

(2.7) Â(t)=T(t)(t)T(t)-T(t)-1(t), (t)=T(t)-1B(t),

Ĉ(t)=C(t)T(t), (t)=D(t).

Dacã alegem pe post de transformare Liapunov matricea fundamentalã T(t)=ɸ(t,0), atunci Ṫ(t)=A(t)ɸ(t,0), T-1(t)=ɸ(0,t) avem:

Â(t)=T(t)-1A(t)T(t)-T(t)-1Ṫ(t)=ɸ(0,t)A(t)ɸ(t,0)-ɸ(0,t)A(t)ɸ(t,0)=0,

(t)=T(t)-1B(t)=ɸ(0,t)B(t),

(t)=C(t)T(t)=C(t)ɸ(t,0),

(t)=D(t),

adicã sistemul izomorf obținut aratã astfel

Considerarea unui sistem izomorf poate duce de multe ori la simplificarea rețelei interne a unui sistem dat (eliminarea unui numãr important de conexiuni între componentele primare, precum și între aceste intrãri/ieșiri, cu pãstrarea aplicației intrare-ieșire) dupã cum urmeazã în urmãtoarele exemple.

Exemplul 2.3

Se va simplifica urmãtoarea diagramã:

Fig. 2.2

Matricele sistemului sunt:

A=, B=, C=[-1 3], D=[5].

-1 -5

-5 -1

calculând doi vectori proprii corespunzãtori se va obține v1=, v2= (noua bazã a spațiului X=2).

Matricea schimbãrii de bazã este T=[v1 v2] = și calculând matricea inversã se obține

T-1=

Un sistem izomorf este , unde

-1AT=

T-1B=, [2 -4],

Diagrama sistemului este data de urmãtoarea figurã 2.3 și se va observa reducerea numãrului de conexiuni:

Fig. 2.3

2.3 Sisteme discrete

Se considerã un sistem nestaționar

(unde K este un corp), cu reprezentarea în spațiul stãrilor

x(t+1)=A(t)x(t)+B(t)u(t) (2.8)

y(t+1)=C(t)x(t)+D(t)u(t) t=0,1,2,… (2.9)

dacã B(t)=0, avem de-a face cu un sistem omogen

x(t+1)=A(t)x(t). (2.10)

Obținem

x(1)=A(0)x(0), x(2)=A(1)x(1)=A(1)A(0)x(0),…,;

X(t)=A(t-1)A(t-2)…A(1)A(0)x(0) (2.11)

Matricea fundamentalã a sistemului discret (2.10) este matricea

ɸ(t,s)= (2.12)

unde,

A(t-1)A(t-2)…A(s+1)A(s) ,t>s

I ,t=s

Se poate observa ca ɸ(t,t)=I și de asemenea cã

ɸ(t+1,s)= =A(t)

Atunci (2.11) devine x(t)=ɸ(t,0)x(0).

Principalele proprietãți ale matricei fundamentale discrete sunt:

x(t)=ɸ(t,s)x(s) ,

ɸ(t,r)ɸ(r,s)=ɸ(t,s) , (2.13)

Dacã matricele A(k) sunt inversabile pentru orice , atunci ɸ(t,s) poate fi definitã și pentru t<s prin urmãtoarea formulã:

ɸ(t,s)=[A(s-1)A(s-2)…A(t+1)A(t)]-1

și în acest caz, ɸ(t,r)ɸ(r,s)=ɸ(t,s) poate sã rãmânã adevãratã pentru orice și avem și

[ɸ(t,s)]-1=ɸ(s,t).

Putem sã calculãm starea și ieșirea sistemului . Din (2.8) rezultã

x(1)=A(0)x(0)+B(0)u(0)=ɸ(1,0)x(0)+ɸ(0,0)B(0)u(0);

x(2)= A(1)A(0)x(0)+A(1)B(0)u(0)+B(1)u(1)=ɸ(2,0)x(0)+ɸ(2,1)B(0)u(0)+ɸ(2,2)B(1)u(1);

x(3)=A(2)A(1)A(0)x(0)+A(2)A(1)B(0)u(0)+A(2)B(1)u(1)=

=ɸ(3,0)x(0)+ɸ(3,1)B(0)u(0)+ɸ(3,2)B(1)u(1)+ɸ(3,3)B(1)u(2);

………………………

x(t)=A(t)

Deci conform formulei (2.8)

x(t+1)=A(t)

=

Formula pentru x(t) poate fi rescrisã astfel:

x(t)=ɸ(t,0)x(0)+ɸ(t,k+1)B(k)u(k),

deci aplicația de intrare-ieșire este datã de

y(t)=ɸ(t,0)x(0)+ɸ(t,k+1)B(k)u(k)+D(t)u(t). (2.15)

Ca și în cazul sistemelor continue, aceastã formulã generalã de rãspuns are douã componente: componenta liberã (rãspunsul la intrare nulã)

liber=C(t)ɸ(t,0)x(0)

și componenta forțatã (rãspunsul perturbat)

forțat=

unde matricea de rãspuns la impuls R este definitã prin

R(t,k)= (2.16)

Dacã sistemul este staționar, atunci ɸ(t,s)=At-s pentru orice (existã și pentru t<s doar dacã A este ireversibilã), iar (2.14) , (2.15), (2.16) devin

x(t)=Atx(0) (2.17)

y(t)=CAtx(0) (2.18)

(2.19)

Exemplul 2.4

Se va da sistemul liniar

Fig.2.4

și are reprezentarea de stare:

(t+1)=

y(t)=

iar șirul stãrilor este

=

iar aplicația intrare-ieșire este

y(t)=.

Exemplul 2.5

Fie sistemul unde

, B=, C=[1 -1], D=7.

Polinomul caracteristic al matricii A este

, deci valorile proprii ale lui A sunt 1=2, 2=-1.

Folosind teorema Hamilton-Cayley (A2=A-2I, A3=A2-2A=-A-2I,…) obținem matricea fundamentalã:

ɸ(t,0)=At=

deci valorile proprii verificã relațiile

adicã 2t=.

Se obține

deci

ɸ(t,0)=At=T

Starea sistemului și ieșirea acestuia sunt:

Se poate sã se calculeze aplicația intrare-ieșire folosind sistemul izomorf dat de matricea T=[v1 v2] (v1, v2 vectori proprii).

De multe ori unele sisteme discrete se obțin din modele continue, de exemplu atunci când este folositã procesarea digitalã a semnalelor. În acest caz, atât pentru intrãri cât și pentru ieșiri se considerã „eșantioane” obținute cu un pas fix h (de exemplu secunde).

Fig. 2.5

Se va considera o intrare constantã pe porțiuni u(t)=un, , unde este un șir dat (Fig.2.5).

Notãm și

, (2.20)

Starea sistemului staționar

pentru t=(n+1)h și t0=nh este

,

deci sistemul discretizat este

(2.21)

Exemplul 2.6

Fie sistemul dinamic continuu

y(t)=[1 1]

Atunci

Sistemul discret(izat) corespunzãtor are ecuația de stare

în timp ce ecuația de ieșire rãmâne neschimbatã.

2.4 Aplicații în Matlab

lsim

Syntax: lsim(sys,u,t); lsim(A,B,C,D,u,t)

Aceste comenzi simuleazã și acționeazã forțele de rãspuns la LTI (continuu și discret SISO sau MIMO) modelul sys si intrarea arbitrarã u. Vectorul t specificã exemple de timp pentru simulare t = 0:dt:Tfinal.

Matricea u trebuie sã aibã l afel de mulet coloane sau rânduri precum monstre de timp (lungimea t) și l afel de mulet coloane/ rânduri precum intrãrile în sistem. Fiecare rând coloanã u(i,:)/u(:,i) specificã valoarea intrãrii s pentru monstra de timp t(i). Daca sys este în modul discret atunci u trebuie sã fie indicat la aceeși ratã (valoare) ca sistemul (t este atunci reduntant si poat efi omis sau pus pe o matrice goalã). În cazul timpului continuu monstra de timp dt=t(i+1)-t(i) ese folositã ca sã discretizeze (sã compromitã modelul continuu).

lsim(sys,u,t,x0) sau lsim (A,B,C,D,u,t,x0) comutatorul de rãspuns general general produc eîn stare inițiala x0 sau o intrare u.

Example 3.4.1.

A=[-2 0;0 -3]; B=[1 0;0 1]; C=[1 2]; D=[1 1]; sys=ss(A,B,C,D); x0=[1 0];t=0:0.2:2;

u=ones(length(t),2);

lsim(sys,u,t,x0) (or lsim(sys,u,t))

Funcția lsim poate fi folositã cu argumentele din stânga.

[y,t] =lsim(sys,u,t)

[y,t,x] = lsim(sys,u,t) \% doar pentru modelele tipice poziție spațiu.

[y,t,x] = lsim(sys,u,t,x0) \% cu poziție inițialã x0.

Fig.3.4.1

În acest caz comutarea nu este automat indicatã (desenatã) pe ecarn , dar valorile ieșirii y de rãspuns, vectorul timp t, folosit pentru simulare și traiectoriile de stare (poziție) sunt furnizate. Matricea x și y au lungimea t pe rânduri și respectiv n și p pe coloane.

Fig.2.1 (a)

A=[0 1;-2 -3]; B=[0;1]; C=[1 1]; D=0;

x0=[1 0.5]; t=0:.05:4;u=t; %generates a vector u

[y,x]=lsim(A,B,C,D,u,t,x0);

plot(t,x,t,y)

Fig.2.1(a)

Fig.2.1 (b)

A=[0 1;-2 -3]; B=[0;1]; C=[1 1]; D=0;

x0=[1 0.5]; t=0:.05:4;u=t; %generates a vector u

[y,x]=lsim(A,B,C,D,u,t,x0);

plot(t,x,t,y)

text(0.1,1.45,'y'),text(0.05,0.9,'x1'),text(0.05,0.5,'x2')

Fig.2.1(b)

Fig. 2.2 (a)

A=[-2 0;0 -3];B=[1 0;0 1];C=[1 2];D=[1 -1];

x0=[1 0];t=0:0.2:2;u=[t;t];

[y,x]=lsim(A,B,C,D,u,t,x0);

plot(t,x,t,y) % (utomat indicatã (desenatã) pe ecarn , dar valorile ieșirii y de rãspuns, vectorul timp t, folosit pentru simulare și traiectoriile de stare (poziție) sunt furnizate. Matricea x și y au lungimea t pe rânduri și respectiv n și p pe coloane.

Fig.2.1 (a)

A=[0 1;-2 -3]; B=[0;1]; C=[1 1]; D=0;

x0=[1 0.5]; t=0:.05:4;u=t; %generates a vector u

[y,x]=lsim(A,B,C,D,u,t,x0);

plot(t,x,t,y)

Fig.2.1(a)

Fig.2.1 (b)

A=[0 1;-2 -3]; B=[0;1]; C=[1 1]; D=0;

x0=[1 0.5]; t=0:.05:4;u=t; %generates a vector u

[y,x]=lsim(A,B,C,D,u,t,x0);

plot(t,x,t,y)

text(0.1,1.45,'y'),text(0.05,0.9,'x1'),text(0.05,0.5,'x2')

Fig.2.1(b)

Fig. 2.2 (a)

A=[-2 0;0 -3];B=[1 0;0 1];C=[1 2];D=[1 -1];

x0=[1 0];t=0:0.2:2;u=[t;t];

[y,x]=lsim(A,B,C,D,u,t,x0);

plot(t,x,t,y) % (variant plot(t,u,t,x,t,y) text(0.05,1.05,'u'),)

Fig.2.2 (a)

Fig. 2.2 (b)

A=[-2 0;0 -3];B=[1 0;0 1];C=[1 2];D=[1 -1];

x0=[1 0]; t=0:0.2:2;u=[t;t];

[y,x]=lsim(A,B,C,D,u,t,x0);

plot(t,x,t,y) % (variant plot(t,u,t,x,t,y) text(0.05,1.05,'u'),)

text(0.05,1.05,'y'),text(0.02,0.8,'x1'),text(0.05,0.05,'x2')

Fig.2.2 (b)

Fig. 2.2.2(a)

A=[-2 0;0 -3];B=[1 0;0 1];C=[1 2];D=[1 -1]; x0=[1 0]; t=0:0.2:2; u=[t;sign(t-1)];

[y,x]=lsim(A,B,C,D,u,t,x0); plot(t,x,t,y)

Fig. 2.2.2(b)

A=[-2 0;0 -3]; B=[1 0;0 1]; C=[1 2]; D=[1 -1]; x0=[1 0]; t=0:0.2:2; u=[t;sign(t-1)];

[y,x]=lsim(A,B,C,D,u,t,x0); plot(t,x,t,y)

text(0.05,2.05,'y'),text(0.02,0.8,'x1'),text(0.05,0.05,'x2')

Fig. 2.2.3(a)

A=[-2 0;0 -3];B=[1 0;0 1];C=[1 2];D=[1 -1];

x0=[1 0];t=0:0.2:2;u=[t; ones(1, length(t))];

[y,x]=lsim(A,B,C,D,u,t,x0);

plot(t,x,t,y)

Fig. 2.2.3(b)

A=[-2 0;0 -3];B=[1 0;0 1];C=[1 2];D=[1 -1];

x0=[1 0];t=0:0.2:2;u=[t; ones(1, length(t))];

[y,x]=lsim(A,B,C,D,u,t,x0);

plot(t,x,t,y)

text(0.05,0.15,'y'),text(0.02,0.8,'x1'),text(0.1,0.05,'x2')

Fig.2.2.3(c)

A=[0 1;-2 -3];

B=[0;1];

C=[1 1];

D=0;

x0=[1 0.5];

t=0:.05:4;

u=0*t;

sys=ss(A,B,C,D);

lsim(sys,u,t,x0) (or lsim(sys,u,t))

Fig.2.2.3(d)

A=[0 1;-2 -3];

B=[0;1];

C=[1 1];

D=0;

t=0:.05:4;

u=(1/2)*ones(1, length(t));

sys=ss(A,B,C,D);

lsim(sys,u,t)

Fig.2.2.3(e)

A=[0 1;-2 -3];

B=[0;1];

C=[1 1];

D=0;

t=0:.05:4;

u=sign(t-1);

sys=ss(A,B,C,D);

lsim(sys,u,t)

Fig.2.2.3(f)

A=[0 1;-2 -3];

B=[0;1];

C=[1 1];

D=0;

x0=[0 0];

t=0:.05:4;

u=(1/2)*sign(t-2);

sys=ss(A,B,C,D);

lsim(sys,u,t,x0)

Fig.2.2.3(g)

A=[0 1;-2 -3];

B=[0;1];

C=[1 1];

D=0;

x0=[0 0];

t=0:.05:4;

u=t/3;

sys=ss(A,B,C,D);

lsim(sys,u,t,x0)

Fig.2.3.3(h)

A=[0 1;-2 -3];

B=[0;1];

C=[1 1];

D=0;

x0=[1 -1];

t=0:.05:4;

u=t/3;

sys=ss(A,B,C,D);

lsim(sys,u,t,x0)

Fig.2.3.3(i)

A=[0 1;-2 -3];

B=[0;1];

C=[1 1];

D=0;

x0=[0 0];

t=0:.05:4;

u=t/3;

sys=ss(A,B,C,D);

lsim(sys,u,t,x0)

Fig.2.3.3(j)

A=[0 1;-2 -3];

B=[0;1];

C=[1 1];

D=0;

x0=[2 -1];

t=0:.05:4;

u=t/4;

sys=ss(A,B,C,D);

lsim(sys,u,t,x0)

Fig. 3.4.3

A=[2 -1 3;1 0 1;-1 2 0];

B=[0 0;1 -2;0 1];

C=[1 -1 1;1 0 2];

D=0;

Ts=0.05;

sys=ss(A,B,C,D,Ts);

x0=[1 -1 0]; t=0:.05:4;

u=[t/2;sign(t-1)];

lsim(sys,u,t,x0)

Fig. 3.4.4

A=[0 1;-2 -2];

B=[0 1;1 -2];

C=[1 -1;0 2];

D=0;

sys=ss(A,B,C,D);

t=0:.05:4;

step(sys,t) % or step(sys)

Sisteme izomorfe

A=[-7 -4 4;10 6 -8;4 3 -5];

B=[1 -1;0 2;3 -1];

C=[-8 -5 5;3 3 -4];

D=[3 -4;1 0];

[T,A1]=eig(A,'nobalance')

T =[-0.5000 -0.5000 -0.0000;1.0000 1.0000 1.0000;0.5000 0.2500 1.0000];

A1 =[-3.0000 0 0; 0 -1.0000 0; 0 0 -2.0000]

sys1=ss(A1,inv(T)*B,C*T,D)

Capitolul 3. SISTEME LINIARE. CALCULUL FUNCȚIILOR DE MATRICE.

Calculul funcțiilor de matrice este necesar în multe aplicații. Analiza comportãrii dinamice a sistemelor liniare continue face apel la exponențiala matricealã, pentru calculul cãreia a fost elaboratã o întreagã serie de tehnici numerice. În acest capitol se vor prezenta câteva metode generale de calcul pentru funcțiile de matrice, precum și metodele specifice, cele mai apreciate, de calcul pentru funcția exponențialã de argument matriceal.

3.1 Funcții de matrice

Considerãm o funcție definitã pe o mulțime D din planul complex și fie o matrice datã. Vom defini noțiunea de matrice adicã semnificația expresiei

F=f(A) (3.1)

Dacã f este o funcție polinomialã

(3.2)

atunci matricea

(3.3)

este bine definitã și poate fi numitã valoarea polinomului f în punctul (sau pe matricea) A.

Fie acum (z) polinomul minimal al matricei A și i=1:l, zerourile acestuia, zeroul având ordinul de multiplicitate . Avem

(z)=zm+ (3.4)

unde

(3.5)

Conform teoremei împãrțirii cu rest, existã unic douã polinoame q, r cu grad(r)<m astfel încât

f= (3.6)

Polinomul minimal fiind un polinom anulator pentru matricea A, i.e (A)=0, din (3.6) rezultã

(3.7)

Valoarea polinomuluik f în punctul A este aceeași cu cea a polinomului r și același lucru se poate spune și despre orice alt polinom al cãrui rest la împãrțirea prin este r. Polinomul r poate fi determinat prin aplicarea algoritmului de împãrțire cu rest a polinoamelor.

O altã modalitate de calcul a polinomului r se bazeazã pe faptul cã din (3.4) rezultã

, k=0 : mi-1, (3.8)

și în consecințã, cei m coeficienți ai polinomului r se pot determina, conform (3.6), prin rezolvarea sistemului de m ecuații liniare

, i=1 : l, k=0 : mi-1 (3.9)

unde indicele superior indicã ordinul derivatei.

Relațiile (3.7) și (3.9) servesc ca bazã pentru definirea valorii în punctul A a oricãrei funcții f care admite derivatele cerute în (3.9).

Exemplul 3.1

Fie f(z)=ez și

Avem (z)=z2-3z+2 și {(1, 1), (2, 1)}. Dacã r(z)=r1z+r0, sistemul (3.9) se poate scrie asfel

,

de aici rezultã

r0=2e-e2

r1=e2-e

și

F=eA=(e2-e)A+(2e-e2)I3=

Teorema 3.1 Fie spectrul matricei și un domeniu cu frontiera suficient de netedã asfel încât . Dacã f este o funcție analiticã pe atunci

f(A)= (3.10)

Expresia (3.10) poate servi ca definiție pentru funcțiile analitice, iar calculul integralei Cauchy

fij= (3.11)

poate fi efectuat cu ajutorul teoremei reziduurilor de mai jos.

Teorema 3.2 (teorema rezidurilor). Fie f olomorfã în domeniul D exceptând singularitãțile izolate z1, z2,…, zn și fie C un Contur simplu de D omotop cu o constantã (C~0) în D care conține la interior punctele z1,…, zn.

Atunci

(3.12)

Demonstrația de bazeazã pe egalitatea

(3.13)

Exemplul 3.2

Considerãm funcția f și matricea A din exemplul 3.1

Avem {1, 2, 2} și

(zI-A)-1=.

Alegem un domeniu simplu conex oarecare ce satisface Se obține

celelalte elemente ale matricei F=eA, sunt nule.

În continuare se va prezenta cazul funcțiilor analitice , unde se va utilize relația (3.10) ca o relație definitorie a funcțiilor de matrice.

Posibilitatea exprimãrii funcțiilor de matrice prin serii matriceale de puteri dupã mãtoarea teoremã:

Teorema 3.2 Dacã funcția f(z) se poate dezvolta în serie de puteri în jurul punctului z=z0

f(z)= (3.14)

și seria este convergentã în discul |z-z0|<r atunci aceastã dezvoltare rãmâne valabilã dacã argumentul scalar este înlocuit cu argumentul matriceal A

f(A)= (3.15)

oricare ar fi matricea A al cãrei spectru se aflã în interiorul discului de convergențã.

Din aceastã teoremã rezultã, cã dezvoltãrile

eA=

sinA=

cosA=

sunt valabile oricare ar fi matricea

Exemplul 3.3

Considerãm funcția f și matricea A din exemplul 3.1. Prin ucție rezultã

deci,

.

3.2 Calculul exponențialei matriceale

În acest paragram prezentãm principalele metode pentru calculul exponențialei matriceale

ɸ(t)=etA (3.16)

unde t>0 este un parametru scalar. Funcția ɸ(t) este matrice de tranziție a stãrilor sistemului liniar , adicã satisface ecuația diferențialã matricealã cu condiția inițialã ɸ(0)=I.

Având în vedere importanța unui calcul fiabil al exponențialei matriceale pentru studiul sistemelor liniare, o discuție prealabilã a condiționãrii numerice a problemei calculuilui lui ɸ(t)=etA, adicã a sensibilitãții lui ɸ(t) în raport cu variațiile lui A, este necesarã, întrucât erorile de calcul ale lui ɸ(t) se traduc, prin mecanismul de propagare inversã, în perturbații la nivelul datelor d eintrare, adicã al matricei A.

Notând cu exponențiala asociatã matricei perturbate A+E avem =(A+E) cu (0)=I. Punând =, se obține

unde deci Pentru simplitate, neglijãm termenul și, prin integrare, deducem,

Utilizând norma matricealã putem scrie

relație din care rezultã

(3.17)

unde,

(3.18)

este numãrul de condiție al lui A relativ la calculul lui

Obținem o evaluare grosierã a lui observând cã

unde

De asemenea, avem

deci

(3.19)

În concluzie, pentru a limita sensibilitatea lui în raport cu variațiile lui A e necesar sã limitãm , de exemplu astfel încât

, (3.20)

fie alegând t suficient de mic, fie (în cazul în care t este impus) utilizând proprietatea

(3.21)

Într-adevãr, oricare ar fi existã m astfel încât , iar dacã e cunoscut atunci, conform (3.21), etA poate fi calculat printr-un proces de ridicare succesivã la pãtrat. O evaluare a lui m rezultã observând cã dacã atunci trebuie sã avem log2(t, deci

M=1+[log2(t (3.22)

unde este partea întreagã a argumentului.

Aproximația polinomialã (Taylor) de ordin p este de forma

(3.23)

unde coeficienții ck se determinã din condițiile

k=0 : p. (3.24)

Obținem

ck= (3.25)

adicã cea mai bunã aproximare polinomialã localã de ordin p coincide cu polinomul Taylor corespunzãtor.

Exemplul 3.4

Pentru p=1 din (3.23) obținem T1=1+z iar aproximația T1(tA)=I+tA a lui etA corespunde integrãrii numerice a ecuației diferențiale prin metoda Eucler explicitã de ordin 1. În mod similar se poate arãta ca metoda Runge-Kutta de ordin 4 corespunde alegerii p=4 în (3.23)

Putem evalua eraorea de trunchiere comisã prin utilizarea aproximației polinomiale (3.23) observând cã

unde paranteza dreaptã e majoratã prin

sau (dacã |z|<p+2) prin

Punând z=tA și utilizând prima evaluare, obținem eroarea relativã de trunchiere sub forma

(3.26)

Asfel spus, avem

(3.27)

unde tol este precizia de calcul impusã, dacã alegem p astfel încât

(3.28)

Pentru oricare, menționãm cã dacã și tol este de ordinul 10-5, atunci condiția (3.28) impune

Evaluare apolinomului matreiceal Tp(A) are loc dupã schema

Tk+1(tA)=Tk(tA)+Xk, Xk= (3.29)

Unde între doi termeni succesivi Xk și Xk-1 are loc relația

(3.30)

iar inițializãrile sunt, T0(tA)=I și X0=I.

II. Aproximația raționalã (Padé) de grade (p, q) și de ordin p+q este de forma

, (3.31)

unde coeficienții ck, dk se determinã din condițiile

k=0 : p+q (3.32)

Echivalent, se poate utilize metoda coeficienților nedeterminți, scriind formal

și egalând coeficienții gradelor egale din cei membri pânã la gradul p+q necesar. Rezultã

(3.33)

Exemplul 3.5

Pentru q=0 obține, aproximația polinomialã (). Printre aproximații raționale “diagonale”, cu p=q, sunt

cu polul z=2, respectiv

cu polii z1,2 =3 având Re z1,2 =3, plasați pe cercul de razã

Prin urmare, trecând la cazul matriceal, se constatã cã în general existența aproximației raționale Rpq(tA) nu este asiguratã necondiționat, întrucât matricea Dq(tA) poate fi singularã; de exemplu,

D1(tA) =

Este singularã dacã t = 1 și A are o valoare proprie În general, se poate arãta cã toate zerourile lui Dq(z) sunt situate în Re z > 0 și se acumuleazã la pentru q astfel încât Dq(tA) rezultã sigur nesingularã dacã este satisfãcutã cel puțin una din urmãtoarele condiții

i) Re, i.e. matricea A este stabile;

ii) q este suficient de mare;

iii)

În consecințã, condiția de limitare a normei lui tA este utilã pentru a asigura nu numai buna condiționare a problemei calculului lui etA ci și existența aproximației raționale considerate mai sus, independent de caracteristicile de stabilitate ale matricei date A sau de ordinal q adoptat în calcule.

Pentru a conecta exemplele (3.4) și (3.5) se poate observa cã aproximația raționalã R11(tA) corespunde integrãrii numerice a ecuației diferențiale prin metoda Euler implicate de ordin 1.

Eroarea de trunchiere poate fi limitatã alergând convenabil gradele p, q ale aproximãrii (3.30). Se poate arãta ca dacã atunci

Rpq(tA)=et(A+E) (3.34)

unde AE=EA iar

(3.35)

Deoarece AE=EA implicã et(A+E)=, eroarea relative de trunchiere comisã prin utilizarea aproximației Padé Rpq(tA) poat efi evaluate prin

. (3.36)

și aceastã inegalitate poate fi utilizatã în locul lui (3.26). Evaluarea (3.35) valabilã în cazul general , este simetricã în raport cu p și q deci nu permite alegerea lor univocal.

3.3 Proprietãțile matricei exponențiale

Proprietãțile urmãtoare, sunt proprietãți ale matricei de tranziție ce variazã în sisteme de timp variat:

Proprietatea 1. Funcția are ca soluție unicã

Proprietatea 2. Coloana i a funcției are ca soluție unicã

unde ei al vectorului i are baza canonicã din

Proprietatea 3. Pentru orice t,

Proprietatea 4. Pentru orice este nonsingular și

Pentru sistemele LTI, matricea de tranziție are proprietăți importante, care derivă din teorema Cayley-Hamilton. Se va da în continuare un polinom

și matricea A definitã pe

care este de asemenea o matrice.

Teorema Cayley-Hamilton. Pentru matricea A,

unde,

este polinomul caracteristic al matricei A.

Urmãtoarele proprietãți ale matricei exponențiale sunt consecințe ale teoremei Cayley-Hamilton și sunt specifice cazului de timp invariat.

Proprietatea 5. Pentru fiecare matrice A, existã n al funcțiilor scalare pentru relația

(3.37)

unde ai sunt coeficienții polinomului caracteristic al matricei A. Atunci

Folosind aceste relații, se poate ajunge la

Prin urmare An+1 poate fi scrisã ca o combinație liniarã de forma An-1, An-2,…, A, I. Aplicând proceduri identice pentru creșterea puterii matricei A, pentru fiecare k 0, unde Ak poate fi scris asfel

(3.38)

pentru coeficienții adecvați ai(k). Înlocuind aceastã relație eAt, avem ca rezultat

Relația (3.37) este definitã de

Proprietatea 6. Pentru fiecare a matricii A

Aceasta este o consecințã a proprietãții 5.

3.4 Calculul matricei exponențiale folosind transformãrile Laplace

Am observat în proprietatea (P.1) cã este definitã unic asftel

Luând transformata Laplace de fiecare parte a ecuației diferențiale, ajungem la concluzia urmãtoare

L L

Prin urmare, putem folosi inversa Laplace, transformãrii tabelare pentru a calcula

L [(sI-A)-1].

Notã: Din moment ce folosim transformãrile Laplace unilateral, aceastã metodã dã valori pentru

3.5 Importanța polinomului caracteristic

Folosind relația

unde

este polinomul caracteristic al matricei A, ale cãror rãdãcini are valori proprii ale lui A, iar este adjunctul matricei sI-A, ale cãror intrãri sunt polinoame în s la puterea n-1 sau mai mici.

Pentru a calcula inversa transformãrii Laplace , avem nevoie sã efectuãm fracția parțialã a fiecãrei intrãri a lui .

Acestea sunt de forma

Inversa transformãrii Laplace va fi

L-1

Astfel, atunci când toate valorile proprii ale lui A, au părți reale strict negative, toate intrările converg cãtre 0 ca , ceea ce înseamnă că rezultã

converg la rãspunsul forțat

Notație. Matricea este numitã Hurwitz sau matrice stabilã dacã toate valorile proprii sunt strict reale negative.

3.5 Teoria Dunford-Taylor pentru matricea f(A)

Fie o funcție f:olomorfã în domeniul Dse poate scrie sub forma unei integrale cu parametru

(3.39)

dacã este un contur simplu.

Fie acum și A o matrice pãtratã oarecare, A=[ajk], j, k=1, 2,…,ajk Pentru aceastã matrice se vor defini puterile naturale Am, punând

A0=I, Am=A (3.40)

și matricea p(A), unde este o funcție polinomialã, astfel încât dacã p(z)=a0+a1z+…+amzm, a0, …, am, și vom scrie

p(A)=a0I+a1A+…+akAk+…+amAm (3.41)

unde I este matricea unitate, I= , sunt numerele lui Kronecker. Se poate defini și matricea inversã A- matricii A dar numai dacã proprietatea de definiție fiind AA-1=A-1A=I, iar în acest caz au sens și puterile negative ale lui A-m=(A-1)m, m=1,2,… .

Definiția 1. Matricea rezolventã a matricii A, notatã Rz(A) este matricea inversã a matricii (zI-A), adicã

Rz(A)= (zI-A)-1 (3.42)

(dacã aceastã matrice are sens) consideratã cu aplicație

Observația 1. Matricea rezolventã Rz(A) va avea sens numai pentru acei pentru care det(zI-A; det(zI-A) este o funcție polinom de grad n în z, numitã polinomul caracteristic al matricei A, se deduce, notând prin polinomul caracteristic

P(z)=zn+p1zn-1+…+pn-1z+pn=(z-(z- (3.43)

n1+n2+…+np=n, cã matricea Rz(A) va exista pentru orice

Definiția 2. Numerele complexe definite mai sus vor constitui valorile (numerele) proprii matricii A, iar mulțimea lor va forma spectrul matricii A.

Definiția 3. (Dunford-Taylor). Matricea f(A) se definește prin relația

(3.44)

unde f este olomorfã și este un contur simplu din C care conține la interior spectrul matricii A, dar nu conține singularitãti ale lui f.

Exemplul 3.6

Fie A o matrice pãtratã de ordinul doi

atunci ,

Formula Dunford-Taylor se va scrie asfel:

(3.45)

Pentru calculul matricii (3.45) se va considera douã cazuri:

i) Spectrul matricii A este format din puncte distincte: P(z1)=0, P(z2)=0, .

Rezultã f(A)=

= (3.46)

ii) Spectrul matricii A este dublu: z1=z2.

Rezultã f(A)= =

=

(3.47)

Se observã cã în cazul i f(A) se mai poate scrie

f(A)=f(z1)Z1+ f(z2)Z2, Zj= j=1, 2, ()

unde z1, z2, sunt matrici pãtrate de același ordin ca și A, care nu depind de f ci numai de matricea A. Matricele z1, z2, se pot determina particularizând funcțiile f.

De exemplu, luând f(z)=1 și f(z)=z se dedeuce sistemul matriceal z1+z2=I, z1Z1+z2Z2=A, de unde rezultã Z1 și Z2. Se observã ca nu putem avea Z1=0 sau Z2=0. f(A) se va defini prin (3.47) numai dacã f este definitã pe spectrul matricii A.

Aceleași observații sunt valabile și în cazul z1=z2, când f(A) este determinatã prin (3.47) și avem

()

matricele Z1 și Z2 depinzând doar de matricea A (nu depinde de f) și de accea pot fi obținute direct particularizând f în (3.47)’ (se poate lua iarãși f(z)=1 respectiv f(z)=z), iar f(A) poate fi definitã prin (3.47)’ numai dacã este diferențiabilã în vecinãtatea punctului z1 (f putând fi și valori reale, dacã Dacã z1=z2, se poate întâmpla ca matricea Z2 sã fie nulã, iar în acest caz avem

, iar f(A)=.

Exemplul 3.7

Fie A o matrice pãtratã de ordinul trei,

A=;

zI-A=,

Rz(A)=(zI-A)-1=

= (3.48)

iar,

Pentru matricea f(A) se obține egalitatea

f(A)= (3.49)

unde polinomul caracteristic P este definit prin (3.48)’, iar polinoamele pjk sunt determinate prin elementele matricii rezolvente Rz(A) definitã prin (3.48).

Calculând o sã avem

f(A)= ()

integralele din (3.49) calculându-se cu ajutorul teoremei reziduurilor (f nu are singularitãți in interiorul lui unde se vor afla punctele spectrului lui A).

Calculul matricei din (3.49)’ se realizeazã în trei cazuri:

Spectrul matricii A este format din trei puncte distincte z1, z2, z3, , Conform teoremei reziduurilor, se deduce

f(A)=

(3.50)

unde matricele de ordinul al treilea Z1, Z2, Z3 nu depind decât de matricea A:

Zj=, j=1, 2, 3 … ()

O valoare proprie a matricii A este dublã: z1=z2, iar . Teorema reziduurilor conduce la formula

f(A)= (3.51)

termenii notați prin * din ultima matrice fiind asemãnãtori celor scriși, dar cu schimbarea indicilor polinoamelor pjk. Din (3.51) se deduce posibilitatea scrierii matricii f(A) sub forma

f(A)=f(z1)Z10+f’(z1)Z11+f(z2)Z2 ()

în care matricile Z10, Z11 și Z2 sunt de aceleași ordin cu matricea A și nu depind decât de aceastã matrice (nu depinde de f).

Spectrul matricii A este triplu, adicã polinomul caracteristic () se scrie P(z)=(z-z1)3.

Matricea f(A) va avea forma simplã

f(A)=

= (3.52)

unde matricele de ordinul trei Z10, Z11, Z12 nu vor depinde decât de matricea A, nu și de f, iar din (3.48) se poate observa cã

deoarece deg pjj=2, deg dacã j, k=1, 2, 3.

Exemplul 3.8

Fie A=

atunci P(z)=z3-11z2+36z-36, iar spectrul matricii A este {2, 3, 6}.

Se deduce

B(z)=

Se scrie polinomul caracteristic și matricea B sub forma

P(z)=zn+p1zn-1+…+pn-1z+pn, (3.53)

B(z)=B1zn-1+B2zn-2+….+Bn-1z+Bn ()

Mn), j=1, 2,…, n; deoarece va avea loc egalitatea P(z)Rz(A)=B(z) din care, prin înmulțirea la stânga cu (zI-A) se obține

P(z)I=(zI-A)B(z).

Înlocuind P și B prin (3.53) și () și egalând coeficienții acelorași puteri ale lui z se deduc egalitãțile

B1=I, Bj-ABj-1=pj-1, j=2, 3,…, n, -ABn=pnI

din care rezultã

B1=I, B2=A+p1I, B3=A2+p1A+p2I, B4=A3+p1A2+p2A+p3I, …

Aceste egalitãți permit determinarea matricii Rz(A) cu ajutorul coeficienților polinomului caracteristic și al puterilor succesive ale matricii A. Înlocuind matricile B1, B2,… se gãsește o nouã matrice de formã a matricii B:

B(z)=b0(z)I+b1(z)A+…+bn-1(z)An-1 (3.53’’)

unde b0, b1, …bn-1 sunt funcții polinomiale de grad cel mult n-1, n-2,… 0.

Ca o concluzie, pentru orice funcție f definitã pe spectrul matricii A (dacã polinomul minimal al matricii A are numai zerouri simple) sau pentru orice funcție f care admite derivate în numãr convenabil, se poate construi matricea f(A).

Capitolul IV

Metode de calcul pentru determinarea comportãrii dinamice a sistemelor liniare

Capitolul de față este rezervat prezentării unor metode moderne de calcul a comportării dinamice a sistemelor liniare, netede și invariante în timp. Subiectul tratat produce un interes distinct în numeroase aplicații tehnice și științifice ce apelează la modele matematice descrise de ecuații diferențiale liniare.

După cum se va putea vedea, metodele de calcul propuse acoperă și domeniul clasic al determinării răspunsului unui sistem liniar invariant în timp la impuls Dirac sau la o intrare treaptă, din teoria sistemelor automate.

4.1. Metode de calcul a matricei de tranziție

În acest paragraf se vor prezenta principalele metode de calcul numeric a matricei de tranziție, accentul punându-se pe calculul matricei de tranziție asociată sistemelor liniare, invariante în timp.

Fie sistemul liniar descris de

(4.1)

unde A(t) este o matrice de forma , cu toate elementele continue pentru orice t, iar este un vector n-dimensional. Matricea de tranziție asociată sistemului liniar (4.1), este definită prin relațiile

(4.2)

unde I este matricea unitate de forma . Soluția ecuației diferențiale vectoriale (4.1) este dată de

(4.3)

Principalele metode pentru calculul matricei de tranziție în cazul sistemelor liniare variabile în timp sînt următoarele:

-metode bazate pe integrarea numerică a ecuației diferențiale (4.1), de exemplu aplicarea metodei Runge-Kutta;

-metodă bazată pe dezvoltarea Peano-Baker (4.1). Secvența matricelor converge uniform către , unde

(4.4)

Blair (4.2) a stabilit următoarea dezvoltare în serie a matricei de tranziție

(4.5)

unde este operatorul d/dt, I este matricea unitate de forma și (k) semnifică un număr de k operații consecutive ale operatorului [A(t)-DI]. Din nefericire, formula (4.5) este nepotrivită pentru calculul numeric al matricei de tranziție

În concluzie, pentru calculul numeric al matricei se recomandă utilizarea metodelor bazate pe integrarea numerică a ecuațiilor diferențiale.

În cele ce urmează se va trata numai cazul matricei A constante, adică vom considera sistemul liniar

(4.6a)

(4.6b)

Matricea de tranziție asociată sistemului liniar (4.6), notată, poate fi definită prin seria

(4.7)

iar soluția sistemului (4.6) este

(4.8)

Metodele de calcul ale matricei de tranziție (4.7) pot fi clasificate astfel:

(a) metode ce utilizează valorile proprii ale matricei A,

(b) metode ce nu utilizează valorile proprii ale matricei A,

(b1) metode bazate pe integrarea ecuației diferențiale (4.6) ,

(b2) metode bazate pe dezvoltarea în serie trunchiată,

(b3) alte metode.

După prezentarea tuturor acestor metode, în încheierea paragrafului, se vor face unele observații generale privind calculul numeric al matricei de tranziție.

Metode bazate pe cunoașterea valorilor proprii

În cele ce urmează se vor prezenta principalele metode ce utilizează direct sau indirect valorile proprii ale matricei A pentru calculul matricei de tranziție.

Trebuie subliniat că determinarea matricei de tranziție poate fi redusă la calculul matricei de tranziție prin utilizarea relației

unde,

,

iar matricea transformării U este nesingulară. Evident, utilizarea acestei metode se recomandă când se cunoaște matricea transformării U și, în același timp, calculul matricei este, în raport cu calculul direct al matricei de tranziție , mai ușor de efectuat (de exemplu, matricea este forma Jordan sau forma companion matriceal).

Dificultatea determinării valorilor proprii ale unei matrice de ordin superior face ca aceste metode să aibă un domeniu restrâns de aplicare în calculele numerice.

După cum se știe pentru orice matrice A de forma , avem

în care

este forma canonică Jordan a matricei A, iar

Mărimile sunt valorile proprii de multiplicitate , ale matricei A, iar matricele , sunt de forma . Evident avem

unde

iar

În cazul matricei A cu valori proprii distincte (sau, mai general, în cazul matricei A diagonalizabile), coloanele matricei P sunt tocmai vectorii proprii (liniar independenți) ai matricei A. Din acest punct de vedere calculul efectiv al matricei de tranziție se reduce la determinarea valorilor și vectorilor proprii ale matricei A, ceea ce echivalează cu aducerea matricei A la forma Jordan.

Putzer (4.10) a demonstrat o formulă explicită pentru calculul matricei de tranziție în funcție de valorile proprii ale matricei A. Avem relația

(4.9)

unde

k=1, 2, …, n-1

(4.10)

și

– matrice constante de forma

-valorile proprii ale matricei A, care pot să fie și confundate.

Funcțiile scalare satisfac următoarele ecuații diferențiale

(4.11a)

și

i=2, 3, …,n

(4.11b)

Concluzii

Bibliografie

1. Jora B., Popeea C., Barbulea S., Metode de calcul numeric în Automatică. Sisteme liniare. Ed. Enciclopedică, pp 35-63, 1996.

2. Parlett B. N. A Recurrence among the Elements of Functions of Triangular Matrices, Lin. and Its Applic., vol. 14, pp. 117–121, 1976.

3. Ward R. C. Numerical Computation of the Matrix Exponential with Accuracy J. Numer. Anal., vol. 14, pp. 600–610, 1977.

4. Van Loan C. F. The Sensitivity of the Matrix Exponential, J. Numer. Anal., vol. 6, pp. 971–981, 1977.

5. Moler C. B., Van Loan C. F. Nineteen Dubious Ways to Compute the Exponential of a Matrix, Review, vol. 20, pp. 801–836, 1978.

6. Van Loan C. F. A Note on the Evaluation of Matrix Polynomials, IEEE Trans. Automat. Contr., vol. AC-24, pp. 320–321, 1978.

7. Van Loan C. F. Computing Integrals Involving the Matrix Exponential, IEEE Trans. Automat. Contr., vol. AC-24, pp. 395–404, 1978.

8. Valter Olariu – Valeriu Prepelițã, Teoria Distribuțiilor, funcții complexe și aplicații. Ed. Științificã și Enciclopedicã, pp.164-182, 1986.

9. Hespahna J.P., Linear Systems Theory, Press, pp.46-53, 2009.

10. Ionescu V., Lupaș L., Tehnici de calcul în teoria sistemelor. I. Sisteme liniare, Editura Tehnică, 1974.