Carte_metode de Calcul Numeric in Automatic A

8/3/2019 Carte_metode de Calcul Numeric in Automatic A

1/177

B. Jora C. Popeea S. Barbulea

METODE DE CALCULNUMERIC IN

AUTOMATICA

SISTEME LINIARE

EDITURA ENCICLOPEDICABucuresti 1996


2/177

Cuprins

Cap. 1. Rezolvarea ecuatiilor matriciale liniare1.1 Ecuatii matriciale liniare 15

1.2 Rezolvarea ecuatiilor matriciale de tip AX = C 17

1.3 Rezolvarea ecuatiilor Sylvester 20

1.4 Rezolvarea ecuatiilor Liapunov 27

Exercitii 31

Bibliografie 34

Cap. 2. Calculul functiilor de matrici. Exponentialamatriciala

2.1 Functii de matrici 352.2 Calculul functiilor de matrici 46

2.3 Calculul exponentialei matriciale 51

Exercitii 60

Bibliografie 63

Cap. 3. Tehnici de procesare a modelelor sistemiceliniare

3.1 Modele sistemice liniare 65

3.2 Conexiuni 69

3.3 Realizari 74

3.4 Conversii de modele 84

3.5 Algoritmi de calcul polinomial 88

Exercitii 97

1


3/177

2

Bibliografie 100

Cap. 4. Calculul raspunsului n timp al sistemelorliniare

4.1 Raspunsul liber al sistemelor liniare 101

4.2 Raspunsul la intrari liniar generate 103

4.3 Raspunsul la intrari etajate 109

4.4 Raspunsul stationar 112

4.5 Calculul caracteristicilor de frecventa 115

4.6 Raspunsul sistemelor liniare discrete 117

Exercitii 118

Bibliografie 122

Cap. 5. Proceduri de analiza sistemica

5.1 Stabilitatea sistemelor liniare 123

5.2 Descompunerea spectrala 125

5.3 Controlabilitate si observabilitate 130

5.4 Teste elementare de controlabilitate 134

5.5 Forma Hessenberg controlabila 139

5.6 Descompunerea controlabila 153

5.7 Stabilizabilitate si detectabilitate 155

5.8 Realizari minimale 157

5.9 Gramieni de controlabilitate si observabilitate 158

5.10 Echilibrarea sistemelor liniare 162

Exercitii 168

Bibliografie 174

Cap. 6. Proceduri de alocare a polilor

6.1 Formularea problemei de alocare 175

6.2 Proceduri de alocare pentru sisteme cu o singura intrare 184

6.3 Proceduri de alocare pentru sisteme cu mai multe intrari 193

Proceduri de alocare suboptimala 194

Proceduri de alocare robusta 214


4/177


5/177

4


6/177

Cuvant introductiv

Lucrarea de fata reprezinta o prima parte a cursului universitar MetodeNumerice, destinat studentilor anului III al facultatii Automatica si Cal-culatoare din Universitatea POLITEHNICA din Bucuresti si, evident, poa-te fi utila tuturor studentilor de la facultatile de acelasi profil din tara.Desigur, ea poate interesa si pe alti utilizatori potentiali (studenti, cadredidactice, cercetatori, ingineri, economisti, etc.) ai metodelor moderne decalcul n analiza si proiectarea sistemica.

Cititorul ideal al lucrarii trebuie sa posede cunostinte de baza privindmetodele de calcul numeric matricial enumerate n finalul acestui cuvantintroductiv precum si unele informatii sigure asupra principalelor notiunide teoria sistemelor liniare. De asemenea, el trebuie s a aiba obisnuint a uneigandiri limpezi, algoritmice, precum si pasiunea lucrului la calculator.

Desigur, cititorul real, dispunand ntr-un grad mai mare sau mai mic de

calitatile mentionate mai sus, si le poate perfectiona prin studiul lucrarii,confirmand, si pe aceasta cale, viabilitatea principiilor sistemice.

Lucrarea este structurata n 8 capitole. Ea este nsotita de o bibliografiede baza (cu lucrari referite prin utilizarea cifrelor romane), utila n generalpentru abordarea tuturor temelor propuse. Bibliografii specifice selective,dar care pot fi extinse cu usurinta de catre nsusi cititorul interesat princumularea bibliografiilor lucrarilor citate, sunt atasate fiecarui capitol. Deasemenea, lucrarea este nsotita de o lista de notatii uzuale iar cuprinsulpoate constitui, n acelasi timp, indice de notiuni.

Capitolul 1 completeaza cunostintele cititorului privind metodele de re-zolvare a sistemelor de ecuatii liniare cu cazul matricial n care necunos-cutele sunt structurate ntr-o matrice X

Rmn, m > 1, n > 1, fapt care

permite dezvoltarea unor tehnici specifice. In particular, ecuatiile matri-ciale Liapunov vor fi ntalnite aproape n toate capitolele lucrarii.

Capitolul 2 trateaza problema importanta a calculului functiilor de ma-trici. Un accent special este pus pe calculul exponentialei matriciale, uti-

5


7/177

6

lizate intensiv n capitolele 4 si 7.

Capitolul 3 expune principalele tehnici de procesare a modelelor sis-temice liniare (conexiuni, realizari, conversii de modele, algoritmi de calculpolinomial) fara de care nu este posibila abordarea niciunei probleme deanaliza sau sinteza sistemica asistate de calculator.

Capitolul 4 descrie metodele numerice de simulare, i.e. de calcul alraspunsului n timp al sistemelor liniare la stimuli externi caracteristicipentru diverse regimuri de functionare. Aceste tehnici de simulare sunt ab-solut necesare pentru validarea demersului teoretic si calculatoriu de analizasi sinteza sistemica si reprezinta o punte de legatura esentiala catre lumeareala a experimentului fizic.

Capitolul 5 prezinta procedurile numerice de analiza a proprietatilorsistemice fundamentale (stabilitate, controlabilitate, observabilitate, etc.)precum si principalele metode de constructie a realizarilor minimale siechilibrate. In acest fel continutul acestui capitol se constituie ntr-o colectiede algoritmi fundamentali, componente de baza ale oricarei activitati deconceptie asistate de calculator a unor sisteme de conducere complexe.

Capitolul 6 realizeaza o trecere n revista a principalelor proceduri nu-merice de alocare a polilor sistemelor liniare. In contextul general al tehni-cilor de sinteza, alocarea polilor reprezinta un instrument necesar pentruasigurarea unei dinamici dorite sistemelor automate elaborate. Daca ncazul sistemelor simple reactia dupa stare ce realizeaza alocarea este unic

determinata si, deci, libertatile de calcul sunt strict procedurale, n cazulsistemelor multiple apar posibilitati suplimentare de optimizare sau robus-tificare a spectrului alocat.

Capitolul 7 abordeaza, din punctul de vedere al metodelor numericede calcul, problema de sinteza liniar-patratica, care vizeaza optimizareaunui sistem liniar cu indice de calitate patratic. In aplicatii problemaliniar-patratica generala apare sub forma unor probleme cu finalitati sau cutehnici de tratare specifice. In acest context sunt prezentate metodele decalcul cele mai performante pentru rezolvarea ecuatiilor matriciale Riccati,nsotite de analize comparative ale eficientei si stabilitatii lor numerice.

Capitolul 8 considera unele generalizari ale problemei de sinteza liniar-patratice si ilustreaza aplicarea algoritmilor de calcul stabiliti la constructia

compensatoarelor H2 si H (sub)optimale.

Principalele rezultate ale expunerii sunt concretizate sub forma de al-goritmi de calcul direct implementabili iar fiecare capitol este nsotit de unset de probleme.


8/177

7

Ne exprimam opinia ca pentru nsusirea materialului prezentat este ab-solut necesara rezolvarea problemelor propuse si, mai ales, implementareaalgoritmilor, urmata de experimentarea lor pe exemple numerice concretesi semnificative.

Autorii multumesc colegului lor prof. Paul Flondor pentru comentariileconstructive facute pe marginea lucrarii. De asemenea, multumesc domnu-lui Marcel Popa, directorul Editurii Enciclopedice, pentru atentia plina desolicitudine cu care a urmarit aparitia prompta si n cele mai bune conditii alucrarii. Autorii apreciaza n mod deosebit interesul manifestat de domnuling. Sabin Stamatescu, directorul firmei ASTI Bucuresti, fata de continutullucrarii si posibilele aplicatii ale procedurilor de calcul prezentate n con-ducerea cu calculator a proceselor industriale. In sfarsit, dar, desigur, nu n

ultimul rand, multumiri se cuvin studentilor Siana Petropol, Simona Chir-vase, Stefan Murgu, Catalin Petrescu, Simona Ruxandu, Adrian Sandu,Laura Stan, Florin Roman, Madalina Zamfir, Ion Vasile care au realizatredactarea computerizata a lucrarii cu mentiunea expresa ca fara muncalor plina de daruire nimic din ce am realizat mpreuna nu ar fi fost posibil.

Enumeram n continuare principalele probleme de calcul numeric ma-tricial, mpreuna cu succinte referiri la modul de rezolvare al acestora, pecare cititorul interesat de cuprinsul lucrarii de fata este bine sa le aiba nvedere.

1. Pentru rezolvarea unui sistem de ecuatii liniare Ax = b, n care ma-tricea A Rnn este inversabila, se utilizeaza procedura de triangularizareprin eliminare gaussiana cu pivotare partiala

A R = MA, b Mb,unde R rezulta superior triunghiulara inversabila iar M este o secventa detransformari elementare stabilizate, urmata de rezolvarea prin substitutienapoi a sistemului superior triunghiular rezultat. Efortul de calcul nece-sar este Nop n3/3.

In cazul simetric A = AT, se recomanda conservarea acestei proprietatiprin utilizarea schemei de eliminare simetrizate A D = M AMT, propusede Parlett [ III ], n care D rezulta cvasidiagonala.

In cazul simetric pozitiv definit A = AT > 0, se impune aplicarea fac-torizarii Cholesky A = LLT.

2. Pentru rezolvarea, n sensul celor mai mici patrate (CMMP), asistemului supradeterminat Ax = b, unde matricea A Rmn este monica,


9/177

8

i.e. are coloanele liniar independente (rangA = n), se utilizeaza procedurade triangularizare ortogonala

A R = U A, b U b,unde R rezulta superior triunghiulara iar U este o secventa de transformariortogonale elementare (reflectori sau rotatii). Triangularizarea ortogonalaeste, evident, echivalenta cu factorizarea (sau descompunerea) A = QR,unde factorul ortogonal Q = UT se obtine n forma factorizata. Efortulde calcul necesar este Nop mn2 n3/3 iar eventuala formare explicitaa matricii Q, prin acumularea transformarilor elementare partiale, necesitanca 2n3/3 operatii n cazul cel mai defavorabil (m = n).

Daca matricea A nu este neaparat monica, atunci, pentru asigurareastabilitatii numerice este necesara utilizarea unei strategii adecvate de piv-otare a coloanelor, i.e. A R = U AP, unde P este o matrice de permutareiar R rezulta superior trapezoidala, vezi [III].

3. Pentru rezolvarea, n sens CMMP, a sistemului subdeterminatAx = b, unde matricea A Rmn este epica, i.e. are liniile liniar indepen-dente (rangA = m), se utilizeaza procedura de triangularizare ortogonala ladreapta

A L = AV,unde L rezulta inferior triunghiulara. (Aceasta procedura este echivalentacu procedura precedenta reformulata prin dualitate, i.e. aplicata lui AT.)Efortul de calcul necesar este Nop

nm2

m3/3.

4. Pentru calculul valorilor proprii ale unei matrici A Rnn seutilizeaza algoritmulQR care, n esenta, construieste iterativ forma SchurS a lui A, i.e.

A S = QTAQ, (Nop 6n3),unde S rezulta cvasisuperior triunghiulara iar Q este o secventa de trans-formari ortogonale. In caz de necesitate valorile proprii pot fi ordonatepe diagonala principala a lui S, conform oricarui criteriu de ordonare im-pus, cu ajutorul unei secvente suplimentare de transformari ortogonale deasemanare. Calculul matricii Q prin acumularea transformarilor partiale(necesara, e.g. pentru calculul vectorilor proprii) implica aproape dublareanumarului de operatii mentionat.

In cazul simetric A = A

T

, se utilizeaza avantajos versiunea simetrica aalgoritmului QR, i.e. A = QTAQ, unde rezulta diagonala.5. Pentru rezolvarea problemei CMMP generale n care matricea A

Rmn este de rang nu neaparat maximal r def= rangA min(m, n), precum


10/177

9

si a altor probleme importante de calcul numeric matricial (determinarearangului, a pseudoinversei, operatii cu subspatii liniare, etc.) se utilizeazadescompunerea valorilor singulare DVS

UTAV = ,

unde matricea contine pe diagonala principala valorile singulare nenule1 2 . . . r > 0 ale lui A (restul elementelor fiind nule) iar ma-tricile U, V sunt ortogonale. (Algoritmul DVS utilizat pentru calcululdescompunerii de mai sus constituie o adaptare ingenioasa a algoritmu-lui QR simetric). Din nou, formarea factorilor ortogonali U, V necesitaacumularea (relativ costisitoare) a transformarilor partiale.

6. Pentru calculul valorilor proprii generalizate ale unui fascicol BA,unde A, B Rnn, se utilizeaza algoritmul QZ care, n esenta, aduceperechea (A, B) la forma Schur generalizata

A S = QTAZ, B T = QTAZ,

unde S rezulta cvasisuperior triunghiulara, T rezulta superior triunghiularaiar Q si Z sunt secvente de transformari ortogonale.

In cazul simetric A = AT, B = BT > 0 se utilizeaza factorizareaCholesky B = LLT si se aplica algoritmul QR simetric matricii A =L1ALT, ceea ce conduce la diagonalizarea simultana a ambelor matriciA, B.

Pentru detalii privind aspectele teoretice si procedurale ale problemelorde calcul matricial enumerate mai sus recomandam consultarea referintelorbibliografice [V], [VI] sau [IX], [X] .


11/177

10

Bibliografie

Pentru programe de calcul si indicatii de utilizare:

[ I ] Smith B.T., Boyle J.M., Ikebe Y., Klema V.C., Moler C.B. Ma-trix Eigensystem Routines: EISPACK Guide, 2-nd ed., Springer-Verlag, New York, 1970.

[ I I ] Garbow B.S., Boyle J.M., Dongarra J.J., Moler C.B. MatrixEigensystem Routines: EISPACK Guide Extension, Springer - Ver-lag, New York, 1972.

[III ] Dongarra J.J., Bunch J.R., Moler C.B., Stewart G.W. LINPACK

Users Guide, SIAM Publications, Philadelphia, PA, 1978.[ I V ] Moler C.B., Little J.N., Bangert S. PC-MATLAB Users Guide,

The Math Works Inc., 20 N. Main St., Sherborn, Mass., 1987.

Pentru algoritmi de calcul matricial:

[ V ] Stewart G. W. Introduction to Matrix Computations, AcademicPress, New York and London, 1973.

[ V I ] Golub G. H., Van Loan Ch. F. Matrix Computations, Secondedition, The Johns Hopkins University Press, Baltimore and London, 1988.

Pentru chestiuni teoretice de calcul matricial:

[VII] Gantmaher F.R. Teoriia matrit (editia a 2-a), Ed. Nauka, Moscova,1966. (The Theory of Matrices, vols. 1-2, Chelsea, New York, 1959).

Lucrari n limba romana:

[VIII] Ionescu V., Lupas L. Tehnici de calcul n teoria sistemelor,vol.1. Sisteme liniare, vol.2. Sisteme optimale, E.T., Bucuresti, 1973.

[ I X ] Bucur C.M., Popeea C.A., Simion Gh.Gh. Matematici speciale.Calcul numeric, E.D.P., Bucuresti, 1983.

[ X ] Ionescu V., Varga A. Teoria sistemelor. Sinteza robusta. Me-tode numerice de calcul., Ed. ALL, Bucuresti, 1994.

[ X I ] Iorga V., Jora B., Nicolescu C., Lopatan I., Fatu I., Programarenumerica, Ed. Teora, Bucuresti, 1996.

Alte titluri de uz general:


12/177

11

[XII] Astrom K.J., Wittenmark B. Computer Controlled Systems,Prentice Hall Inc., Englewood Cliffs, NJ, 1984.

[XIII] Jamshidi M., Herget C.J. (ed.), Computer-Aided Control Sys-tems Engineering, North-Holland, Amsterdam, 1985. (Masinostroienie,Moscova, 1989.)

[XIV] Solodovnicov V.V. (ed.), Avtomatizirovannoe proektirovaniesistem upravleniia, Masinostroienie, Moscova, 1990.

[ X V ] Jamshidi M., Tarokh M., Shafai B. Computer-Aided Analysisand Design of Linear Control Systems, Prentice Hall Inc., Engle-wood Cliffs, NJ, 1992.

Lista de notatii

Notatii generaleN multimea numerelor naturale.Z multimea numerelor ntregi.R multimea numerelor reale.C multimea numerelor complexe.

C

semiplanul stang deschis al planului variabilei complexe s.

D1(0) discul unitate deschis al planului variabilei complexe z.Rn spatiul liniar n-dimensional al vectorilor (coloana) x cu n componente

reale xi R, i = 1 : n.Cn spatiul liniar n-dimensional al vectorilor (coloana) x cu n componente

complexe xi C, i = 1 : n.ek, k = 1 : n baza standard a spatiului liniar Rn, respectiv Cn.Rmn spatiul liniar al matricilor cu m linii si n coloane cu elemente reale

aij R, i = 1 : m, j = 1 : n.

Cmn spatiul liniar al matricilor cu m linii si n coloane cu elementecomplexe aij C, i = 1 : m, j = 1 : n. 1

1 In calcule, vectorii se identifica cu matricile cu o singura coloana iar scalarii seidentifica cu matricile (sau vectorii) cu un singur element.


13/177

12

AT transpusa matricii (reale sau complexe) A.

AH conjugata hermitica a matricii (complexe) A, i.e. AH = AT, unde Aeste conjugata complexa a lui A.

A+ pseudoinversa normala (Moore-Penrose) a matricii A; daca A este

monica A+ = (ATA)1

AT, daca A este epica atunci A+ = AT(AAT)1

.

i(A), i = 1 : p, p = min(m, n) valorile singulare ale matricii A ordonateastfel ncat 1 2 . . . p.

(A) multimea {1(A), 2(A), . . . , p(A)} a valorilor singulare ale ma-tricii A.

r = rangA rangul matricii A, i.e. numarul valorilor singulare nenule.

In matricea unitate de ordinul n.

A1 inversa matricii patrate nesingulare A, i.e. AA1 = A1A = In.

AT = (A1)T

= (AT)1

AH = (A1)H

= (AH)1

trA urma matricii patrate A, i.e. suma elementelor diagonale.

detA determinantul matricii patrate A.

i(A), i = 1 : n valorile proprii ale matricii patrate A de ordin n.

(A) spectrul (de valori proprii) {1(A), 2(A), . . . , n(A)} al matricii A.(A) = maxi=1:n|i(A)| raza spectrala a matricii A.cond(A) = A A1 numarul de conditie la inversare al matricii A (

este o norma matriciala consistenta, vezi mai jos).

(x, y) = yTx produsul scalar standard a doi vectori reali; n cazul complexprodusul scalar este (x, y) = yHx.

x = (x, x)1/2 norma euclidiana a vectorului x; se noteaza x2Q = xTQxunde Q este o matrice simetrica (Q = QT).

xp = (ni=1 |x|p)1/p p-normele vectorului n-dimensional x, p 1; ncalcule se utilizeaza n special x1,x2 =x si x=maxi=1:n|xi|.

(A, B) = tr(BTA) (tr(BHA)) produsul scalar a doua matrici reale (com-plexe).


14/177

13

AF = (A, B)1/2 norma Frobenius a matricii A,A2F =

mi=1

nj=1 |aij|2 sau A2F =

ri=1 i

2.

|A|p = (r

i=1 ip)

1/p p-normele Schatten, p 1; n calcule se utilizeaza

n special norma-urma |A|1 =r

i=1 i, norma Frobenius |A|2 = AFsi norma spectrala |A| = 1(A).

Ap = maxxp=1Axp p-normele induse; n calcule se utilizeaza nspecial norma A1 = maxj=1:n

mi=1 |aij |, norma spectrala A2 =

1(A) si norma A = maxi=1:mn

j=1 |aij |.

Notatii sistemice generalex Rn vectorul (marimilor) de stare.

u Rm vectorul (marimilor) de intrare.

y Rl vectorul (marimilor) de iesire.

S = (A,B,C,D) reprezentarea de stare a unui sistem liniar (A Rnn,B Rnm, C Rln si D Rlm).

T = (N, p) reprezentarea de transfer a unui sistem liniar (N este ma-tricea coeficientilor polinoamelor ce definesc numaratorul matricii detransfer iar p este vectorul coeficientilor polinomului numitor comuni.e. T(s) = C(sI A)1B + D = N(s)/p(s)).

(t) ((k)) impulsul unitate continuu (discret).

1(t) (1(k)) functia treapta unitara continua (discreta).

Transformari(1) MAN (MAN1 sau MANT) transformare de echivalent a (bilate-

rala) a matricii A Rmn (M si N sunt matrici patrate nesingulare;transformarea de echivalenta conserva rangul iar daca M, N suntortogonale atunci conserva si valorile singulare).

(2) N AN1 transformare de asemanare a matricii A Rnn (transfor-marea de asemanare conserva valorile proprii).


15/177

14

(3) N ANT transformare de congruent a a matricii A Rnn (N estenesingulara; aplicata unei matrici A simetrice, transformarea de con-gruenta conserva rangul si inertia i.e. numerele de valori proprii neg-ative, nule si, respectiv, pozitive).

Daca N este ortogonala atunci transformarile (2) si (3) coincid si definesctransformarea de asemanare ortogonala.

PrescurtariSISO sigla pentru sisteme simple avand o singura intrare si o singura

iesire (Single-Input Single-Output).

SIMO, (MISO, MIMO) sigla pentru sistemele cu o singura intrare simai multe iesiri. Celelalte prescurtari au semnificatii evidente.

FSR(G) forma Schur reala (generalizata).

FSC(G) forma Schur complexa (generalizata).

DVS descompunerea valorilor singulare.

FSH forma (bloc-)superior Hessenberg.

EM(A)L ecuatie matriciala (algebrica) Liapunov continua.

DEM(A)L ecuatie matriciala (algebrica) Liapunov discreta.

EM(A)R ecuatie matriciala (algebrica) Riccati continua.

DEM(A)R ecuatie matriciala (algebrica) Riccati discreta.


16/177

Capitolul 1

Rezolvarea ecuatiilormatriceale liniare

Acest capitol este consacrat tehnicilor de rezolvare a unor sisteme liniare, ngeneral de mari dimensiuni, structurate n exprimari matriceale care permitdezvoltarea unor metode specifice de calcul.

1.1 Ecuatii matriceale liniare

Ecuatiile matriceale liniare sunt sisteme de ecuatii liniare care pot fi scrisecompact ntr-o forma matriceala de tipul

f(X) = C (1.1)

unde X Rmn este matricea necunoscutelor, C Rpq o matrice data,iar f : Rmn Rpq o aplicatie liniara (i.e. f satisface f(X1 + X2) == f(X1) + f(X2) si f(X1) = f(X1) pentru orice matrice X1, X2 Rmn si orice scalar ). Se poate arata ca orice aplicatie liniara de argu-ment matriceal poate fi scrisa sub forma

f(X) =k

i=1

AiXBi, Ai Rpm, Bi Rnq (1.2)

pentru un anumit k si, n consecinta, (1.1) devine

ki=1

AiXBi = C. (1.3)

15


17/177

16 CAPITOLUL 1. ECUATII MATRICEALE LINIARE

Evident, ecuatia (1.3) poate fi scrisa ntr-o forma desfasurata ca sis-tem de pq ecuatii liniare cu mn necunoscute care poate servi ca baza pentruexprimarea conditiilor de existenta/unicitate ca si pentru elaborarea pro-cedurilor de rezolvare ntr-un sens sau altul. O astfel de abordare eludeazastructura matriceala a sistemului iar aplicarea tehnicilor clasice de re-zolvare pe sistemul desfasurat, neexploatand structura interna a datelor deintrare, este, de cele mai multe ori, ineficienta (vezi exercitiul 1.9).

Din acest motiv, pentru diferite cazuri particulare ale ecuat iei matriceale(1.3), conditiile de existenta/unicitate ale solutiilor, ntr-un sens precizat, seexprima n raport cu matricele Ai, Bi, i = 1 : k, iar metodele de rezolvarefac apel la tehnici speciale.

Ecuatiile matriceale liniare (1.3) cele mai ntalnite se obtin pentru k = 1

AXB = C (1.4)

a carei rezolvare se reduce imediat la rezolvarea unor ecuat ii matricealeavand forma particulara

AX = C (1.5)

XB = C (1.6)

respectiv, pentru k = 2,

A1XB1 + A2XB2 = C (1.7)

cu particularizarile

AX+ XB = C, (1.8)AXB + X = C. (1.9)

Ecuatiile (1.8), (1.9) sunt cunoscute sub denumirea de ecuatii matricealetip Sylvester. Datorita frecventei utilizari a ecuatiilor matriceale de forma(1.8) n teoria sistemelor dinamice continue, respectiv a ecuatiilor de forma(1.9) n domeniul sistemelor dinamice discrete, n continuare ecuatia (1.8) vafi referita ca ecuatie Sylvester continua, respectiv (1.9) ca ecuatie Sylvesterdiscreta.

In sfarsit, considerand n (1.8) B A, A AT si n (1.9) B A,A AT, C C, obtinem ecuatiile matriceale liniare cunoscute subdenumirile de ecuatie Liapunov continua pentru

ATX+ XA = C (1.10)

respectiv ecuatie Liapunov discreta pentru

ATXA X = C. (1.11)


18/177

1.2. ECUATII DE FORMA AX = C 17

In cele mai multe aplicatii matricele de intrare A, B, C sunt reale, ma-tricea solutie X rezultand, la randul sau, reala. Totusi, algoritmii prezentatiraman aplicabili si n cazul matricelor de intrare complexe. Mai mult, chiarn cazul unor date initiale reale, unele din metode de calcul, apeland ladeterminarea valorilor si vectorilor proprii, necesita, tranzitoriu, utilizareaunei aritmetici complexe. Atunci cand pentru date de intrare reale se poateutiliza exclusiv o aritmetica reala (vezi exercitiul 1.14) vor fi facute pre-cizarile cuvenite.

In cele ce urmeaza vom fi interesati n calculul solutiilor ecuatiilor ma-triceale liniare cu semnificatie pentru domeniul teoriei sistemelor n general,si al sistemelor automate n special, respectiv al sistemelor de forma (1.5),(1.6), (1.8), (1.9) si (1.10), (1.11). Conditiile de existenta si unicitate ale

solutiilor se exprima n mod specific pentru fiecare tip de ecuatie, singurulrezultat comun fiind dat de

Propozitia 1.1 Daca ecuatia matriceala liniara (1.3) admite o solutie X Rmn, atunci aceasta solutie este unica daca si numai daca ecuatia ma-triceala omogena (obtinuta din ecuatia initiala pentru C = 0) admite dreptunica solutie X = 0.

Demonstratia este propusa ca exercit iu pentru cititor sau poate fi gasita n[VII].

In cazul ecuatiilor matriceale liniare (1.8), (1.9) si (1.10), (1.11) o analizasimpla arata ca matricea necunoscuta X si membrul drept C au aceleasidimensiuni (m

n). In consecinta, enuntul de mai sus capata urmatoarea

forma mai precisa.

Propozitia 1.2 Ecuatiile matriceale liniare (1.8), (1.9) si (1.10), (1.11)sunt global solubile, i.e. admit o solutie X Rmn oricare ar fi membruldrept C, daca si numai daca ecuatiile matriceale omogene corespunzatoareadmit drept unica solutie X = 0.

Altfel spus, are loc urmatoarea alternativa (a lui Fredholm): fie ecuatiileconsiderate sunt global solubile, fie ecuat iile omogene admit o solutie ne-triviala X = 0.

1.2 Rezolvarea ecuatiilor matriceale de tip AX=

C

Ecuatiile matriceale de tipul AX = C, unde A Rpm, X Rmn,C Rpn,reprezinta o colectie de ecuatii vectoriale obisnuite, toate avand


19/177


aceeasi matrice A a coeficientilor. Intr-adevar, partitionand matricele X siC pe coloane

X = [ x1 x2 . . . xn ], C = [ c1 c2 . . . cn ] (1.12)

cu xj = Xej , cj = Cej , j = 1 : n, ecuatia matriceala AX = C devineechivalenta cu setul de sisteme liniare

Axj = cj , j = 1 : n (1.13)

care se rezolva cu mijloacele clasice, cu precizarea ca transformarile nece-sare asupra matricei A a sistemelor (1.13) (cum ar fi, de exemplu, n cazul

p = m si A inversabila, triangularizarea prin eliminare gaussiana cu piv-

otare partiala) se efectueaza o singura data. Avand n vedere ca rezolvareasistemelor (1.13) nu ridica probleme deosebite, scrierea algoritmilor faceobiectul exercitiilor 1.1 si 1.3.

Observatia 1.1 Ecuatiile de tipul XB = C se reduc la ecuatii de tipulAX = C prin transpunere.

Observatia 1.2 Daca matricea A Rpm este monica(i.e. are coloaneleliniar independente) iar sistemele (1.13) sunt rezolvate n sensul celor maimici patrate (CMMP), atunci ansamblul pseudosolutiilor

xj = A+cj = (A

TA)1

ATcj , j = 1 : n (1.14)

respectiv

X = A+C, A+ = (ATA)1AT (1.15)

are o semnificatie asemanatoare si anume aceea ca matricea X Rmnminimizeaza norma Frobenius a reziduului matriceal

R = C AX (1.16)adica

RF = C AXF = minXRmn

C AXF. (1.17)

Intr-adevar,

RF = n

j=1

pi=1

r2ij =

nj=1

cj Axj2

este minima daca si numai daca cj Axj sunt minime pentru toti j =1 : n.


20/177

1.2. ECUATII DE FORMA AX = C 19

Observatia 1.3 Daca matricea A Rpm este epica (adica are lini-ile liniar independente), atunci ansamblul solutiilor normale ale sistemelor(1.13)

xj = A+cj = A

T(AAT)1

cj , j = 1 : n (1.18)

respectiv

X = A+C, A+ = AT(AAT)1

(1.19)

reprezinta solutia de norma Frobenius minima a sistemului AX = C :

XF = minAX=C

XRmn

XF. (1.20)

Intr-adevar, este evident ca X dat de (1.19) este o solutie a ecuatiei ma-triceale AX = C. In plus, daca A = LQ este o factorizare a matricei A,unde Q Rmm este ortogonala iar L Rpm are structura

L = [ L 0 ] (1.21)

cu L Rpp inferior triunghiulara nesingulara, atunci ecuatia matricealaAX = C se scrie

LQX = C. (1.22)

Cu notatia

QX = Y = YY , (1.23)cu Y Rpn, Y R(mp)n, avand n vedere structura (1.21) a matriceiL, rezulta

Y = (L)1C.

In consecinta, toate solutiile sistemului AX = C se scriu sub forma

X = QT

(L)1C

Y

(1.24)

cu Y R(mp)n arbitrar. Dar

X2F = QT (L)1C

Y 2F = (L)1CY 2F =

= (L)1C2F + Y2F


21/177


minimul obtinandu-se pentru Y = 0. Prin urmare solutia de norma Frobe-nius minima este

X = QT

(L)1

C0

(1.25)

verificandu-se imediat, pe baza factorizarii LQ, ca aceasta coincide cusolutia (1.19). Observatia 1.4 In cazul n care matricea A este de rang nemaximal

rang A < min(p,m) ,

notand cu xj , j = 1 : n, pseudosolutiile normale, n sensul celor mai micipatrate, ale sistemelor (1.13), matricea X = (xj)j=1:n poate fi interpre-

tata drept pseudosolutia normala n sensul normei Frobenius a ecuatieimatriceale AX = C. Cu alte cuvinte X Rmn este matricea de normaFrobenius minima dintre toate matricele X Rmn care minimizeazanorma Frobenius a reziduului matriceal R = C AX :

XF = minXRmn

CAXF=minim

XF. (1.26)

Justificarea afirmatiilor de mai sus face obiectul exercitiului 1.2.

1.3 Rezolvarea ecuatiilor matriceale

SylvesterI. Vom ncepe cu prezentarea modalitatilor de rezolvare a ecuatiei ma-triceale Sylvester continue (1.8)

AX+ XB = C

unde n general A Cmm, B Cnn, C Cmn. Conditiile de existentasi unicitate ale solutiei acestei ecuatii sunt date de

Propozitia 1.3 Ecuatia matriceala (1.8) admite o solutie X Cmn unicadaca si numai daca

i + j = 0 (1.27)oricare ar fi i

(A) si oricare ar fi j

(B) sau, altfel spus, daca si

numai daca(A) (B) = . (1.28)

Daca matricele A, B, C sunt reale iar conditiile (1.27) sunt satisfacuteatunci solutiaX rezulta si ea real a.


22/177

1.3. ECUATII SYLVESTER 21

Demonstratie. Exista matricele unitare U Cmm si V Cnn (i.e.UHU = U UH = Im, V

HV = V VH = In, indicele superior H semnificanddubla operatie de transpunere si conjugare) astfel ncat

UHAU = S Cmm, (1.29)

VHBV = T Cnn (1.30)sunt formele Schur complexe ale matricelor A respectiv B si, deci, au ostructura superior triunghiulara. Din (1.29), (1.30) avem relatiile

A = U SUH (1.31)

B = V T VH (1.32)

cu care ecuatia (1.8) devine

U SUHX+ X V T V H = C

de undeSUHXV + UHXV T = UHCV . (1.33)

Notand cuY = UHXV , C = UHCV (1.34)

(1.33) se scrieSY + Y T = C. (1.35)

Avand n vedere natura nesingulara a transformarilor efectuate, ecuatiamatriceala (1.8) admite o solutie unica X Cmn daca si numai dacaecuatia (1.35) admite o solutie unica Y Cmn. Conditiile de existentasi unicitate ale solutiei ecuatiei (1.35) rezulta din aplicarea unei tehnici derezolvare directa. Intr-adevar, tinand seama de structura superior triunghi-ulara a matricelor S si T rezulta ca, scriind pe coloane, (1.35) devine

Syj + Y tj = cj, j = 1 : n (1.36)

unde yj = Y ej , tj = T ej , cj = Cej . Dar

tj = [ t1j t2j . . . tjj 0 . . . 0 ]T

si, prin urmare, (1.36) devine

Syj +

jk=1

tkjyk = cj , j = 1 : n (1.37)


23/177


sau

(S+ tjjIm)yj = cj j1k=1

tkjyk, j = 1 : n. (1.38)

Se vede acum limpede ca ecuatiile (1.38) sunt sisteme liniare triunghi-ulare avand membrul drept calculabil daca ordinea de rezolvare a acestorsisteme este j = 1, 2, 3,...,n. Mai mult, existenta si unicitatea solutiei esteconditionata de nesingularitatea matricelor (S + tjjIm), j = 1 : n. Ma-tricea S fiind superior triunghiulara, aceasta conditie este satisfacuta dacasi numai daca

sii + tjj

= 0,

i = 1 : m,

j = 1 : n,

conditie identica cu (1.27) ntrucat (A) = (S) si (B) = (T), observatiecare ncheie demonstratia propozitiei.

Tehnicile de rezolvare ale ecuatiei Sylvester continue fac apel la trans-formari care conduc la ecuatii matriceale n care, n locul matricelor Asi B, apar matrice cu o structura simplificata (Hessenberg, triunghiulara,diagonala) fapt care permite aplicarea metodelor consacrate de la sistemeleclasice. O prima alternativa, numita ad-hoc varianta Schur-Schur, arela baza procedura evidentiata de demonstratia propozitiei 1.3. In cazuluzual, n care datele de intrare precum si solutia X sunt reale, redactareaalgoritmului de calcul este urmatoarea.

Algoritmul 1.1 (Date matricele A Rmm, B Rnn,C Rmn cu (A) (B) = , algoritmul calculeaza solutiaX Rmn a ecuatiei Sylvester continue AX + XB = C prinaducerea matricelor A si B la forma Schur complexa prin trans-formari unitare de asemanare. Algoritmul utilizeaza functia fsccare calculeaza forma Schur complexa si matricea unitara detransformare corespunzatoare pentru o matrice de intrare data.)

1. [ U, S] = fsc(A)

2. [ V, T ] = fsc(B)

3. C C = UHCV4. Pentru j = 1 : n

1. Daca j > 1 atunci1. cj = cj

j1k=1 tkjyk.

2. Se rezolva sistemul superior triunghiular(S+ tjjIm)yj = cj .


24/177


5. X = Re(U Y VH)

Comentarii. Referitor la algoritmul de mai sus se cuvin urmatoareleprecizari:

1. Aducerea unei matrice reale la forma Schur complexa prin trans-formari de asemanare se face n doua etape: ntr-o prima etapa cu algo-ritmul QR se obtine forma Schur reala dupa care, ntr-o a doua etapa,utilizand rotatii complexe sau reflectori complecsi (exercitiile 1.4 : 1.7)blocurile 2 2 sunt reduse la forma superior triunghiulara complexa.

2. Efectuarea calculelor n format virgula mobila cu numere complexeconduce la un rezultat afectat de erori de rotunjire complexe desi rezultatulteoretic este o matrice reala. De aceea, n ultima instructiune, este inclusaeliminarea partilor imaginare.

3. Evident, efortul de calcul cel mai important se consuma n execut iainstructiunilor 1 si 2 de aducere la forma Schur complexa a matricelor A siB si de acumulare a transformarilor. Din motive de eficienta se impune analiza alternativelor n care se renuntala aducerea la forma Schur a ambelor matrice A si B.

Astfel, n asa numita varianta Hessenberg-Schur numai matricea Beste adusa la forma Schur complexa aparand urmatoarele diferente n raportcu algoritmul de mai sus:

In instructiunea 1 matricea A este adusa printr-un algoritm de calculdirect (neiterativ) algoritmul HQ la forma superior Hessenberg. In

acest fel se evita faza iterativa a algoritmului QR si trecerea de la formaSchur reala la forma Schur complexa.

In compensatie, la instructiunea 4.2, se rezolva, n loc de n sistemetriunghiulare, n sisteme de tip Hessenberg incluzand eliminare gaussianacu eventuala pivotare.

Scrierea explicita a algoritmului face obiectul exercitiului 1.10.

In cazul matricelor A si B simple (adica al matricelor pentru careexista un set complet de vectori proprii liniar independenti, n particular,al matricelor cu valori proprii distincte) metoda bazata pe calculul valorilorsi vectorilor proprii poate conduce la o precizie superioara a rezultatelor

simultan cu un timp de executie rezonabil. Fie

(A) = {1, 2, . . . , m} C(B) = {1, 2, . . . , n} C (1.39)


25/177


spectrele de valori proprii ale matricelor A respectiv B si

U = [ u1 u2 . . . um ],V = [ v1 v2 . . . vn ]

(1.40)

matricele corespunzatoare de vectori proprii (adica uj = U ej este un vectorpropriu al matricei A asociat valorii proprii j (A), etc.). Atunci

A = U LU1, L = diag(1, 2, . . . , m)B = V MV1, M = diag(1, 2, . . . , n)

(1.41)

si ecuatia Sylvester (1.8) se scrie sub forma

U LU1X+ X V M V 1 = C (1.42)

echivalenta cuLU1XV + U1XV M = U1CV . (1.43)

NotandU1XV = YU1CV = Z

(1.44)

din (1.43) rezulta ecuatia

LY + Y M = Z (1.45)

cu necunoscuta matriceala Y ale carei elemente, n condit iile propozitiei1.1, se calculeaza imediat cu relatiile

yij = zij/(i + j) , i = 1 : m, j = 1 : n. (1.46)

Matricea necunoscuta X din ecuatia initiala rezulta din (1.44)

X = U Y V1. (1.47)

Presupunand ca dispunem de o functie numita vvp care, pentru o ma-trice de intrare simpla T, furnizeaza o matrice V ale carui coloane sunt vec-torii proprii ai matricei T si un vector p al valorilor proprii corespunzatoricu sintaxa

[ V, p ] = vvp(T) (1.48)

putem prezenta urmatorul algoritm.

Algoritmul 1.2 (Date matricele simple A R

mm, BRnn cu (A) (B) = si matricea C Rmn, algoritmul

calculeaza solutia X Rmn a ecuatiei matriceale AX+ XB =C pe baza calculului valorilor si vectorilor proprii cu ajutorulfunctiei vvp.)


26/177


1. [ U, l ] = vvp(A)

2. [ V, m ] = vvp(B)

3. Se rezolva sistemul matriceal U Z = CV cu necunoscutaZ.

4. Pentru i = 1 : m

1. Pentru j = 1 : n

1. yij = zij/(li + mj).

5. Se rezolva sistemul matriceal XV = U Y n raport cu X.

6. X = Re(X).

Rezolvarea sistemelor matriceale de la instructiunile 3 si 5 ale algorit-mului se realizeaza conform recomandarilor sectiunea 1.2.

II. Modalitatile de calcul ale solutiei ecuatiei matriceale Sylvester dis-crete (1.9)

AXB + X = C

unde A Cmm, B Cnn, C Cmn sunt asemanatoare. Conditiile deexistenta si unicitate ale solut iei sunt date de

Propozitia 1.4 Ecuatia matriceala (1.9) are o solutie X Cmn unicadaca si numai daca

ij = 1 (1.49)oricare ar fi i (A) si oricare ar fi j (B). Daca matricele A, B, Csunt reale iar conditiile (1.49) sunt satisfacute atunci solutia X rezulta siea reala.

Demonstratie. Fie U Cmm si V Cnn matrice unitare care definescreducerea matricelor A si B la formele Schur complexe S, respectiv T, cain (1.29) si (1.30). Cu (1.31) si (1.32), ecuatia (1.9) devine

U SUHX V T V H + X = C (1.50)

care, la randul ei, este echivalenta cu

SY T + Y = C (1.51)

unde Y si C sunt date de (1.34). Ecuatia (1.9) are o solutie unica X dacasi numai daca ecuatia (1.51) are o solutie unica Y. Daca yj = Y ej , tj =


27/177


T ej , cj = Cej , j = 1 : n atunci, tinand seama de structura superiortriunghiulara a matricelor T si S, din (1.51) obtinem

(St11 + Im)y1 = c1,

(Stjj + Im)yj = cj Sj1

k=1 tkjyk, j = 2 : n.(1.52)

Ecuatiile (1.52) sunt sisteme liniare superior triunghiulare care admit solutiiunice daca si numai daca matricele Stjj + Im, j = 1 : n sunt nesingulare,respectiv

siitjj + 1 = 0, i 1 : m, j 1 : n. (1.53)Cum nsa

{sii

|i = 1 : m

}= (A) si

{tjj

|j = 1 : n

}= (B) si, n

plus, rezolvarea sistemelor (1.52) n conditiile (1.53) este posibila n ordineaj = 1, 2,...,n n care termenul liber este calculabil, rezulta ca propozitia estedemonstrata.

Demonstratia de mai sus conduce imediat la algoritmul de rezolvare aecuatiei matriceale Sylvester discrete n varianta Schur-Schur.

Algoritmul 1.3 (Fiind date matricele A Rmm, B Rnn astfel ncat ij = 1 pentru toti i (A) si j (B)si matricea C Rmn algoritmul calculeaza solutia X a ecuatieiSylvester discrete AXB + X = C prin reducerea matricelor Asi B la forma Schur complexa. Algoritmul utilizeaza functia fscpentru calculul formei Schur complexe si a matricei de transfor-mare corespunzatoare.)

1. [ U, S] = fsc(A)

2. [ V, T ] = fsc(B)

3. C C = UHCV4. Pentru j = 1 : n

1. Daca j > 1 atunci

1. cj cj S (j1

k=1 tkjyk )

2. Se rezolva sistemul superior triunghiular(Stjj + Im)yj = cj .

5. X = Re(U Y VH).

Avand n vedere similitudinea algoritmilor 1.1 si 1.2, similitudine carese extinde la celelalte variante, lasam n sarcina cititorului elaborarea vari-antei HessenbergSchurca si justificarea faptului ca urmatorul algoritm, cunotatiile din algoritmul 1.2, rezolva problema n cazul matricelor de intrareA si B simple (diagonalizabile).


28/177

1.4. ECUATII LIAPUNOV 27

Algoritmul 1.4 (Date matricele simple A Rmm , B Rnn cu ij = 1 pentru toti i (A), j (B) si ma-tricea C Rmn, algoritmul calculeaza solutia X a ecuatieiAXB + X = C pe baza determinarii valorilor si vectorilor pro-prii ale matricelor A si B cu ajutorul functiei vvp.)

1. [ U, l ] = vvp(A)

2. [ V, m ] = vvp(B)

3. Se rezolva sistemul matriceal U Z = CV n raport cu Z.

4. Pentru i = 1 : m

1. Pentru j = 1 : n

1. yij

= zij

/(lim

j+ 1)

5. Se rezolva sistemul matriceal XV = U Y n raport cu X.

6. X = ReX.

1.4 Rezolvarea ecuatiilor matricealeLiapunov

Ecuatiile matriceale Liapunov continua ATX+XA = C si discreta ATXAX = C sunt cazuri particulare ale ecuatiilor Sylvester corespunzatoare si,prin urmare, rezultatele paragrafului precedent sunt aplicabile n ntregime.Concret, avand n vedere ca (AT) = (A), obtinem urmatoarele consecinte

ale propozitiilor 1.3 si 1.4.

Corolar 1.1 Ecuatia matriceala Liapunov continua (1.10) admite o solutieunica daca si numai daca matriceaA nu are nici o pereche de valori propriiopuse, i.e.

i + j = 0, i, j (A), (1.54)In particular, 0 (A) i.e. matricea A este nesingulara.

Corolar 1.2 Ecuatia matriceala Liapunov discreta (1.11) admite o solutieunica daca si numai daca matriceaA nu are nici o pereche de valori propriiinverse, i.e.

ij = 1, i, j (A), (1.55)In particular, 1 (A) .

Algoritmii 1.1, 1.2, 1.3 si 1.4 se simplifica n sensul ca este necesarao singura reducere la forma Schur, respectiv un singur calcul de valori sivectori proprii.


29/177


I. Intr-adevar, sa consideram mai ntai cazul ecuatiei Liapunov continue(1.10) si fie reducerea (1.31) a matricei A Rnn la forma Schur complexa.Rezulta

AT = AH = U SHUH. (1.56)

Cu (1.31), (1.56) ecuatia Liapunov (1.10) devine

U SHUHX+ XUSUH = C (1.57)

respectivSHY + Y S = C (1.58)

unde

Y = U

H

XU,

C = U

H

CU. (1.59)Intrucat SH are o structura inferior triunghiulara, scriind ecuatia ma-triceala (1.58) pe coloane obtinem sistemele inferior triunghiulare

(SH + sjjIn)yj = cj j1k=1

skjyk, j = 1 : n. (1.60)

In conditiile corolarului 1.1, aceste sisteme sunt nesingulare si pot fi rezol-vate n ordinea j = 1, 2 ..., n. Cu observatia ca solutia X a ecuatiei initialese calculeaza apoi din prima relatie (1.59) rezulta urmatorul algoritm.

Algoritmul 1.5 (Date matricele A Rnn, satisfacand(1.54), si C

Rnn oarecare, algoritmul calculeaza solutia X a

ecuatiei Liapunov ATX + XA = C pe baza reducerii matriceiA la forma Schur complexa cu ajutorul functiei fsc.)

1. [ U, S] = fsc(A)

2. C UHCU3. Pentru j = 1 : n


1. cj = cj j1

k=1 skjyk2. Se rezolva sistemul inferior triunghiular

(SH + sjjIn)yj = cj.

4. X = Re(U Y UH).

In multe aplicatii matricea C din (1.10) este simetrica. Intr-un astfelde caz apar simplificari suplimentare datorita faptului ca, n condit iile deexistenta si unicitate mentionate, si matricea solut ie este simetrica. Intr-adevar, daca X este solutia, se verifica imediat ca XT satisface ecuatia


30/177


matriceala (1.10) si, cum solutia este unica, rezulta X = XT. Acest faptpoate fi exploatat observand ca, datorita simetriei matricelor C si X, ma-tricele Y si C din (1.58) sunt hermitice si n consecinta este suficient sa secalculeze, de exemplu, numai triunghiul superior al matricelor C, Y si, nfinal, X (exercitiul 1.11).

II. Algoritmul de rezolvare a ecuatiei Liapunov discrete (1.11), bazatpe reducerea matricei A la forma Schur complexa, urmareste ndeaproapeschema de mai sus.

Algoritmul 1.6 (Date matricele A Rnn, satisfacand(1.55), si C

Rnn oarecare, algoritmul calculeaza solutia X a

ecuatiei Liapunov ATXA X = C pe baza reducerii matriceiA la forma Schur complexa cu ajutorul functiei fsc.)

1. [ U, S] = fsc(A)

2. C UHCU3. Pentru j = 1 : n


1. cj = cj SH(j1

k=1 skjyk)

2. Se rezolva sistemul inferior triunghiular(SHsjj In)yj = cj .

4. X = Re(U Y UH).

Tratarea cazului n care matricele C si, drept consecinta, X sunt simetriceface obiectul exercitiului 1.11.

Variantele bazate pe calculul valorilor si vectorilor proprii ale matriceiA sunt, de asemenea, propuse cititorului (exercitiul 1.12).

Cazul important, n care matricea A este stabila (vezi cap. 5) iar ma-tricea C este simetrica si negativ (semi)definita este abordat n capitolul 7(se pot consulta [ X ] si [ 6 ], [ 7 ]).

Incheiem prezentarea metodelor de rezolvare a ecuatiilor matriceale li-niare prin evidentierea posibilitatii de a reduce rezolvarea unei ecuatii Lia-punov continue la rezolvarea unei ecuatii corespondente discrete si reciproc.Pentru nceput sa observam ca daca ecuatia Liapunov discreta (1.11) ad-mite o solutie unica atunci, conform (1.55), 1 (A) respectiv matriceaA In este nesingulara si, deci, putem defini transformarea (omografica)

B = (A In)1(A + In). (1.61)


31/177


Apoi, din (1.61) rezulta ca 1 (B), i.e B In este nesingulara (exercitiul1.13), deci

A = (B + In)(B In)1, (1.62)respectiv,

AT = (BT In)1(BT + In) (1.63)expresii care, introduse n ecuatia Liapunov discreta (1.11), conduc, dupacateva prelucrari elementare care conserva solutia, la

BTX+ XB = D (1.64)

undeD =

1

2(BT In)C(B In). (1.65)

In consecinta, ecuatia Liapunov continua definita de (1.64), (1.61) si (1.65)are aceeasi solut ie cu ecuatia Liapunov discreta (1.11) si poate servi casuport pentru o modalitate de rezolvare a acesteia din urma. Reciproc,data ecuatia Liapunov continua (1.64), respectiv date matricele n n B siD cu valorile proprii ale matricei B satisfacand conditia i + j = 0 oricarear fi i, j (B) ( n particular 0 (B)) putem gasi solutia X a acesteiecuatii rezolvand ecuatia Liapunov discreta (1.11) cu A dat de (1.62) si

C = 2(BT In)1D(B In)1 (1.66)

obtinut din (1.65).

Proprietati numerice

Algoritmii prezentati utilizeaza n exclusivitate transformari ortogonale (u-nitare) ceea ce le confera remarcabile proprietati numerice. Toti algoritmiiinclud cel putin o executie a algoritmului QR de aducere iterativa a un-eia din matricele de intrare la forma Schur reala sau complexa urmata derezolvarea unor sisteme triunghiulare sau de tip Hessenberg. AlgoritmulQR este admis a fi un algoritm numeric stabil ca si algoritmii de rezolvareprin substitutie a sistemelor triunghiulare. Mai mult, se poate arata caalgoritmii din prezentul capitol sunt numeric stabili n ntregul lor, i.e.solutia calculata reprezinta solutia exacta a problemei respective cu datelede intrare perturbate nesemnificativ. In consecinta, preciziile efective ce seobtin depind de conditionarile numerice ale matricelor initiale si de nivelulpracticat al tolerantelor.


32/177


Programe MATLAB disponibile

Pentru rezolvarea ecuatiilor matriceale Sylvester si Liapunov continue estedisponibila functia lyap, care implementeaza algoritmii 1.1 si 1.5 (se ape-leaza functiile schur pentru aducerea la forma Schur reala si rsf2csfpentruobtinerea formei Schur complexe). In cazul discret se poate utiliza functiadlyap, care impementeaza metoda transformarii omografice (1.61).

Exercitii

E 1.1 Fie date matricele A

Rmm nesingulara si C

Cmn. Scrieti

algoritmii de rezolvare a sistemului matriceal AX = C folosind eliminareagaussiana cu pivotare partiala si completa. Evaluati numarul asimptoticde operatii aritmetice. Cum procedati daca matricea A este simetrica,eventual pozitiv definita ?

E 1.2 Se considera date matricele A Rpm, C Rpn si fieA = U11V1

T, 1 = diag(1, 2, . . . , r), 1 2 . . . r > 0descompunerea valorilor singulare a lui A. Aratati ca X = V1

11 U

T1 C

este matricea de norma Frobenius minima dintre toate matricele X Rmncare minimizeaza norma Frobenius a reziduului matriceal R = C AX.

E 1.3 Admitand ca dispuneti de o functie care calculeaza valorile singulare

si matricele de transformare corespunzatoare pentru o matrice de intrareA Rpm data, scrieti algoritmul de rezolvare, n sensul celor mai micipatrate, a ecuatiei matriceale AX = C.

E 1.4 Data o matrice A R22 avand (A) C\R, determinati o matricede rotatie bidimensionala complexa P C22, definita de

P =

c s

s c

, c2 + |s|2 = 1, c R

astfel ncat matricea S = PHAP C 22, PH def= PT sa fie superior tri-unghiulara.

E 1.5 Se da o matrice A R

nn n forma Schur reala. Se cere unalgoritm de calcul al unei matrice unitare 1 Q Cnn si al matricei supe-

1O matrice P Cnn se numeste unitara daca PHP = P PH = In si joaca, ncalculul matriceal complex, acelasi rol cu cel jucat de matricele ortogonale n calcululmatriceal real.


33/177


rior triunghiulare S Cnn astfel ncat QHAQ = S. (Matricea S astfelobtinuta se numeste forma Schur complexa a matricei A).

Indicatie. Matricea Q poate fi un produs de matrice unitare de forma

Qi =

I 0 00 Pi 0

0 0 I

Cnn

cu Pi C22 o rotatie complexa bidimensionala.

E 1.6 Fie A C22 superior triunghiulara. Sa se calculeze o matriceunitara U

C22 astfel ncat

UHAU =

a22

0 a11

,

i.e. sa se obtina permutarea elementelor diagonale.

E 1.7 Daca B Cnn este o matrice superior triunghiulara data, se cere unalgoritm de calcul al unei matrice unitare V Cnn si a matricei superiortriunghiulare S Cnn astfel ncat VHBV = S si |s11| |s22| . . . |snn|.

Observatie. Daca B este forma Schur complexa a unei matrice A Rnn, atunci matricea S de mai sus poarta numele de forma Schur com-plexa ordonata a matricei A. Evident, se poate utiliza orice alt criteriu de

ordonare a elementelor diagonale.Indicatie. Se foloseste o secventa de transformari de forma

Vi =

I 0 00 Ui 0

0 0 I

unde Ui sunt matrice 2 2 determinate n exercit iul 1.6.

E 1.8 Admitand ca numarul asimptotic de operatii aritmetice n formatvirgula mobila necesar pentru executia algoritmului QR de aducere a uneimatrice reale n n la forma Schur reala, (inclusiv calculul matricei detransformare) este, n medie statistica, conform [ VI ], N1 = 12n3, evaluati

a) numarul asimptotic de operatii impus de executia algoritmului de laexercitiul 1.7;

b) numarul asimptotic de operatii necesar pentru rezolvarea, cu ajutorulalgoritmului 1.2 (varianta SchurSchur), a ecuatiei Sylvester continue (1.8);


34/177


35/177


a) Daca B = (A In)1(A + In) atunci pentru orice i (A) avemi = (i + 1)/(i 1) (B);

b) {0, 1} (B) = ;c) (B) (B) = .

E 1.14 Elaborati algoritmii de tip Schur-Schur pentru rezolvarea ecuatiilorSylvester AX+ XB = C si AXB + X = C care sa utilizeze formele Schurreale ale matricelor A si B si exclusiv o aritmetica reala 2.

Indicatie. Utilizand reducerea matricelor A si B la forma Schur realaecuatiile (1.35), respectiv (1.51), au matricele S si T cvasisuperior triunghi-ulare. Partitionand matricea Y n blocuri conform dimensiunilor blocurilor

diagonale ale matricelor S si T, scriind (1.35), respectiv (1.51), pe bloc-coloane si utilizand o tehnica de bloc-substitutie napoi se obtine ca ecuatiile(1.35), respectiv (1.51), se reduc la un set de ecuatii Sylvester continue,respectiv discrete, definite de blocurile diagonale ale matricelor S si T.Aceste ecuatii se pot rezolva, n ordinea sugerata mai sus, ntr-o aritmeticareala, prin desfasurarea lor n sisteme uzuale de ordin cel mult patru (veziexercitiul 1.9).

E 1.15 Elaborati algoritmii de tip Schur pentru rezolvarea ecuatiilor Lia-punov ATX+ XA = C si ATXA X = C care sa utilizeze exclusiv formaSchur reala a matricei A si o aritmetica reala 3.

Indicatie. Adaptati indicatiile de la exercitiul precedent. Pentru detaliiputeti consulta [ 4-6] si [ X ].

Bibliografie

[ 1 ] Davison E.J., Man F.T. The Numerical Solution ofATQ+QA = C,IEEE Trans. Automat. Contr., vol. AC-13, pp. 448449, 1968.

[ 2 ] Berger C.S. A Numerical Solution of the Matrix EquationP = PT+S,IEEE Trans. Automat. Contr., vol.AC-16, pp. 381-382, 1971.

[ 3 ] Bartell R.H., Stewart G.W. Solution of the Matrix Equation AX +XB = C, Commun. Ass. Comput. Mach., vol 15, pp. 821-826, 1972.

[ 4 ] Barraud A.Y. A Numerical Algorithm to Solve ATXA X = Q, IEEE

Trans. Automat. Contr., vol.AC-22, pp. 883885, 1977.2Varianta Schur-Schur reala pentru rezolvarea ecuatiei Sylvester continue este referita

n literatura de specialitate cu denumirea de algoritmul Bartels-Stewart [ 3 ] .3Varianta Schur reala pentru rezolvarea ecuatiei Liapunov discrete este referita n

literatura de specialitate cu denumirea de algoritmul Kitagawa-Barraud[ 4 , 5 ] .


36/177


[ 5 ] Kitagawa G. An Algorithm for Solving the Matrix EquationX=FTXF+S, Int. J. Control, vol. 25, pp. 745753, 1977.

[ 6 ] Sima V. Comparison of some Algorithms for Solving the Lyapunov-typeEquations, Rev. Roum. Sci.Techn. Electrotechn. Energ., vol. 25, pp.625632, 1980.

[ 7 ] Hammarling S.J. Numerical Solution of the Stable, Non-negative Def-inite Lyapunov Equation, IMA J. Numer. Anal., vol. 2, pp. 303323,1982.


37/177



38/177

Capitolul 2

Calculul functiilor dematriceExponentiala matriceala

Calculul functiilor de matrice este necesar n multe aplicatii. In particu-lar, analiza comportarii dinamice a sistemelor liniare continue face apel laexponentiala matriceala, pentru calculul careia a fost elaborata o ntreagaserie de tehnici numerice mai mult sau mai put in performante. In con-

secinta, prezentam n continuare unele metode generale de calcul pentrufunctiile de matrice, precum si metodele specifice, cele mai apreciate, decalcul pentru functia exponentiala de argument matriceal.

2.1 Functii de matrice

Consideram o functie f : D C C definita pe o multime D din planulcomplex si fie A Cnn o matrice data. Ne propunem mai ntai sa definimnotiunea de functie de matrice adica semnificatia expresiei

F = f(A). (2.1)

Pentru nceput observam ca daca f este o functie polinomiala

f(z) =

Ni=0

cizi, ci C, i = 0 : N (2.2)

37


39/177

38 CAPITOLUL 2. FUNCTII DE MATRICE

atunci matricea

F =Ni=0

ciAi (2.3)

este bine definita si poate fi numita valoarea polinomului f n punctul (saupe matricea) A.

Fie acum A(z) polinomul minimal al matricei A1 si i C, i = 1 : l,

zerourile acestuia, zeroul i avand ordinul de multiplicitate mi. Avem

A(z) = zm +

m1k=0

izi =

li=1

(z i)mi (2.4)

unde

mdef= grad(A) =

li=1

mi. (2.5)

Conform teoremei mpartirii cu rest, exista unic doua polinoame q, r cugrad(r) < m astfel ncat

f = A q + r. (2.6)

Polinomul minimal A fiind un polinom anulator pentru matricea A, i.e.A(A) = 0, din (2.6) rezulta

Fdef= f(A) = r(A). (2.7)

In consecinta, valoarea polinomului f n punctul A este aceeasi cu cea apolinomului r si acelasi lucru se poate spune despre orice alt polinom alcarui rest la mpartirea prin A este r. Polinomul r poate fi determinatprin aplicarea algoritmului de mpartire cu rest a polinoamelor.

O alta modalitate de calcul a polinomului r se bazeaza pe faptul ca din(2.4) rezulta

(k)A (i) = 0, k = 0 : mi 1, (2.8)

si, n consecinta, cei (cel mult) m coeficienti ai polinomului r se pot deter-mina, conform (2.6), prin rezolvarea sistemului de m ecuatii liniare

r(k)(i) = f(k)(i), i = 1 : l, k = 0 : mi 1, (2.9)

unde indicele superior indica ordinul derivatei.1Amintim ca A este polinomul monic de grad minim cu proprietatea A(A) = 0.

Daca p(z) = det(zI A) este polinomul caracteristic al matricei A iar d(z) este cel maimare divizor comun (monic) al tuturor minorilor de ordinul n1 ai matricei caracteristicezI A atunci A(z) = p(z)/d(z).


40/177

2.1. FUNCTII DE MATRICE 39

Relatiile (2.7) si (2.9) servesc ca baza pentru definirea valorii n punctulA a oricarei functii f care admite derivatele cerute n (2.9). Pentru oabordare formala introducem

Definitia 2.1 Fie A Cnn si = {(i, mi) | i = 1 : l, i C, mi N} (2.10)

multimea zerourilor polinomului minimal A al matricei A mpreun a cumultiplicitatile respective. Daca functia f : D C C este analitica peo multime deschisa D ce contine punctele i , i = 1 : l atunci spunem caf este definita pe spectrul matricei A iar multimea valorilor functiei f pespectrul matricei A este

f() =

f(k)(i) | i = 1 : l, k = 0 : mi 1

. (2.11)

In particular, o functie ntreaga f (i.e. analitica pe D = C) este definitape spectrul oricarei matrice A Cnn.

In conditiile definitiei 2.1 introducem

Definitia 2.2 Fie date o matrice A Cnn si o functie f definita pespectrul lui A. Daca polinomul minimal A al matricei A are gradul matunci polinomulr de grad cel multm1, unic determinat de sistemul liniar(2.9), se numeste polinomul de interpolare Lagrange-Sylvester al functieif pe spectrul matricei A.

Putem acum sa introducem notiunea de functie de matrice extinzandrelatia (2.7) de la polinoame la orice functie definita pe spectrul matriceiargument.

Definitia 2.3 Daca f este o functie definita pe spectrul unei matrice Aatunci valoarea functiei f n punctul A este

f(A) = r(A), (2.12)

unde r este polinomul de interpolare Lagrange-Sylvester al functiei f pespectrul lui A.

In scopul evidentierii cazurilor cand o functie de matrice este reala, n

continuare vom utiliza urmatoarea terminologie

Definitia 2.4 O functie f : D C C este reala pe spectrul matriceiA Rnn daca

f(i) = f(i), i , i (A),


41/177


unde z este notatia pentru conjugatul numarului complex z.

In particular f(i) R pentru orice i (A) R.

In acest context avem

Propozitia 2.1 Daca functia f este reala pe spectrul matricei reale A Rnn atunci polinomul de interpolare Lagrange-Sylvester are coeficientiireali si, n consecint a, matricea F = f(A) este reala.

Exemplul 2.1 Fie f(z) = ez si

A = 1 1 00 2 00 0 2

.Avem A(z) = z

2 3z + 2 si = {(1, 1) , (2, 1)} . Daca r(z) = r1z + r0,atunci sistemul (2.9) se scrie

r1 + r0 = e2r1 + r0 = e

2

de unde rezulta r0 = 2e e2 , r1 = e2 e si, n consecinta,

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

e e2 e 00 e2 0

0 0 e2

.

Definitia 2.3 a functiilor de matrice pe baza polinomului de interpolareLagrange-Sylvester pune n evidenta urmatoarele aspecte.

a) Evaluarea oricarei functii de matrice se reduce la evaluarea unuipolinom matriceal. Totusi, desi sunt stabilite expresii analitice pentrucoeficientii polinomului de interpolare Lagrange-Sylvester, aceasta cale nueste recomandata pentru calculul numeric al functiilor de matrice din con-siderente de eficienta si stabilitate numerica.

b) Valoarea functiei f n punctul A este determinata exclusiv de multi-mea valorilor functiei f pe spectrul matricei A.

c) Daca matricea A nu are valori proprii multiple, atunci polinomulminimal coincide cu polinomul caracteristic iar sistemul liniar (2.9) a caruisolutie furnizeaza coeficientii polinomului de interpolare Lagrange-Sylvesterdevine

r(i) = f(i), i = 1 : n. (2.13)


42/177


Urmatoarele rezultate (ale caror demonstratii pot fi gasite, de exemplu,n [ VII ]), permit evidentierea unor proprietati utile n elaborarea unorproceduri de calcul efectiv.

Teorema 2.1 Fie (A) spectrul matricei A Cnn si D C un domeniucu frontiera suficient de neteda astfel ncat (A) D. Daca f este o

functie analitica pe D atunci

f(A) =1

2i

(zI A)1f(z)dz . (2.14)

Expresia (2.14) poate servi ca definitie pentru functiile analitice (pe un

domeniu) iar calculul integralei Cauchy

fij =1

2i

eTi (zI A)1ejf(z)dz (2.15)

poate fi efectuat cu ajutorul teoremei reziduurilor.

Exemplul 2.2 Consideram functia f si matricea A din exemplul 2.1.Avem (A) = {1, 2, 2} si

(zI A)1 =

1z1

1(z1)(z2) 0

0 1z2 0

0 0 1z2

.

Prin urmare, f fiind analitica n tot planul complex, putem alege un dome-niu simplu conex oarecare ce satisface (A) D. Obtinem

f11 =1

2i

ez

z 1 dz = Rez

ez

z 1

z=1

= e,

f22 = f33 =1

2i

ez

z 2 dz = Rez

ez

z 2

z=2

= e2,

f13 =1

2i

ez

(z 1)(z 2) dz = Rez ez

(z 1)(z 2)z=1 ++Rez

ez

(z 1)(z 2)

z=2

= e + e2,


43/177


celelalte elemente ale matricei F = eA fiind, evident, nule. In cele ce urmeaza vom aborda n exclusivitate cazul functiilor analitice

si vom putea utiliza, n consecinta, relatia (2.14) ca o relatie definitorie afunctiilor de matrice.

Mentionam, de asemenea, posibilitatea exprimarii functiilor de matriceprin serii matriceale de puteri. In acest sens avem urmatoarea

Teorema 2.2 Daca functia f(z) se poate dezvolta n serie de puteri n jurul punctuluiz = z0

f(z) =

k=0k(z z0)k (2.16)

si seria este convergenta n discul |z z0| < r atunci aceasta dezvoltareramane valabila daca argumentul scalar este nlocuit cu argumentul ma-triceal A

f(A) =k=0

k(A z0I)k (2.17)

oricare ar fi matricea A al carei spectru se afla n interiorul discului deconvergent a.

Din aceasta teorema rezulta, printre altele, ca dezvoltarile

eA =

k=0

1

k!Ak,

sin A =k=0

(1)k(2k + 1)!

A2k+1, cos A =k=0

(1)k(2k)!

A2k

sunt valabile oricare ar fi matricea A Cnn.Exemplul 2.3 Reluam matricea A si funct ia f din exemplul 2.1. Prininductie rezulta imediat

Ak =

1 2k 1 00 2k 0

0 0 2k

si deci

F = eA =k=0

1

k!Ak =

k=0 1k! k=0 2

k1k! 0

0

k=02k

k! 0

0 0

k=02k

k!

=


44/177


=

e e2 e 00 e2 0

0 0 e2

.

In continuare, prezentam cateva proprietati ale functiilor de matrice,utile n dezvoltarile procedurale care fac obiectul metodelor de calcul reco-mandate ca fiind cele mai bune n momentul actual.

Propozitia 2.2 Daca functiaf este definita pe spectrul matricei A si

B = T AT1 (2.18)

unde T este o matrice nesingulara, atunci

f(B) = T f(A)T1. (2.19)

Demonstratia este imediata pentru functiile analitice pe baza relatiei defi-nitorii (2.14) ntrucat din (2.18) rezulta

(zI B)1 = (T(zI A)T1)1 = T(zI A)1T1.

Transferul transformarilor de asemanare (2.18) de la nivelul argumentelormatriceale la nivelul functiilor are o importanta decisiva n elaborarea unortehnici de calcul adecvate care se pot focaliza asupra unor structuri ma-triceale remarcabile.

Propozitia 2.3 Daca matricea A este (bloc-) diagonala

A = diag(A11, A22, . . . , App) (2.20)

atunci f(A) este (bloc-) diagonala si

Fdef= f(A) = diag (f(A11), f(A22), . . . , f (App)) . (2.21)

Demonstratie. Din faptul ca

(zI A)1 = diag

(zI A11)1, (zI A22)1, . . . , (zI App)1

(2.21) rezulta direct din (2.14). Observatia 2.1 Conform propozitiilor 2.2 si 2.3, daca

T1AT = diag(J1, . . . , J p)


45/177


este forma Jordan a lui A, atunci avem

f(A) = T diag(f(J1), . . . , f (Jp))T1,

si calculul lui f(A) se reduce la calculul matricelor f(Jk), unde Jk suntblocurile Jordan. Considerand un bloc Jordan generic

J =

1 0 00 1 0

. . .. . .

. . . 10 0 0

si utilizand definitia 2.3 obtinem

f(J) =

f() f() . . .f(m1)()(m 1)!

0 f() . . .f(m2)()(m 2)!

......

. . ....

0 0 . . . f ()

,

(exercitiul 2.4).

Propozitia 2.4 Daca matriceaA

Cnn este superior (inferior) triunghi-

ulara, iar f este o functie definita pe spectrul lui A atuncia) F = f(A) este superior (inferior) triunghiulara;

b) fii = f(aii), i = 1 : n.

Demonstratie.

a) Se stie ca inversa unei matrice superior (inferior) triunghiulare nesin-gulare este superior (inferior) triunghiulara. In consecinta (zI A)1 estesuperior (inferior) triunghiulara si, din (2.14), rezulta ca aceeasi structurao are si matricea F = f(A).

b) Daca A este o matrice triunghiulara atunci elementele diagonale alematricei (zI A)1 sunt

1z aii , i = 1 : n.

Functia f fiind definita pe spectrul matricei A, rezulta ca i = aii, i =1 : n nu sunt puncte singulare pentru f si, n consecinta, utilizand formula


46/177


integrala a lui Cauchy pentru elementele diagonale din (2.14), rezulta

fii =1

2i

f(z)

z aii dz = f(aii), i = 1 : n.

q.e.d.

Propozitia 2.5 Fie (A) = {1, . . . , n} spectrul matricei A Cnn.Atunci pentru orice functie f definita pe spectrul lui A avem

(f(A)) = {f(1), . . . , f (n)} . (2.22)

Demonstratie. Fie B = QH

AQ forma Schur complexa a matricei A. AtunciB este superior triunghiulara si bii = i, i = 1 : n. Aplicand propozitiile2.2 si 2.4 rezulta

f(A) = Qf(B)QH

si, deci, valorile proprii ale matricei f(A) sunt f(bii) = f(i), i = 1 : n.

Propozitia 2.6 Matricele A Cnn si F = f(A) comuta i.e.

Af(A) = f(A)A. (2.23)

Demonstratie. Existenta matricei F = f(A) presupune ca functia f estedefinita pe spectrul matricei A. In acest caz nsa si functia h definita de

h(z) = zf(z) = f(z)z este definita pe spectrul matricei A. Din faptul cah(A) are semnificatie rezulta (2.23). Teorema 2.1 mpreuna cu proprietatea (2.23) stabileste un izomorfism

f(z) ; f(A) ntre algebra comutativa a functiilor analitice f(z) pe spectrullui A si algebra (de asemenea comutativa) a functiilor de matrice f(A).

Observatia 2.2 Din propozitiile 2.2 si 2.4 rezulta ca daca QHAQ =B este forma Schur a lui A atunci f(A) = Qf(B)QH, unde G = f(B)este superior triunghiulara cu gii = f(bii). Tinand seama ca forma Schurprezinta proprietati numerice net superioare formei canonice Jordan si capropozitia 2.6 sta la baza unui algoritm fiabil de calcul al functiilor dematrice triunghiulare (vezi sectiunea urmatoare) deducem ca rezultatele demai sus fundamenteaza o categorie importanta de proceduri de calcul al

functiilor de argument matriceal. Proprietatile urmatoare, care ilustreaza izomorfismul de algebre comu-

tative amintit anterior, permit extinderea unor identitati care pun n relatiefunctii de argument scalar cu functii de matrice.


47/177


Propozitia 2.7 Fie g(z1, z2, . . . , zp) un polinom n nedeterminatele z1,z2, . . . , zp si fi, i = 1 : p functii definite pe spectrul matricei A Cnn.Definim functia

h(z) g(f1(z), . . . , f p(z)). (2.24)Atuncih este o functie definita pe spectrul matricei A si daca

h(k)(i) = 0, i = 1 : l, k = 0 : mi 1 (2.25)

unde (i, mi), i = 1 : l sunt zerourile polinomului minimal A cu multi-plicitatile respective, atunci

H

def

= h(A) = g(f1(A), . . . , f p(A)) = 0. (2.26)

Demonstratie. Daca ri, i = 1 : p sunt polinoamele de interpolare Lagrange-Sylvester ale functiilor fi pe spectrul matricei A, atunci consideram functiau definita de

u(z) = g(r1(z), . . . , rp(z))

care este, evident, o functie polinomiala de argument scalar. In acesteconditii, din (2.25), (2.26) si modul de definire a polinoamelor de interpolareLagrange-Sylvester (2.9) rezulta

u(k)(i) = 0, i = 1 : l, k = 0 : mi 1,

respectiv polinomul de interpolare al functiei u este identic nul, de undeobtinem

u(A) = g(r1(A), . . . , rp(A)) = g(f1(A), . . . , f p(A)) = 0

q.e.d. Cateva consecinte imediate ale propozitiei 2.7, care evidentiaza trans-

ferul unor formule scalare n cazul matriceal sunt urmatoarele relatii utile.

a) Fie polinomul g(z1, z2) = z21 +z

221 si funct iile f1(z) = sin z, f2(z) =

cos z. Atunci functia definita n (2.25)

h(z) g(sin z, cos z) = sin2z + cos2z 1

este identic nula si prin urmare (2.26) este satisfacuta pentru orice matriceA. De aici rezulta ca, pentru toate matricele patrate A, avem

sin2A + cos2A = I.


48/177


b) Similar, considerand polinomul g(z1, z2) = z1z21 si functiile f1(z) =ez, f2(z) = e

z obtinem eAeA = I, adica

eA =

eA1

.

c) De asemenea, cu g(z1, z2, z3) = z1 + z2 + z3 si f1(z) = eiz, f2(z) =

cos z, f3(z) = i sin z rezulta ca pentru orice matrice patrata A esteadevarata formula

eiA = cos A + i sin A.

Propozitia 2.8 Daca

(z) =g(z)

h(z)(2.27)

este o functie rationala ireductibila, definita pe un domeniu care nu contineradacinile polinomului h(z) iar A Cnn atunci functia este definita pespectrul matricei A daca si numai daca

h(i) = 0 i (A). (2.28)Daca (2.28) sunt satisfacute atunci

(A) = g(A)[h(A)]1 = [h(A)]1g(A). (2.29)

Demonstratie. Functia este indefinit derivabila n toate punctele n careh(z) nu se anuleaza. In consecinta (k)(i) sunt definite daca si numai

daca (2.28) sunt adevarate. Pe de alta parte din (2.28) si propozitia 2.5rezulta ca 0 (h(A)) i.e. matricea h(A) este nesingulara. Aplicand acumpropozitia 2.7, din identitatea

h(z)(z) = (z)h(z) = g(z)

rezultah(A)(A) = (A)h(A) = g(A)

care mpreuna cu nesingularitatea matricei h(A) conduc la (2.29).

Propozitia 2.9 Daca functia compusa h = g f este definita pe spectrulmatricei A Cnn atunci

h(A) = g(f(A) (2.30)

i.e. h(A) = g(B) unde B = f(A).

Demonstratia este propusa cititorului.


49/177


2.2 Calculul functiilor de matrice

Tehnicile numerice de evaluare a functiilor de matrice, recomandate deexperienta numerica acumulata, evita determinarea efectiva a polinomuluide interpolare Lagrange-Sylvester, n primul rand din motive de eficienta[ 1 ]. Metodele care s-au impus si pentru care exista si unele rezultateteoretice privind problemele legate de stabilitatea numerica se pot mpartin doua categorii:

a) metode bazate pe calculul valorilor proprii;

b) metode aproximative bazate pe trunchierea unor dezvoltari n serie.

Prezentam, n continuare, procedurile bazate pe calculul valorilor pro-

prii; metodele aproximative, bazate pe trunchierea unor dezvoltari n serie,vor fi expuse n paragraful urmator, dedicat cazului particular dar impor-tant pentru teoria sistemelor liniare, al calculului exponentialei matriceale.

Pentru matricele simple (i.e. diagonalizabile) calculul functiilor de ma-trice prin evaluarea valorilor si vectorilor proprii se bazeaza pe propozitiile2.2 si 2.3. Intr-adevar, n acest caz, daca V = [ v1 v2 . . . vn ] este matriceavectorilor proprii ai matricei date A Cnn atunci A = VV1 unde = diag(1, 2, . . . , n).

Prin urmare, pentru orice functie definita pe spectrul matricei A, apli-cand (2.19) si (2.21), rezulta

F = f(A) = V f()V1 = V diag(f(1), f(2), . . . , f (n))V1. (2.31)

Introducand acum functia vvp de calcul a valorilor si vectorilor propriiai unei matrice (desigur, prin aplicarea algoritmului QR) cu sintaxa

[ V, p ] = vvp(A)

unde p este vectorul valorilor proprii pentru matricea A, relatia (2.31) sta labaza urmatorului algoritm de calcul al functiilor de matrice diagonalizabile.

Algoritmul 2.1 (Date matricea simpla A Cnn si functiaf : D C C definita pe spectrul matricei A, algoritmulcalculeaza F = f(A) prin determinarea valorilor si vectorilorproprii.)

1. [ V, p ] = vvp(A)2. Pentru i = 1 : n

1. i = f(pi)

3. D = Vdiag(1, 2, . . . , n)


50/177

2.2. CALCULUL FUNCTIILOR DE MATRICE 49

4. Se rezolva sistemul matriceal liniar nesingular F V = D nraport cu F.

Evident, efortul de calcul principal este destinat calculului valorilor sivectorilor proprii si, ntr-o oarecare masura, rezolvarii sistemului matricealliniar.

In cazul functiilor f reale pe spectrul matricei reale A Rnn (vezidefinitia 2.4), conform propozitiei 2.1 matricea rezultat F este reala. Intru-cat matricele V, D si vectorul p din algoritmul 2.1 nu sunt n mod necesarreale, calculul (aproximativ) n aritmetica complexa conduce la aparitiaunor componente parazite imaginare n solutia calculata . De aceea, ntr-

un astfel de caz, algoritmul 2.1 se completeaza cu o instructiune de eliminarea acestor componente:

5. F = Re(F).

Daca matricea A nu este diagonalizabila, o extindere a algoritmului 2.1apeleaza la calculul formei Jordan a matricei A, iar aplicarea propozitiilor2.1 si 2.2 reduce calculul functiei de matrice la calculul functiilor avand caargumente blocurile Jordan. Aceasta tehnica nu este nsa recomandata nprimul rand datorita unei sensibilitati ridicate a structurii Jordan n raportcu perturbatiile ce apar la nivelul datelor de intrare.

Proprietati numerice mult mai bune, n toate cazurile, se obtin dacan locul formei canonice Jordan se utilizeaza forma Schur complexa (sau

reala) a matricei A. Aceasta alternativa este conditionata de existentaunui algoritm eficient de calcul al functiilor de matrice triunghiulare (saucvasitriunghiulare). Un astfel de algoritm, propus de B.N. Parlett [ 1 ],are la baza propozitia 2.3 si proprietatea de comutativitate formulata npropozitia 2.6. Vom deduce algoritmul pentru cazul generic al matricelorcu valori proprii distincte.

Fie T Cnn o matrice superior triunghiulara cu tii = tjj pentru oricei = j si f o functie definita pe spectrul lui T. Daca F = f(T) atunci,conform (2.23), avem

F T = T F. (2.32)

Tinand cont, conform propozitiei 2.4, ca matricea F este, de asemenea,

superior triunghiulara si scriind (2.32) pe elemente obtinem

jk=i

tikfkj =

jk=i

fiktkj , j = 1 : n, i = 1 : j, (2.33)


51/177


de unde, n ipotezele mentionate, avem

fij =1

tjj tii [tij(fjj fii) +j1

k=i+1

(tikfkj fiktkj)], j = 1 : n, i = 1 : j.

(2.34)

Relatia (2.34) este utila numai n masura n care se poate gasi o ordine decalcul a elementelor matricei F astfel ncat, la fiecare moment al procesuluide calcul, n membrul drept al expresiei (2.34) sa apara numai elemente dejacalculate. O astfel de ordine exista si ea poate fi evidentiata observand caelementele diagonale sunt calculabile, conform (2.22), cu formula

fii = f(tii) (2.35)

si ca n membrul drept al relatiei (2.34) apar elementele F(i, i : j 1)din stanga elementului fij si F(i + 1 : j, j) de sub elementul fij (vezidiagrama din figura 2.1).

j

i F(i, i : j 1) fij

' F(i + 1 : j, j)

dd

dd

dd

dd

dd

dd

Figura 2.1: Matricea F = f(T).

De exemplu, se poate adopta o ordine diagonala de calcul al elementelortriunghiului superior, dupa cum este indicat n diagrama din figura 2.2

pentru n = 5, numarul nscris n matrice marcand numarul de ordine pentrucalculul elementului de pe pozitia respectiva.

Aceasta nu este singura ordine posibila. Intr-adevar, daca efectuamcalculele pe coloane n ordinea j = 1, 2, . . . , n, pe fiecare coloana calculele


52/177

2.2. CALCULUL FUNCTIILOR DE MATRICE 51

1 6 10 13 152 7 11 14

3 8 124 9

5

Figura 2.2: Ordinea diagonala de calcul a elementelor matricei F = f(T).

efectuandu-se de jos n sus (vezi diagrama a) din figura 2.3) sau pe liniin ordinea i = n, n 1, . . . , 1, pe fiecare linie ordinea de calcul fiind de lastanga la dreapta (diagrama b) din figura 2.3), atunci n momentul calculu-lui unui element curent oarecare, elementele din stanga si de sub elementulcurent sunt deja calculate.

1 3 6 10 152 5 9 14

4 8 137 12

11

11 12 13 14 157 8 9 10

4 5 62 3

1

(a) (b)

Figura 2.3: Ordinea pe coloane (a) si ordinea pe linii (b) de calcul aelementelor matricei F = f(T).

Prezentam n continuare algoritmul corespunzator ordinii diagonale decalcul, celelalte variante facand obiectul unor exercit ii. Pentru a urmarimai usor indexarea, observam ca, atribuind indicele q directiilor paralelecu diagonala principala a matricei, rezulta urmatoarea schema de calcul.

1. Se calculeaza fii, i = 1 : n cu relatia (2.35).

2. Pentru q = 2 : n

1. Pentru i = 1 : n q + 11. Se calculeaza fi,i+q1 cu relatia (2.34).


53/177


Pentru a nu modifica indexarea folosita n expresia (2.34) a elementuluifij facem urmatoarele schimbari de indici:

p = q 1, j = i +p,cu care schema de mai sus ne conduce la urm atorul algoritm de calcul alfunctiilor de matrice triunghiulare cu elementele diagonale distincte.

Algoritmul 2.2 (Parlett)(Date o matrice superior triun-ghiulara T Cnn cu tii = tjj pentru orice i = j si o functief : D C C definita pe spectrul matricei T, algoritmul cal-culeaza, n ordine diagonala, elementele matricei superior tri-

unghiulare F = f(T).)1. Pentru i = 1 : n

1. fii = f(tii)

2. Pentru p = 1 : n 11. Pentru i = 1 : n p

1. j = i +p

2. s = tij(fjj fii)3. Daca p > 1 atunci

1. Pentru k = i + 1 : j 11. s = s + tikfkj fiktkj

4. fij =s

tjj tiiSchema recurenta care sta la baza algoritmului 2.2 este deosebit de

eficienta, n afara celor n evaluari de functii scalare de la instructiunea 1,

algoritmul necesitand un numar asimptotic de Nop 2n33 operatii n virgulamobila.

Cu acest algoritm, completat cu procedura de aducere a unei matricedate la forma Schur complexa (superior triunghiulara) prin transformariortogonale de asemanare, se obtine o procedura cu bune calitati numericepentru calculul functiilor de matrice. Introducem functia fsc de calcul aformei Schur complexe avand sintaxa

[ S, U ] = fsc(A)

unde datele de iesire sunt matricea unitara de transformare U si formaSchur complexa S a matricei A legate prin relatiile

S = UHAU, A = U SUH, (2.36)


54/177

2.3. CALCULUL EXPONENTIALEI MATRICEALE 53

indicele superior H semnificand dubla operatie de transpunere si conjugare.Notam cu funmt (functie de matrice triunghiulara) functia de calculrealizata de algoritmul 2.2 de mai sus, cu sintaxa (neformala)

F = funmt(f, T)

avand ca parametrii de intrare functia f si matricea superior triunghiu-lara T, iar ca parametru de iesire matricea F = f(T). Cu aceste notatiiprezentam

Algoritmul 2.3 (Date o matrice A cu valori proprii dis-tincte si o functie f definita pe spectrul matricei A, algoritmul

calculeaza matricea F = f(A) prin metoda aducerii matricei Ala forma Schur complexa.)

1. [ S, U] = fsc(A)

2. T = funmt(f, S)

3. F = U T UH

In cazul n care matricea A are valori proprii multiple, algoritmul 2.3nu mai este functional pentru ca formula (2.34) nu mai este aplicabila.De asemenea, n situatia existentei unor valori proprii foarte apropiate,aplicarea formulei (2.34) conduce la fenomene de instabilitate numerica.Solutia problemei ntr-un astfel de caz poate consta n utilizarea unei vari-ante bloc a algoritmului Parlett dupa o prealabila grupare a valorilorproprii apropiate n cadrul asa numitei forme Schur complexe ordonate.Pentru detalii recomandam [1], [VII].

2.3 Calculul exponentialei matriceale

Desi metodele de calcul prezentate n sectiunea precedenta sunt aplicabilepentru toate functiile definite pe spectrul matricei argument, pentru anu-mite functii, de un interes aplicativ deosebit, au fost dezvoltate procedurialternative, cu calitati numerice superioare.

In acest paragraf prezentam principalele metode pentru calculul expo-nentialei matriceale

(t) = etA, (2.37)

unde t > 0 este un parametru scalar. Dupa cum se stie, funct ia (t) estematricea de tranzit ie a starilor sistemului liniar x = Ax, adica satisfaceecuatia diferentiala matriceala liniara = A cu conditia initiala (0) = I.


55/177


Avand n vedere importanta unui calcul fiabil al exponentialei ma-triceale pentru studiul sistemelor liniare, o discutie prealabila a conditionariinumerice a problemei calculului lui (t) = etA, adica a sensibilitatii lui (t)n raport cu variatiile lui A, este necesara, ntrucat erorile de calcul ale lui(t) se traduc, prin mecanismul de propagare inversa, n perturbatii lanivelul datelor de intrare, adica al matricei A.

Notand cu (t) exponentiala asociata matricei peturbate A + E avem

= (A + E) cu (0) = I. Punand = + , obtinem

+ = A + A + E + E ,

unde (0) = I, deci (0) = 0. Pentru simplitate, neglijam termenul

E si, prin integrare, deducem

(t) =

t0

e(t)AEeA d = (t)

t0

eAEeA d .

Utilizand norma matriceala = 2

putem scrie

(t) (t) t0

eAEeA d

relatie din care rezulta

(t)

(t)

cond((t))E

A

, (2.38)

unde

cond((t))def= max

E = 1t0

eAEeA d A (2.39)

este numarul de conditie al lui A relativ la calculul lui (t) = etA.

Obtinem o evaluare grosiera a lui cond((t)) observand ca

to

eAEeA d t0

eA E eA d

unde

eA

eA , > 0.

De asemenea, avem t0

e2A dt =1

2A (e2tA 1) ,


56/177

2.3. CALCULUL EXPONENTIALEI MATRICEALE 55

deci

cond(t) 12

(e2tA 1). (2.40)

In concluzie, pentru a limita sensibilitatea lui (t) n raport cu variatiilelui A e necesar sa limitam tA, de exemplu astfel ncat

tA 1, (2.41)

fie alegand t suficient de mic, fie (n cazul n care t este impus) utilizandproprietatea

etA = (et2m

A)2m

, m 1. (2.42)

Intr-adevar, oricare ar fi tA exista m astfel ncat t2m A 1 iar dacae

t2m

A e cunoscut atunci, conform (2.42), etA poate fi calculat printr-unproces de ridicare succesiva la patrat. O evaluare a lui m rezulta observandca daca tA 1 atunci trebuie sa avem log2(tA) < m, deci

m = 1 + [log2(tA) ], (2.43)

unde [] este partea ntreaga a argumentului.In practica m se poate alege printr-un procedeu de njumatatire suc-

cesiva. Schema de calcul, numita de njumatatire si ridicare la patrat(scalling and squaring), pe scurt (1/2)2, este urmatoarea:

Schema (1/2)2

1. Se calculeaza A.2. Se initializeaza m = 0.

3. Cat timp tA 11. t t/22. m m + 1

4. Se calculeaza F = etA.

5. Cat timp m 11. F F22. m

m

1

Observatia 2.3 Avand n vedere necesitatea limitarii normei lui tA lacalculul exponentialei, schema (1/2)2 se aplica ntotdeauna indiferent demetoda utilizata pentru calculul efectiv (Taylor, Pade, etc.). Trebuie nsa


57/177


mentionat ca ridicarea la patrat (adica nmultirea a doua matrice) intro-duce erori de rotunjire suplimentare, mai precis calculul lui F2 se face cuerori mici fata de F2, dar nu neaparat fata de F2. Din aceasta cauzautilizarea schemei (1/2)2 ridica probleme de precizie daca norma etA emica 2. Pe de alta parte, daca exponentiala e utilizata pentru simulareaevolutiei unui sistem liniar, o conditie de tipul tA < 1 orienteaza alegereapasului de discretizare si, n consecinta, poate fi asigurata de la nceput.

In ipoteza ca tA e limitata (de exemplu cu schema (1/2)2), metodeleuzuale de calcul al exponentialei matriceale constau n constructia uneiaproximari locale n jurul lui z = 0 a functiei ez 3. De regula aproximarile

utilizate sunt de tip polinomial (Taylor) sau rational (Pade) si permit cal-culul lui etA prin evaluarea aproximarii considerate pentru z = tA.

I. Aproximatia polinomiala (Taylor) de ordin p este de forma

Tp(z) =

pk=0

ckzk (2.44)

2Acest lucru e posibil chiar daca tA 1 deoarece marginea etA etA utilizatamai sus este n general supraestimata (de exemplu A = 5, t = 1). De aceea este util sadispunem pentru evaluarea lui cond((t)) de estimari mai precise pentru etA, eventualcapabile sa tina seama de repartitia spectrului lui A si n special de dispersia acestuia(stiffness). Se stie ca daca A are valorile proprii i, atunci etA are valorile proprii eit

(propozitia 2.5); pe de alta parte pentru orice matrice A si orice norma matriceala

avem A (A), unde (A) def= max | i |, deci etA max | eit |. Se poate arataca daca Re i < atunci exista M astfel ncat e

tA M et , t 0. In mod analog,daca < Re i atunci exista N astfel ncat etA N et , t 0. Prin urmare dacaspectrul lui A e cuprins n banda < Re i < atunci

t0

etAEetA dt M NE

t0

e() d = M

Carte_metode de Calcul Numeric in Automatic A

Documents

Transcript of Carte_metode de Calcul Numeric in Automatic A