LIBRARY OF THE
UNIVERSITY OF ILLINOIS
AT URBANA-CHAMPAIGN
510.84
T£63c *o. ai-30
is
AUG. 5 1.976
ihe person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below.
Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University.
UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN
reftl
f
R
n
JUL
12
DEC * & 8BT&
1980
L161 — O-1096
510.814
I^63c
no.2U
Engin.
1
00
GINGERING LfBRARy y OF ILLINOIS
MS
BIBLIOGRAPHIC DATA SHLET
1. K, port No.
[JIUC-CAC-.DN-72-2U
3. Recipient's Accession Ni
•J. 1 it le ami Suht it Ic
A NUMERICAL SOLUTION OF THE MATRIX RICCATI EQUATIONS
5- Report Date
January] 20, 1972
7. Autlu>r(s)
Killion Noh
8. Performing Organ i /.at ion Kept.
No
CAC-2U
9. Performing Organization Name and Address
Center for Advanced Computation University of Illinois at Urb ana- Champaign Urbana, Illinois 61801
10. Project/Task^ Work Unit No.
11. Contract .'Grant No.
DAHC0i+ 72-C-OOOl
12. Sponsoring Organization Name and Address
U.S. Army Research Office-Durham
Duke Station
Durham, North Carolina
13. Type of Report & Period Covered
Research
14.
15. Supplementary Notes
16. Abstracts
The eigenvector solution of the time-invarient matrix Riccati equation is discussed.- The coefficient matrix of the canonical equation is allowed to have multiple eigenvalues, namely, the matrix could be either derogatory or defective. The solution of matrix Riccati equation is then calculated from a part of similar- ity transformation which should reduce the coefficient matrix to the Jordan canonical form.
|
1 |
||
|
17. Key Words and Document Analysis. 17a. Descriptors Ordinary and Partial Differential Equations ; Matrix Algebra Nonlinear and Functional Equations / 17b. ldcntif iers/Open-Ended Terms 17c. COS ATI Field/Group |
||
|
18. Availability Statement .-. . . , . ,...,,. No restriction on distribution. Available from National Technical Information Service, Springfield, Virginia 22151 |
19. Security Class (This Report ) UNri.ASSIFM-n |
21. No. of Pages 38 |
|
20. Security Class ( 1 his Page UNri.ASSIFIF.D |
22. Price - ______ |
FOI'M NTlS-iS INCV. J-//I
THIS FORM MAY 1)1- KFPRODUCHP
USCOMM-OC I4832-P72
CAC Document No. 2k UIUCDCS-R-72-U98
A NUMERICAL SOLUTION OF THE MATRIX RICCATI EQUATIONS
By Killion Noh
Center for Advanced Computation University of Illinois at Urb ana-Champaign Urbana, Illinois 6l801
January 20, 1972
Submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science in the Graduate College of the University of Illinois at Urbana-Champaign and supported in part by the Advanced Research Projects Agency of the Department of Defense and "was monitored by the U.S. Army Research Office-Durham under Contract No. DAHCO^ 72-C-OOOl.
Digitized by the Internet Archive
in 2012 with funding from
University of Illinois Urbana-Champaign
http://archive.org/details/numericalsolutioOOnohk
&NGINEEJUNG U8RAA*
ii
ABSTRACT
The eigenvector solution of the time-invariant matrix Riccati equation is discussed. The coefficient matrix of the canonical equation is allowed to have multiple eigenvalues, namely, the matrix could be either derog- atory or defective. The solution of matrix Riccati equation is then cal- culated from a part of similarity transformation which should reduce the coefficient matrix to the Jordan canonical form.
Ill
ACKNOWLEDGEMEIW
The author would like to express his thanks to his advisor. Professor Daniel L. Slotnick, and to the ILLIAC IV Project and the Center for Advanced Computation under whose employ this research was conducted. A note of appreciation is also extended to the Department of Computer Science for its cooperation and support during the time of the author's studies at the University of Illinois.
Deepest gratitude also goes to Professor Ahmed H. Sameh for suggesting this thesis topic and for his guidance throughout its preparation. Without his encouragement and patience, this work would have never been done.
IV
TABLE OF CONTENTS
Page
1. INTRODUCTION 1
2. FORMULATION OF THE PROBLEM 2
3. SOLUTION OF THE MATRIX RICCATI EQUATION 6
k. COMPUTATIONAL METHODS AND NUMERICAL RESULTS 13
LIST OF REFERENCES 18
APPENDIX I. ALGOL PROGRAMS 20
APPENDIX II. ILLIAC IV PROGRAMS 28
1 . INTRODUCTION
It is well known [1,2] that the linear regulator problem with quadratic cost functional is reduced to the problem of the matrix Riccati equation.
Although the solution of the matrix Riccati equation exists and is unique [2,3], it is not easy to obtain a numerical solution for the high-order system.
One of the numerical methods is the eigenvector solution which has been studied by several authors [U , 5 »6 ] . In their studies, however, the coeffi- cient matrix of the canonical equation is assumed to have distinct eigen- values, which is a restriction to the method of eigenvector solutions. Recently, Martensson [7] discussed the case in which the coefficient matrix is nondiagonable. Similarly, in this paper, the assumption of distinct eigenvalues will also be removed so that the coefficient matrix is allowed to have the multiple eigenvalues, namely, the matrix is either derogatory or defective. To do this, we explicitly construct the similarity transfor- mation which should reduce the coefficient matrix to the Jordan canonical form, either diagonal or nondiagonal, and then the solution of the matrix Riccati equation is calculated from a part of the similarity transformation matrix.
Since matrix Riccati equations of fairly large sizes appear in many engineering applications and since the solution algorithm can be easily constructed such that it becomes suitable for a parallel computer [8], we provide in Appendix II the basic codes written in ILLIAC IV language, GLYPNIR.
2. FORMULATION OF THE PROBLEM
In this section, the source of the matrix Riccati equation and its relation to optimal control problems are summarized in short for our further discussions [1,2,3].
Let us consider the linear time- invariant dynamical system
x(t) = Ax(t) + Bu(t)
y(t) = Cx(t) (2.1)
x(tQ)=x0 where A, B, and C are n x n, n x r, and m x n constant matrices, respec- tively, and x(t) is an n-dimensional state vector, u(t) is an r-dimensional control vector, and y(t) is an m-dimensional output vector, respectively. Further, we assume that the system (2.1) is completely observable. The optimal output regulator problem is then to determine the control u(t) which minimizes the quadratic cost functional,
V = |yT (t ) IV(t ) + \ Jf (yT(t) Qy(t) + uT(t)Ru(t))dt
t0 (2.2)
where P and Q are m x m symmetric positive semidefinite matrices and R is an r x r symmetric positive definite matrix. Substituting y(t) = Cx(t) into (2.2), V is rewritten as
V = \ xT(tf ) CTPCx(tf) + | / (xT(t)CTQC x(t) + uT(t)Ru(t))dt
t0 (2.3)
T T Since the system (2.1) is completely observable, C PC and C QC are symmetric
and positive semidefinite [l].
The optimal feedback control, therefore, is given by
u*(t) = -R_1BTK(t)x(t) - (2.U)
where K(t) is an n x n symmetric positive definite matrix which is the
solution of the matrix Riccati equation
K(t) + K(t) A + ATK(t) - K(t)BR_1BTK(t) + CTQC = 0 (2.5) with the boundary condition
K(t ) = CTPC ' (2.6)
Furthermore, the optimal trajectory is the solution of the system of differential equations
x(t) = (A - BR"1BTK(t)) x(t) (2.7)
x(tQ) = xQ
The minimum cost is given by
V • = |xT(t) K(t) x(t) (2.8)
In addition to the assumption of complete observability we further assume that the system is completely controllable. Setting P = 0 and t = °°, the output regulator problem turns out to be the following: minimize
V = \ /~(yT(t)Qy(t) + uT(t)Ru(t)) dt (2.9)
subject to
x(t) = Ax(t) + Bu(t)
y(t) = Cx(t) (2.10)
x(0) = xQ Then the optimal control is given by
u(t) = -R_1BTK x(t) (2.11)
where K is the n x n symmetric positive definite matrix which is the solu- tion of algebraic Riccati equation
KA + ATK - KBR_1BTK + CTQC = 0 (2.12)
The optimal trajectory is the solution of the system
x(t) = Gx(t)
-1 T G = A - BR B K
(2.13)
x
(0) =
Note that for an optimal control system the matrix G is diagonable and has' eigenvalues with negative real parts. Therefore, the symmetric positive definite solution K of (2.12) must be such' that G has eigenvalues with negative real parts.
Kalman [2] has shown that if the system (2.1) is completely observable and completely controllable, then
lim K(t ) = K (2.1U)
t -*»
f
We shall observe this equation in the following.
Let K(t) = Z(t)X-1(t) (2.15)
where X(t) and Z(t) are n x n matrices and X(t) is non-singular. Substituting
K(t) = (Z(t) - K(t) X(t)) X_1(t) (2.16)
into (2.5) and setting
X(t) = AX(t) - BR-1BT Z(t) (2.17)
we obtain
Z(t) = -CTQCX(t) - ATZ(t) (2.18)
Therefore (2.5) becomes a system of 2n-dimensional linear homogeneous differential equations
x(t)
"i(t) 1= \\
z(t) [-c1
-1 T -BR B
-C QC -A
Z(t)
(2.19)
Since we are interested in increasing the time t, let
t = tf - t
(2.20)
Then
where
|
X(r) |
= ¥ |
"*(*) |
||
|
_Z(t)_ |
Z(t) |
|||
|
-A |
-it" BR B |
|||
|
¥ = |
T C QC |
AT |
||
|
„ |
—* |
(2.21)
with
X(0) = I
Z(0) = CTQC (2.22)
Therefore, the solution will be given by the following matrix exponential
X(f)"
Z(t)
= e
¥t
T CTQC
(2.23)
Consequently, the solution of the matrix Riccati equation will be obtained by (2.15) as t-> °°. In the next section, we shall show that two matrices X(t) and Z(t) are obtained from the similarity transformation which reduces ¥ into the Jordan canonical form. Throughout . the remaining discussions, the increasing time variable t will be used instead of t.
3. SOLUTION OF THE MATRIX RICCATI EQUATION
We begin with two lemmas for our main result.
Lemma 1. The eigenvalues of the matrix ¥ of (2.21) are symmetric with respect to the imaginary axis.
Proof. Let I be a 2n x 2n matrix given by
xi =
-I 0
(3.1)
where I is an n x n identity matrix. Then it is obvious that
-1 T
I = I
1 1
T T
I WI = -¥
(3.2) (3.3)
Since the eigenvalues of a matrix are unchanged by either similarity trans- formations or transposing, the eigenvalues of W occur in pairs with opposite
signs.
Lemma 2. Let
T =
All A12
A A LA21 A22J
(3.10
where A , A , A , and A are n x n, m x m, nxm, and m x n matrices, respectively. Assume that A.^ is non-singular. Then
-1
T T
11 12
T T
L 21 22J
-1
(3.5)
exists iff A - A A " A is non-singular, and is given by
-1
-1
Tll : All " + All " ^2 T22 A21 All T12 = "All ' A12 T22
-1
(3.6)
T = -T A A 21 22 21 11
-1
T22 = (A22 " A21 All A12)
7
where submatrices T , T , T , and T are the same sizes as A , A22' A12' and A21' resPectively-
Proof. Suppose A - A Ann" ^o is non-sin8'ular' From TT " = I.
All Tll + A12 T21 = Xn
(3.7)
lll T12 + ^2 T22 " °
A21 Tll + A22 T21 " °
A21 T12 + A22 T22 = Im Premultiplying the second equation by A A and subtracting from the
fourth,
T22 " (A22 " A21 All " A12) '
And hence
T12 = ^1 \2 T22
Similarly, from the first and third, T21 = "T22 A21 ^1 '
-1
-1
Tll == All " + All " h2 T22 A21 All
-1
Hence T is a right inverse of T. In a similar way we can show that T is also a left inverse of T. Conversely, suppose T is non-singular. Then
0 ? det
A,
^1 "12
LA21 A22J
= det
A22~A21 hi A12.
-1
-det (A^) det(A22-A21A11- A^)
(3.8)
Therefore, A - A A A is nonsingular.
Assuming that W has distinct eigenvalues, we now turn to constructing the similarity transformation with which W of (2.2l) is reduced to the
8
diagonal form. By lemma 1, W has n eigenvalues with positive real parts and n eigenvalues with negative real parts. If we construct the similarity- transformation S, whose columns are the right eigenvectors of ¥, such that the first n columns are the right eigenvectors of W corresponding to the n eigenvalues of W with positive real parts, then
s_1ws =
1
X.
*A
2n
(3.9)
,-1
where, of course, the rows of S " are the left eigenvectors of W. O'donnell [5] and Potter [6] have shown that if we partition S into four n x n submatrices,
S =
Sll S12
.321 S22J
(3.10)
then the solution of the matrix Riccati equation is given by
-1
K - S21 S±1 provided that S is non-singular.
(3.11)
We shall now remove the restriction that W has distinct eigenvalues and present an alternative proof to that of Martensson [7] to show that (3.11 ) holds for the case in which ¥ has nondiagonal Jordan form.
Let us consider the general Jordan form J with the similarity transfor- mation S, which will be constructed later, S-1¥S = J
where, J =
"i ,
(3.12)
'J
m
and in which
m.
A. 1
1
A. 1
l
". A.
, (m. x m. )
(3.13)
If we denote S by its column vectors s , . . . , s , then
¥ Ls1, s2,
> S2n] = [S1' S2' •••• S2n] J
(3.1M
For any Jordan block J , if i. > 1 then, dropping subscripts of s for
m. l
l
convenience,
(1) . _(D
Ws
= sv 'A.
l
w«<2> ■ *(l) + > .(2) ws = s + A . s
(3.15)
or
TT (m. ) _ (m.-l) (m. )
¥si=si +A.si
(W - A.I) s(l) = 0 i
(W - A.I)2 s(2) = 0
(3.16)
/ vin. (m. )
(W-A.IJis l =0
Thus, s , 1 < j <m., is a principal vector of degree j associated with
, 0 . 0 . ... . , n , (l ) (m. )
A.. Since S is nonsmgular, the m. principal vectors s , ..., s l
are independent. If we consider all the Jordan blocks J , 1 < i < £, we
m. -
l
have the 2n principal vectors s , s ,
, s^ which constitute the
2n
similarity transformation S in (3.12). Further, we rearrange the columns of S such that the first n vectors are the principal vectors corresponding to the n eigenvalues with positive real parts of W.
10
Let us now split J of (3.12) into
J =
: J2^
|
• \ |
(3.17)
where J and J are diagonal with the eigenvalues of W with positive and negative real parts respectively, and E and E are matrices with one's on the super-diagonals and zeros elsewhere. To associate with (3.6), let us partition S and S as
|
sll |
S12 |
|
|
s = |
||
|
S21 |
S22 |
|
|
- |
||
|
_ |
||
|
s-1- |
Tll |
T12 |
|
T21 |
T22 |
(3.18)
where all submatrices are n x n. Note that, by lemma 2, it is necessary to assume that S and S - S S S are nonsingular.
Then
|
wt e = e |
SJS 1t |
|||||||
|
= Se^s'1 |
||||||||
|
0m Sll S12 |
e x 0 |
e^ 0 |
— T 11 |
T 12 |
||||
|
S21 S22 |
n J2t 0 e ^ — — |
0 ^ |
T 21 |
T 22 |
S11eJlteEltT11 + S12eJ2teE2tT21
S11eJlteEltT12 + S12eJ2teE2tT22
q eJlteEltT + S eJ2teE2tT S eJlteEltT + S eJ2teE2tT b21 L 111 22e d 6 ^ 121 b21 ± X ^2 b22 d & * X22
(3.19)
11
Since the elements of Jp have negative real parts, e 2 -^ 0 as t + » and thus the terms having e 2 in (3.19) vanish as t -*■ °° . Therefore,
X(t) Z(t)
¥t
T C QC
„ J_ t iii_ tm ^ J. t E., t_
Slie X 6 1 Tll Sll6 X e 2 T12
21 X e X 111 b21 ^ X12
ri
T C QC
Su eV eV (T1X + T12 CTQC ) S21 eV eV (T^ + T CTQC)
(3.20)
Therefore the solution of the matrix Riccati equation is given by
K = lim Z(t) X_1(t) t-*»
= lim [S eJlt e¥ (T +T CTQC ) ] [S eJlW (T +T CTQC) I"1
t-K»
■ S21 Sll
-1
(3-21)
provided that the indicated inverse exists. However, we shall show that only the nonsingularity of S and S - S S S p is needed for the
stationary solution of the matrix Riccati equation. From (3.12),
T T
11 12
T T
21 22
-1 T -A BR B
T T
C QC A
Sll S12
21 22
= J
(3.22)
Equating the (2,1) elements on each side of (3.22),
(-T21 A - T2£1 CTQC) Su + (T BR^B1 + ^ AT) 8^ = 0
12 or
■ -T21 ASn + T22 AT S21 + T21 BR^ S21+ T22 CTQC S^ = 0 (3.23)
Substituting T21 = - T22 S21 S^"1 into (3.23),
T22S21Sll'1+ T22ATS21 " T22S21Sll"lBR"lBTs21 + ^2°^° Sll=°
(3.21+)
Premultiplying by T22" and post multiplying by S _1 ,
S21Sll"lA + ^^l3!!"1 " ^l5!!"1^"1^^!3!!"1 + <?<* = ° .0.25)
which is to be proved.
13
k. COMPUTATIONAL METHODS AND NUMERICAL RESULTS
The major three steps for computing the solution of the matrix Riccati equation are as follows:
1) Compute eigenvalues and eigenvectors of W,
2) Compute principal vectors of ¥ if W is defective,
3) Compute S S " from S.
For computing the eigensystem, we used two methods: (i) the QR algorithm [9] with an inverse iteration routine [10], and (ii) Jacobi-like method [ll] for finding the eigenvalues and eigenvectors. Their perfor- mances are satisfactory and two numerical test examples are given in Test Result 1 and Test Result 3. When the matrix W is nondiagonable, Test Result 2, the first method requires the evaluation of the principal vectors as well [12]. This poses some numerical difficulties as we are essentially trying to solve the ill-conditioned system of equation
(W - Al)k x = x., k > 1 where ^ is an approximation of true eigenvalues A. Thus, the Jacobi-like method is more attractive in this case for providing an approximate eigen- system that is suitable for engineering purposes.
Appendix I contains the listings of the ALGOL programs used on the B-6500 Appendix II contains the listings of two ILLIAC IV programs, namely, the inverse iteration routine and a routine for solving a system of complex linear equations, both written in GLYPNIR. ILLIAC IV programs for the QR algorithm [13, 1^] and the Jacobi-like algorithms [15, l6] are already available.
Test Result 1.
ik
A =
|
1 |
0 |
0 |
0 |
0 |
|
1 |
0 |
-1 |
0 |
0 |
|
0 |
0 |
-1 |
0 |
0 |
|
0 |
0 |
1 |
0 |
-1 |
|
0 |
0 |
0 |
0 |
-1 |
p =
|
0 |
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
0 |
0 |
B =
|
1 |
0 |
0 |
|
0 |
0 |
0 |
|
0 |
1 |
0 |
|
0 |
0 |
0 |
|
0 |
0 |
1 |
Q =
|
0 |
0 |
0 |
0 |
0 |
|
0 |
10 |
0 |
0 |
0 |
|
0 |
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
10 |
0 |
|
0 |
0 |
0 |
0 |
0 |
c =
|
1 |
0 |
0 |
0 |
0 |
|
0 |
1 |
0 |
0 |
0 |
|
0 |
0 |
1 |
0 |
0 |
|
0 |
0 |
0 |
1 |
0 |
|
0 |
0 |
0 |
0 |
1 |
R =
|
1 |
0 |
0~ |
|
0 |
1 |
0 |
|
0 |
0 |
1 |
The numerical example was tested by the QR algorithm and the inverse iteration method. The ten distinct eigenvalues of W are
± 1.
± 1.728760U77185 ± 1.577533767573i,
± 1.35319578U0U6 ± 1.1537U989928Ui. We choose the five eigenvectors corresponding to the five eigenvalues with positive real parts and computing K = S S from them gives the solution K as
15
1.262782609
2.1+9^009759
-0.819173651
0.668267901
-O.UU3608958
2.1+9^009759 7.^ 3 5^51161+
-1.8257^1858 1.1229101+32
-0.668267901
-0.819173651 -1.82571+l858 1.63831+7303 1. 8257^1858 -0.819173651
0.668267901 1.1229101+32 1.8257^1858 7.1+ 35^51161+ -2.1+91+009759
-O.Ul+3608958 -0.668267901 -0.819173651 -2.1+91+009759 1.262782609
-1 T All the eigenvalues of G = A - BR B K have negative real parts and hence
the above solution K is acceptable.
Test Result 2.
A =
-3 '2
-2 1
B =
C =
1 0 0 1
P=Q=
0 0 0 0
|
The mat |
rix — |
|||
|
3 |
-2 |
0 |
0 |
|
|
2 |
-1 |
0 |
1 |
|
|
¥ = |
||||
|
0 |
0 |
-3 |
-2 |
|
|
0 |
0 |
2 |
1 |
R =
CO
is defective, that is, W has four eigenvalues ±1 (double) and each +1 and -1 correspond to a quadratic elementary divisor, respectively. This example was tested both by Jacobi-like method and by a combination of the QR algorithm and the inverse iteration method.
16
The Jacobi-like method gives K as -10.0000081* 6.00000U3
6.000003^ -k. 0000027
and the QR and inverse iteration method give K as -9. 9999917 5.99999^7
5.9999953 -3.9999983 The exact solution is
~Vlo 6 6 -k
For computing the principal vector by the QR and inverse iteration, we chose the perturbation for the true eigenvalues as much as 10
Test Result 3.
A =
|
6 |
k |
k |
1 |
|
h |
6 |
1 |
k |
|
h |
1 |
6 |
k |
|
1 |
h |
h |
6 |
B =
|
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
1 |
C=R=
|
1 |
0 |
0 |
0 |
|
0 |
1 |
0 |
0 |
|
0 |
0 |
1 |
0 |
|
0 |
0 |
0 |
1 |
P=Q=
|
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
0 |
|
0 |
0 |
0 |
0 |
The exact eigenvalues of W are
±15, ±5, ±5, ±1 Since the matrix W is derogatory, +5 and -5 correspond to a two-dimensional subspace, respectively. The computed eigenvalues by Jacobi-like method are
IT
lU. 99999999930 k. 9999999991+9
-U. 99999999982 0.99999999997
-1^.9999999984 4. 99999999981 _1+. 999999999^8.
-0.99999999995
The computed solution K is given by
-1.999999999811 1.9999999999^2 1.999999999913 1.999999999884 -2.000000000015 -2.000000000000 2.0000000002T6 -2.000000000^07 -2.000000000378
-2.00000000036U 2.000000000509 2.000000000I180
-2.000000000015 2.000000000087 2.000000000^80
-2.000000000568
in -which the absolute error is less than 10
-9
18
LIST OF REFERENCES
[I] M. Athans and P. L. Falb, Optimal Control, McGraw Hill, New York, 1966
[2] R. E. Kalman, Contributions to the Theory of Optimal Control, Boll. Soci. Mat. Mexicana, 5 (i960), pp. 102-119
[3] R. S. Bucy and P. D. Joseph, Filtering for Stochastic Processes with Applications to Guidance, Interscience, New York, 1968
[h] A. G. J. MacFarlane, An Eigenvector Solution of the Optimal Linear Regulator Problem, J. Electronics and Control, 1^4- (1963), pp. 61*3-6^5
[5] J. J. O'donnell, Asymptotic Solution of the Matrix Riccati Equation of Optimal Control, Proc . Fourth Annual Allerton Conf . , 1966, pp. 577-586
[6] J. E. Potter, Matrix Quadratic Solutions, J. SIAMAppl. Math., ll+ (1966), pp. U96-5OI
[7] K. Martensson, On the Matrix Riccati Equation, Inf. Sci., 3 (l97l), pp. 17-^9
[8] D. L. Slotnick, et al, The ILLIAC IV Computer, IEEE Trans, on Computer. C-17 (1968), pp. 7^6-757
[9] R. S. Martin, G. Peters, and J. H. Wilkinson, The QR Algorithm for Real Hessenberg Matrices, Num. Math., lk (1970), pp. 219-231
[10] J. M. Varah, The Calculation of the Eigenvectors of a General Complex Matrix by Inverse Iteration, Math. Comp. , 22 (1968), pp. 785-791
[II] P. J. Eberlein and J. Boothroyd, Solution to the Eigenproblem by a Norm Reducing Jacobi Type Method, Num. Math., 11 (1968), pp. 1-12
[12] H. J. Bowdler, R. S. Martin, G. Peters, and J. H. Wilkinson,
Solution of Real and Complex Systems of Linear Equations, Num. Math., 8 (1966), pp. 217-234
[13] M. Ogura, A Parallel Scheme for Reducing a Real Matrix to the Upper Hessenberg Form, Center for Advanced Computation Document No. 11, University of Illinois, 1971
[ill] M. Ogura, The QR-Algorithm, Center for Advanced Computation Document No. 12, Univ. of 111., 1971
19
[15] W. Bernhard, ILLIAC IV Codes for Jacobi and Jacobi-like Algorithms, Center for Advanced Computation Document No. 19, Univ. of 111., 1971
[l6] A. H. Sameh, On Jacobi and Jacobi-Like Algorithms for a Parallel Computer, Math. Comp. , vol. 25 (l9Tl), pp. 579-590
20
APPENDIX I
ALGOL PROGRAMS
21
I TWO ALGflL PROCEDURES' RlCCATU AND RICCATI2. FOR COMMUTING THE OOOOOlOO?
I! SOLUTION OF THE MATRIX RICCATI EOUATION ARE LISTEO BtLOWt FOR 00000200?
% "RICCATlU. THE PROCEDURE "EIGENw (BY EflERLEIN AND bUOTHKoyD) WAS 00000300?
% USEU WHILE A COMBINATION OF "HQR" (BY MaRTINi ET AL) AND 00000400?
I "INVITErATION" WAS APPLIED TO «RICCATI2% LINEAR EqUATIUn SOLVER OOOOO5OO? % (BY BOWDLER* ET AL ) WAS ALSO USED TO PROVIDE THE INytRSE OF A MATRIX 00000600?
% ANO TO CALCULATE THE PRINCIPAL VECTORS fROM ( ( W-LAMflUA*! )**K>*X«0. 00000700?
00000800?
00000900?
PROCEDURE RlCCATIl(N.L»M»A»B#C»P»Q»R«TMX»W»K)l 00001000?
VALUE N»L#M,TMXI 00001100?
INTEGLR N,L.M.TMX| 00001200?
ARRAY A.e(C.P»Q»R»W.K£i»n| 00001300?
COMMENT SOLVES THE MATRIX RlccATI EQUATION . OOOOUOO?
KA*(AT)K-KR(RINV)(BT)K*(CT)QC«o, 00001500?
A Is N*N» R IS N*L» C IS M*N. P AND 0 ARE M*M. R IS L*L OOOOUOO?
MATRICES* RESPECTIVELY* FORM 2N*2n MATRIx W ANU COMpUTE 00001700?
EIGENVALUES ANU EIGENVECTORS OF W* WHERE FOUR N*N BLOCKS 00001800?
Op w ARE C 1 • 1 >«-A» Cl»2)«8(RlNV>(BT)» ( 2» 1 )■( C ' >QC» AND 00001900?
(2»2)"(AT). REARRANGE EI GEN VEcTOH ( Pr INC I PAL VEcT0R> MATRIX 00002000?
SO THAT THE pIHST N COLUMNS CORRESPOND TO TO THE EIGENVALUES 00002100?
WlTH P0SITIVE REAL rArTS. dENoTINq ThE UPP^R HaLf uF N COLUMN 00002200?
VECTORS BY SMTXl AND THE LOWER HALF BY SMTX2* I HE SOLUTION 00002300?
K IS THEN GlvEN BY «■( SMTx2) ( SMTxl INv ) . THE PRUCEOUrE 00002400?
"FI6EN" IS USED For EIGENPrOBlEM qF Wi 00002500?
00002600?
REGIN COMMENT FIRST OECLARE PROCEOuRESl 00002700?
00002800?
PROCEUURE SETUP(NtL»M.A»B»C»Q»RlNV.W)| 00002900?
VALUE N.L.M, INTEGER N.L»M) 00003000?
REAL ARRAY A»B»C»Q.RINV»W[l#l]j 00003100?
COMMENT THE PROCEDURE SETS UH ThE 2N*2N mATRIX W, rInV Is INVERSE OF RI00003200?
«EGIN 00003300?
REAL ARRAY RINVBTC1 lL»l IN].QCC1 IM»1 IN] » 00003400?
REAL SUMj 00003500?
INTtGER I.J.KJ 00003600?
FOR Ij"l STEP 1 UNTIL N DO 00003700?
FOR J l ■! STEP 1 UNTIL N DO 00003800?
BtGIN 00003QOO?
W[I,J]I«-A[ I. J]} 00004000?
W[N*I,N*J]IbA( J.I1I 00004100?
ENO) 00004200?
FOR 1 1 » 1 STEP 1 UNTIL L 00 00004300?
FOR J»«l STEP 1 UNTIL N DO 00004400?
BEGIN StJMl.O.Ol 00004500?
FOR K,«l STEP 1 UNTIL L DO 00004600?
S|jM|«SUM*RINVCl»K]*BC J#KJ1 00004700?
RINvBTC I.J]l«SUM» 00004800?
ENOi 00004900?
FOR I|M STEP 1 UNTIL N 00 00005000?
FOR Jl»l sTEp 1 UNTIL N 00 00005100?
BtttIN SuHlaO.Ol 00005200?
FOR Ki«l STEP 1 UNTIL L DO 00005300?
SuH«"SUH*B[IsK]*RINVBTtK#J]| 00005400?
Wt I,N*J]l«SUMl 00005500?
END| 00005600?
FOR I|»l STEP I UNTIL M DO 00005700?
22
FOR Jl«l STEP 1 UNTIL N DO
BEGIN SuMl.O.Oi
for K|«i step i until m do
SuMl"SUH*OtItK]*CCK.JJI QCC I # J) I -SUM l ENO| TOR I|«l STEP 1 UNTIL N DO FOR J|«l STEP 1 UNTIL N DO BtG]N SuNlsO.Ol
FOR K|»l STEP 1 UNTIL N DO StjMI»SUM*CtK. i J*«C t K. J3 1 W[N*I.J]|"SUMI ENDI ENO OF SETUP|
COMMENT INSERT THE PROCEDURES INNPRODtDECOMPOSE. AND ACCSOLyEi
PROCEUURE lNVERSE(N»AtX)l
VALUE N| INTEGER Nj
DOUBLE ARRAY A»XC 1 • 1 1 I
COMMENT THE PROCEDURE COMPUTES THE INVERSE OF AN N*N rUl
UNSyMMcTRIc MATRIX A. USINq ThE PROCEOURgS "DECQMpUsE* AND "ABSOLVE". WRITING ON X| BEGIN
INTtGER I.J.IT.D2I REAL Oil
ARRAY AA'BB.BU I N* 1 iNJ.P[llNll FUR I • * 1 STEP 1 UNTIL N DO 1 UNTIL N Oo
BCl»J]l«ltO ELSE B(I»J]l«O.OI 1 UNTIL N DO 1 UNTIL N DO
FuR Jl«l
IF J«I
FUR I,«l
FUR J t -1
STEP THEN STEP STEP AACT.Jlt-ACI.J1l DtCOMpOSE(N»AA.Dl.02»P)l AtCSULVE(N.N.A,AA.P»B.X»BB»IT)| ENO INVERsEl
COMMENT iNjSfRT THE PRnCEDURE EIGENi
PROCEUURE MuLT(Nl.N2iNJ»A.B»C)l VALUE N1»n2.N3I INTEgEH N1*N2»n3> DOUBLt ARRAY A.B.Ctl»lJI REGIN
INTtGER I.J.KI
REAL SUMI
FOR I l » 1 STEP i
FOR J t • 1 STEP 1
BEGIN SUM|«0.0|
FUR K,"t STEP
SUMI»SUM*ACI»K]*BCK,J]| CI I .J] l«SUMI ENO» FNO MULTl
ARRAY HiNv»SMTXl,SMTX2lllN.llN].WW.TEMP,VECT[llN*N»llN*NJ| INTLGEH ArRAy search(1»nii INTtGER N?» I. J- INDEX'
UNTIL UNTIL
Nt N3
DO OO
1 UNTIL N2 OO
00005800? 00005900? 00006000? 00006100? 0000*200? 00006300? OOOO64OO? 00006500? 00006600? 00006700? 00006800? 00006900? 00007000? 00007100? 00007200? 00007300? 00007400? 00007500? 00007600? 00007700? 00007800?
00007900? 00008000? 00008100? 00008200? 00008300? 00008400? 00008500? 00008600? 00008700? 00008800? 00008900? 00009000? 00009100? 00009200? 00009300? 00009400? 00009500?
00009600? 00009700? 00009800? 00009900?
00010000? 00010100? 00010200?
00010300?
00010400? 00010500?
00010600? 00010700? 00010800? 00010900?
00011000? oooiuoo?
00011200? 00011300? 00011400?
23
N2I«N+N| 00011500?
INVtRSE(L.R»RlNV)l 0001U00?
SETUP <N,L#M»A»B»C»Q»RINV»W)J 00011700?
TOR It*l STEP 1 UNTIL N2 DO 00011800?
FOR Jl "1 STEP 1 UNTIL N2 DO 00011900?
TEMPI I. JJ t«M[ I.J] I 00012000?
EIQtN(N2»TMX,TEMP,VECT)l 00012100?
COMMENT SEARCH N COLUMNS OF VECT CORRESPONDING TO N EIGENVALUES 00012200?
WITH POSITIVE REAL PARTS. SMTXl AND SMTX2 ARi IT* UPPER 00012300?
AND LOWER HALVES* RESPECTI VELY| 00012400?
INDLX|.0l 00012500?
FOR I » » 1 sTEp 1 UNTIL N2 DO 00012*00?
I* TEMP[I»I1 GTR 0 THEN 00012700?
BLGIN 00012800?
INDEX|«INDEX*1 , 00012900?
SEArChI INDEX] tail 00013000?
ENO| 00013100?
FOR J|»1 STEP 1 UNTIL N DO 00013200?
FOR I, al STEP 1 UNTIL N DO 00013300?
BEGIN 00013400?
SMTxltliJ) |«VECTC I » SEARCH! J]]j 00013500?
SMTX2tI.J]l«VEcTtN*I»SEARcHCJ}]l 00013600?
ENDI 00013^00?
INVLRSE(N, SMTXl. TEMP)I 00013800?
MULHN»N»N»SMTX2.TEMP»K)| 00013900?
END RiCCATH» 00014000?
00014100?
00014200?
PR0CEUURE RICCatI2(n.L»M»a»B»C»P»0»R»W»K»0EFECt»DER0(»AIE)I 00014 300?
VALUE N.L.Mf 00014400?
INTEGtR N.L.Ml 00014500?
ARRAY A»3.C.PfQ»R*H*Ktl'l)l 00014600?
INTEGtH ARRAY DEFECT, DEKUGATEt 1 J J 00014700?
COMMENT SnLvES THE MATRIX RICCATI EQUATION 00014800?
KA*(AT)K-KB(HINV)(BT)K*(CT)QC«0. 00014900?
A Is N*N. b IS N*L» C IS M*N. P AND Q ARE M*M, R IS L*L 00015000?
matrices* respectively. form 2n*2n matrix w and compute 00015100?
eigenvalues anu eigenvectors of w, where four n*n blocks 00015200?
OF W ARE (l.l)«-A. (1.2)«B(RINV)(BT)* (2*1)«(CU0C AND 00015300?
(2»2)«(AT). REARRANGE El GEnVEcTOR( PRlNC I PAL VECTOR* MATRIX 00015400?
SO THAT THE FIRST N COLUMNS CORRESPOND TO TO THE EIGENVALUES 00015500? WITH POSITIVE HEAL PARTS. DENOTING THE UPPER H*LF Up N COLUMN 00015600?
VECTORS BY SMTXl AND THE LOWER HALF BY SMTX2, THE SOLUTION 00015700?
K IS THEN GIVEN BY K«( SMTX2) < SMTXl INV ) . THE PROCEDURES 00015800?
"hShLD"* "hOR** AnD "INVITEKATION" ARE USEq FOK EH»ENPRORlEM 00015900?
OF Wl 00016000?
00016100?
REGIN COMMENT FIRST DECLARE PROCEDURES! 00016200?
00016300?
PROCEUURE SETUP(N»L*MtA*B,C»Q*RINV*W)| OOOI64OO?
VALUE N.L.Mf INTEGER N.LtM* 00016500?
REAL ARRAY A » B . C » Q . RIN V • W [ 1 i'l ) I 00016600?
COMMENT THE PROCEDURE SETS UP ThE 2N*2N MATRIX W. RINV IS INVERSE Of R! 00016700?
QEGIN 00016800?
REAL ARRAY R I NVBT[ 1 I L« 1 I N J . QC [ 1 I M» 1 1 N ] I 00016900?
REAL SUM1 00017000?
INTLGER I.J.KI 00017100?
2k
FOR I t -1 STEP 1 UNTIL N DO TOR J|»1 STEP 1 UNTIL N DO BtGIN
M[I.J]I«-ACI»J]|
WtN+I.N*J]l«ACJ#ni
ENDI
Fo« i»«i step i until l Do
FOR Jial STEP 1 UNTIL N DO BtGIN SuM|aO»0|
FOR K|«l STEp 1 UNTIL L DO SuM|«SUM*RJNVt I»KJ*B[ J»K] J
RINvBTtI»J]l«SUMI
END|
FOR I|»l STEP 1 UNTIL N 00 FOR J t"l STEp 1 UNTIL N DO BtGIN SuMl-0,01
FOH Kt«l STEP 1 UNTIL L 00
SuM,«SUM*BU.K)*RlNVBT[K.JJ| W[I.N*JJ|«SUmI ENO| FOR Iial STEP l UNTIL M 00 FOR J|"l STEP 1 UNTIL N DO BLqIN SUMI.O,0|
FOR Kf.l STEP 1 UNTIL M DO
SUM|.SUM*Q[ I.K]*C[K»Jj I «CC If j] I»SUMJ ENDl FOR Il»l STEP 1 UNTIL N DO FOR J t "1 STEP 1 UNTIL N DO BtGIN SUMl«O.0l
FOR K|b1 STEP I UNTIL M DO SuM|«SUM*C(K.n*UC[K»Jj; W[N*I,j) l«SUMI END>
end of setupi
comment Insert the procedures innprod. decompose* and accsolvei
proceuure inverse(n»a#x)|
VALUE N| INtEgEH N|
nOUBLt ARRAY A.XU.lJI
COMMENT ThE PROCEOURE COMPUTES THE INVERSE OF AN N*N RtAL
unsymmetric hatrix a. using The procedures "DECQmpu$e" waCcsolvE". writing on XI
BEGIN
INTtGER I.J.IT.D2I REAL nl |
ARRAY AA'BB.BUtN.l fNJ»Ptl tNJ| FUR I i - 1 STEP 1 UNTIL N DO
DO
B[I»J]nO«OI
AND
FuR J I -l
IF J-I
FUR I,-!
FUR J i ■ 1
END
step
then step
STEP
AAt !ij]|iA[I,jii ntCnMp0sE(N.AA,0l.02#P)| ACCSULVE(N#N.A.AA.P»B#X»BB»IT)|
i until n
fl[I. J] i-l .0 ELSE 1 UNTIL N DO 1 UNTIL N
00
00017200?
00017300?
00017400?
00017500?
00017600?
00017700?
00017*00?
00017900?
00018000?
00018100?
00018200?
00018300?
00018400?
00018500?
00018600?
00018700?
00018800?
00018900?
00019000?
000i9i00?
00019200?
00019300?
00019400?
00019«00?
00019600?
00019700?
00019800?
00019900?
00020000?
00020100?
00020200?
00020300?
00020400?
00020500?
00020600?
00020700?
00020800?
00020900?
00021000?
00021100?
00021200?
00021300?
00021400?
00021500?
00021600?
00021700?
00021800?
00021900?
00022000?
00022100?
00022200?
00022300?
00022400?
00022500?
00022600?
00022700?
00022800?
25
00022900?
Comment InseRt the procedures hshld and hori ooo23ooo?
00023100?
PROCEDURE MuLT(Nl»N2*Ni'A»B«C)l 00023200?
VALUE N1»N2*N3| INTEGER N1.N2.N3j 00023300?
DOUBLE ARRAy A.B.CU.Ul 00023400?
REGIN 00023500?
!NTtOER I»J»K» 00023600?
REAL SUMS 00023700?
FOR Il«l STEP 1 UNTIL Nl DO 00023800?
FOR Jl«l STEp 1 UNTIL N3 00 00023900?
BEGIN SUM|«0,0I 00024000?
FUR K|«l STEP 1 UNTIL N2 DO 00024100?
SUM i * SUM* A [ I»K]*B[K«J)| 00024200?
CI I • Jl l.SUMI 00024300?
ENOI 00024400?
END MULTl 00024500?
00024600?
PROCEUURE lNVlTERATION<N»w»X»KK»P>| 00024700?
VALUE N, TNTEGER N.KKj 00024800?
DOUBLE ARRAy W.Xtl.lli 00024900?
ARRAy P[l]l 00025000?
CnMMENT SnLVES WX.O BY INVERSE ITERATION* WHERE W IS AN N*N REAL 00025100?
MXTRlX. WRITING SOLUTION ON X. THUS. THE PROCEDURE COMPUTES 00025200?
ThE EIqENvECTOKS CORRESPONDING TO THE GIVEN EIGENVALUES. 00025300?
N Is THE ORDER OF M AND KK IS THE NUMBER OF iTtRATiONS 00025400?
NEEqEO FOR ITERATING EACH EIGENVECTOR! 00025500?
BEGIN 00025600?
INTfcGER I.J»D2»R.L»MAXIN0EX| 00025700?
DOUBLE DUNAXVALi 00025800?
DOUBLE ARRAY MN( 1 1 N. i | NJ »B.B8[ i | N» 1 | U | 00025900?
LABEL INVIT.OUTI 00026000?
FOR Il"l STEP 1 UNTIL N DO 00026100?
FOR Jl ■ l STEP 1 UNTIL N DO 00026200?
ww [ I » J ] t»ri t I » J 3 1 00026300?
DECUMpOsE(N»WW»Dl»02»P)l 00026400?
FOR Il«l STEP 1 UNTIL N DO 00026500?
BII,1]»81,0» 00026600?
KK I "1 1 R| ■ I j 00026700?
INVITI 00026800?
ACCSOLVe(N»R»W»WW»P»B»X.BB*L)I 00026900?
MAXVALUOj 00027000?
FOR Ii"l STEP 1 UNTIL N DO 00027100?
If ABs(X(I*l]> GTR MAXVAL THEN 00027200?
BEGIN 00027300?
MAXvALl«ABS(Xr I#i])l 00027400?
MaXiNDEXI«P 00027500?
ENDi 00027600?
MAXvA|_l«X[MAXlNDEX»ni 00027700?
FOR I|«l STEP 1 UNTIL N DO 00027800?
XlI»l]l«X[Iil]/MAXVALI 00027900?
FOR III} STEP 1 UNTIL N 00 00028Q00?
I* ABS(BCI»1]-Xtl»l]) GTR P"12 THEN 00028100?
BLGIN 00028200?
IF KK GTR 10 THEN GO TO OUTj 00028300?
FOR J|-l STEP 1 UNTIL N DU 00028400?
BtJ,l]l"X[j.i ]l 00028500?
26
KK»»KK+U GO TO INVITI END)
nuTi
ENO INVREsElTi
INTtGER I.J»JJ»CNT.N2»INDEX#ITC0UNT|
ARRAY TEMpfORGNtSHTX*veCTORtt|N«N«l|N4N]*RlNViSHTXUSHTX2tl|N«l|N]»
EIGR»EIGI»!NTCH(1IN4N)I REAL PERTURB. MPUiERl INTtGER ARRAY SEARCHllIN)! N2l"N*N|
INVk.RSE(L.R.RlNV)l SETUP(N,L,M,A»B,C.Q.HINV.H)| TOR I|«l STEP 1 UNTIL N2 00 FOR J|"l STEP 1 UNTIL N2 00
HSHLO(i,N,,TEMP)| HQR(N2»TEMP.I
N2 N2
DO 00
N2 00
H
EIQR.EIGI* CNT)| FOR I | »1 STEP 1 UNTIL N2 00 FOR J I - 1 sTEp 1 UNTIL N2 00
ONGHCli J]l«NCItJ]l COMMENT COMPUTE EIGENVECTORS OR PRINCIPAL VECTORSi FOR Jjl-1 STEP 1 UNTIL N2 00 BEGIN
FUR I t»l STEP l UNTIL
FUR Ji-i STEP 1 UNTIL *[I»J] ««ORGW[I,J]l
FUR I » - 1 STEP 1 UNTIL
•»tIiI]l"WU.i]-ElGR[ JJ]*PERTURB|
JJ«DEFECT[JJ] then BEGIN
FOR Ii-1 STEP 1 UNTIL N2 00
W[I.I J»«Wti,IJ*pERTURB*MpLIERl
MuLTCN2.N2»N2»W»W,TEMP)l
Ff)R Ital STEP 1 UNTIL N2 00 FOR Jlal STEP 1 UNTIL N2 DO *t I«J]i«TEMP[ I#J)I END | IF JJaOfROGATEUJ] THEN
FOH Ii-l STEP 1 UNTIL N2 00 Mf I • 1 1 |»WC 1. 1 1 ♦PERTURB* J Jl
INVITFRATI0N(N2#W.VECT0R»ITC0UNT»INTCH)» FUH I t a.j STEP i UNTIL N2 DO
FUR Jiaj.2 DO
SMTx(l»JJ]|avEcTORtI#l Ji END JJ>
COMMENT SEARCH N COLUMNS OF SMTx CORReSpONOiNG TO WITH POSITIVE REAL PARTS, SMTxl AND SMTx2 AND LOWER HALVES* RESPECTI VELYl INDtXiaQi
FOR ll«l STEP 1 UNTIL N2 00 l> ElGRfl] GTR 0 THEN BLGIN
iNOFXtalNUEX^l | SEArcmC INDEX] I.I > ENDl
N tlGENVALUES ARt ITS UPPER
00028600? 00028700? 00028800? 00028900? 00029000? 00029100? 00029200? 00029300? 00029400? 00029500? 000296O0? 00029700? 00029800?
00029900? 00030000? 00030100? 00030200? OOO3O1OO? 00030400? 00030500? 00030600? 00030/00? 00030800? 00030900? 00031000? OOO3UOO? 00031200? 00031300? 00031400? 00031500? 00031600? 00031700? 00031800? 00031900? 00032000? 00032100? 00032200? 00032300? 00032400? 00032500? 00032.00? 00032700? 00032800? 00032900? 00033000? 00033100? 00033200? 00033300? 00033400? 00033500? 00033600? 00033700? 00033^00? 00033900? 00034000? 00034100? 00034200?
27
FOR Jl«i STEP 1 UNTIL N DQ 00034300?
FOR Ii"l STEP 1 UNTIL N DO 00034400?
BEGIN 00034500?
SMTXltI.J]l«SMTX[I»SEARCH[Jlll 00034600?
SMTX2(I,J)I«SMTX(N+I'SEARCH(J])J 00034700?
CNDI 00034800?
INVtRSE(N,SMTXl»TEMP)l 00034900?
MULT(N»N»N»SMTX2.TEMP»K)| 00035000?
END RlCCATl2' 00035100?
28
APPENDIX II
ILLIAC IV PROGRAMS
29
* TWO ILLIAC IV GLYPNIR PROGRAMS. "ClINSYs" AND "INVITE-RATION" ARE
f LISTED BELOW, "CLlNSTS" IS USED TO OBTAIN THE INVERS* OF A MATRIX
% AND ThE SOLUTION Or A COMPLEX SYSTEM OF UNSYMMETRIC LINEAR EQUATIONS.
% MNVITERATION" IS USED TO CALCULATE THE PRINCIPAL VECTORS OF DEGREE
* K FKOM ((WLAMBDAM)**K)*XnO.
SUBROUTINE cLINSyS(CINT N.ClNT R#PCPOINT A'PCPOINT x»PCPOINt B)J BEGIN x THE SUBROUTINE SOLVES A COMPLEX SySTEM OF UNSyMMETKlC LINEAR
x Equations ax»r. a is n*2n and b is n*2R matrices* respectively.
* ThEir qEnErAi EiEMENTS ArE At I »2J-H*I* At I »2J] AND BCp2j-ll*
* I*B[I,2J], RESPECTIVELY, R IS THE NUMBER OF RIGHT-HAND SIDES,
SUBROUTINE INNpROOf CREAL CWCREAl C2»PREAL AK»PREAL BK»C«EAL Dl#
CREAL D2)»
begin * accumulatEs the sum of products and aods it 10 the initial % value (c1.c2j.
CHEAL Sl2|
S12,»CUC2|
Si2l»Sl2*ROwSUM(AK*BK)i
01 I .Si 2 t
D2I.S12-D1I
END)
SUBKgUTTNE CDEC0MP0SE(CiNT N»PCPnlNT A'CREAL DErR'CRtAL i>etI»
CI BEGAN % ThE COMPLEX UNSy
where l is lower
MATRICES. AND Gv AND ITS GENERAL KEEpS THE RECOSO A, THE OETERMIMA
C1NT I» j.K.L.P.PPj
PREAL VICTOR ATI6AJI
PKEAL Z7I
CKEAL X.Y.Z.V.W.H.HHI
LABEL FAlLl MUDEI«TRUEl
MUDE»«TRUE AND pen lss LUUP II»1.1»N DO
lNNpRoD(0»O.AtI]»AtI DLTRI.ii DETII.Oj DETE LUOP Kl-l.l.N 00 BLGIN
Lf«Kl Pl«K*K| PPI«P-
LOOP II«K»1.N DO BEGIN
MoDEI«TRUEl
MoDE»«TRUE AND PEN
Innprodc-graboneca InnpRod(-h,-hh.rtl Innprodcgraroneu
lNNpROD(H.HH»ACIl. Y|«-HJ
modei«true» m0de1.true and pen mode»«truej
NT DETE»CNP0INT INT) I
MMETRIC MATRIX A IS DECOMPOStD INJO LU.
triangular and u is unit upper triangular erwritten on a, a is stored in N*2N array
ELEMENTS ARE At I*2j«l 1«I*A( I *2 J] t INTIM Of ANy INTERCHANGES mAdE To THE ROWS Op NT (0ETR*I#0ETI)*2**DETE IS ALSO COMPUTFD,
N + NJ
]ilNT(I].H)l 1*0 1
1) Zl»OJ
LSS N*N| tI].PP"i)*0*A[I]»ATtPP].H*HH^ (l»*A(l]).ATtP].X,HH)| tIJ.P-t),O.RTL(l»»ACI]).ATtPPj.H»HH)l
AT(Pl.H»HH)l
EQL PP"ll AC 1 3 f «X|
00000100? 00000200? 00000300? OOOOOaOO? 00000500? 00000600? 00000700? OOOOOBOO? 00000900? 00001000? 00001100? 00001200? 00001300?
oooouoo?
00001500?
oooouoo? 00001700?
00001800? 00001900? 00002000? 00002100? 00002200? 00002300? 00002400? 00002500? 00002600? 00002700? 00002800? 00002900? 0000300^? OC003100? 00003200? 00003300? 00003400? 00003500? 00003f00? 00003700? 00003800? 00003900? 0000*000? 00004100? 00004200? 00004300? 00004400? 00004500?
00004600? 00004700? 00004S00? 00004900? 00005000? 00005100? 00005200? 00005300? 00005400? 00005500? 00005600? 00005700?
30
00006800? 00006900? 0000?OOP?
So°oDI!:TTRRUULAN° PEN E°L P'U AI»-TI 00005800,
HODE„TRUE AND PEN.II SSSSfSSSJ
l' * GTR Z ™EN 00006200?
ENB0riN Z,"X' L"" ENDl oZZlll
lFDL,SKUK,T!SSCUTm AN° ^ LSS "*"* OOOoJjSSJ
BEGIN DETHI-OFTRI XXSXJJaXI
22 j »ACK3 1 AfKll-Afi 1»
iNTtLi'^TNTrKii
E *WT|Uj lalNiCKjl 00007200?
lNTfKli«il 00007300?
x.-JrJboUcackj.pp-u, SSJSfJSS?
li:«™*«t*un» Z°o7rl°o°ol
OETj,;JjDETI*V*OETRl SSSS^SSS?
IF ABS(DETR) GTR ABS(OETI) SSSSl?So?
Then whdetr else hi-detii §§§§si§S?
JErG?N"L°ARTEHLNL?EGlN DETE'"01 G° T° FAIU EN°' 0S0O83OO?
Lli IF ABSCW, GEO 1 THEN 22SS2S222-
BEGIN W!.W*0.0625, 000860
°lVATrTi:°a\06'5t SSSS55SSJ
DETE.-DETEMI 00008800?
EN0| G° T° L1» 00008900?
EN0 00009000?
BEGIN LABEL L2i 00009100?
L2, IF IbjCO LSS 0.0625 THEN SSSSSSSS5
BEGIN S^ETR-H, = °0'
DETKI-0ETR*l6l 00009500?
0ETII.0ET1M6I 00009*00?
GO TO l\] 00009700?
end,
MnOEl.TRUE ANO PEN LSS K*K-2| 00010*00*
lNNPROO(-GRABONE(AtK).PP-n.O.AtKl.ATtPP],H»HH)l 00010600?
INNP«nO(-H.-HH.RTL<J..A[K]).AT[Pl,V.MH)l 00010700?
|NlJpR0D(-6RAHnNE(ACK].P-n,0.RTL(l,.A|RPl).AT[PPl,H.MH)l OOOIOSOO?
INNPR00(H.HHPA[K3.AT[P3.H.HH)| OOOloSooJ
SaScVimucf SSSInSS?
MROEI-TRUE ANO PEN-PP-U OOOllJooJ
Mnon'THun oooiuoo?
0000980O? 00009900?
31
modemtrue ANO »EN«P-1I 00011500?
AfK]|»<W*X-V*Y>/Zl 00011600?
END| 00011700?
ENO| 00011800?
FAIL! 00011900?
ENOl «0r CDECOMPOSE. 00012000?
00012100? SUBHOUTINE CACCS0LVE(CINT N.CINT R»PCP0INT A.PCPOINT AA.CNPOINT P» 00012200?
PCPOINT B.PCPOINT X.PCPOINT BB.CiNT I) | 00012300?
BEqIN * SOLVES AX«B. WhERE A IS AN N*2N cOHPlEX UNSY"NETKIC ANq 00012400?
» b is an n*?r complex matrices* respectively. aa 1$ lu 00012500?
* decomposition proouceo by "cdecompose". the hesiuuals 00012600? % bb«b-ax are also calculated and ad*88 is solved. over- 00012700?
* Writing d on bb. l is the number of iterations. 00012800?
00012900?
SUBROUTINE CS0LVE(CINT N.CINT R.PCPOINT A.CNPOINT iNT. 00013000?
PCPOINT B)| 00013100?
BtGIN 00013200?
CINT I,J,K,KK,P,PPl OOOI33OO?
PREAL VECTOR BT[6«]I 00013400?
CREaL X.Y»Z.21.Z2»H.HH» 00013500?
PREAL VECTOR Ct6«l| 00013600?
BOOLEAN AMOOEI 00013700?
LouP I««1»1.N On 00013800?
IF INTCII NEQ I THEN 00013900?
LOOP Jt«R*R."l»l 00 00014000?
begin oooiaioo?
Xi«QRABONE(B£IJ»J-l)l 000U200?
MODEI-TRUE AND PEN"J-ll 00014300?
B[Ilt»GRABONE<BClNT[U].J-l)l 0001 4400?
Bt INTC 111 l"XI 00014500?
EnOI 00014600?
LOOP K»«R*R»-2.2 DO 00014700?
BEGIN 00014800?
KkIbK-11 00014900?
LOOP Il'l.l.N 00 00015000?
BEGIN 00015100?
MODEl'TRUE ANO PEN LSS 1*1-2 J 00015200?
lNNPROO(-GRABONE(B(I]*KK-l).0»AtI]»BT(KK)»H.HH)| 00015300?
lNNPR0DCH,-HH.RTL(l..AHl).BTtK].X.HH)l 00015400?
INNPR0D(-GRAB0NECBII1»K-1)»0#RTL(1»»A[I]).BTIKK]»H»HH)J 00015500?
INNPR00(H.HH»A[I J.BTtKl.H.HH)) 00015600?
Y|«-HI 00015700?
P|«I«II PP|«P-i| 00015800?
Zi!»GRABONE( ACIJ.PP-DI 00015900?
Z2»»GRABoNF(AU J.p-1 >1 00016000?
ZiaZl*Zl«Z?*42l 00016100?
MoOEl«TRUE AND PEN.KKJ 00016200?
Btl Jl.(X#Zl*Y*Z2)/Zl 00016300?
M00EI"TRUE AND PENnKI 00016400?
Bf I]lB(V*Zl-X*Z2)/Zf 00016500?
END; 00016600?
LOOp Ii«N.-l,l DO 00016700?
BEGIN 00016800?
AMODE«« TRUE AND PEN GEO NJ 00016900?
M0DE«"AM0DEJ 00017000?
lNNpRGD<-GRABoNE<Btl)»KK-l).0»Atn»BT[KK]»H,HH)i 000171 00?
32
lNNPROD(-H,-HH.RTL(l»»A(l])»BTtK)»H,HH)l 00017200?
MoOEl»TRUE ANO PEN»KK-il C C I } I «HI 00017300?
MoOEl«AMOO£| 00017*00?
lNNpR00(*6RAB0NE(BCn»K"l)»0«RTi.(l**ACI])*BTUK]»H«HH)l 00017500?
INNPROO(H,HH»ACI]»BT[KJ,H.HH)I 000l7ft00?
MOOEl«TRuE AND PEN«K"ll C[I]H-HI 00017700?
M0DE|«TRUE AND <PEN«KK-1 OR PEN»K-1)I B[I)|nC[I]l 00017800?
ENOl 00017900?
EnDj 00018000?
ENDi % OF CSOLVE, 00018100?
00018200?
C1NT I»J»K#D2.C»CC» 00018300?
CKEAL E*00*01iXHAX'BBHAX*EPS*El*E2»H»HH| OOOUaOO?
PHEAL VECTOR XTC64 J j 00018500?
BUOtEAN AMODEl 00018600?
LaB£L L3»lLLl 00018700?
EKS|sl(0P«10l 00018800?
AMODE|«pEN GEO N| 00018900?
MUOEUAmOOEI 00019000?
CUO* I»»1»1«N Dq 00019100?
BtLGlN XtI]l>O.OI 00019200?
BB< I 1 l"BC Hi 00019300?
ENDj 00019400?
L»«0» Doinfll 00019500?
MUDEI-TRUEl 00019600?
L3I CJ»OLVE(N#R.AAiP.BB)l 00019700?
L««L*1» D2t«0l D1I»0I 00019800?
MoDEUAMODE» 00019900?
LUQP H«1#1»N 00 00020000?
XtllUXtU + BBCllI 00020100?
LUOP jl.l.l.R 00 00020200?
BLGiN 00020300?
XMAXI.Ol BBMAXl«0l 00020400?
C»«J*J» CCl-C-ll 00020500?
LOOp Il«l#l.N 00 00020600?
BEGIN 00020700?
El I.GRAB0NE(X[IJ#CC-1)I E2 | ■GRABONEC X 1 1 ] # C-l ) I 00020800?
Ei"El*El*E2*E2l 00020900?
Ip E GTR XMAX THEN XMAX|»E| 00021000?
El «»GRftB0N£(BBtl J»CC)I t2 I "GR aBONE ( BB 1 1 ) » C ) 1 00021100?
Ei«El*El*E2*E2* 00021200?
IF E GTR BBMAX THEN BBHAXlaEI 00021300?
0002U00?
00021500?
00021600?
00021700?
MnOEl.TRUEi 00021800?
lNNpR0D(-GRAB0NE(B[I]»C)»0.RTL(l».ACl J)»XT[C]»M»HHJ| 00021900?
lNNpRODCH.HH.ACll.xT(cl»H»MH)l nn07?nnn«>
MnoE««TwuE and pen»c-u bbi u i«-hi
in .
IF BBHAx/XMAx C,T« 01 THEN 01 | .BBHAX/XM Ax I IF BBMAX GTR (2*EPS)*C2*EPS)*XHAX THEN D2l»H 101
Ip E GTR XMAX THEN XMAX|»E| 00021000?
El »»GRABONE(BBtl J»CC)I fc2 » -GR aBONE ( BB t I ) » C ) 1 00021100?
E»»El*El*E2*E2* 00021200?
IF E GTR BBMAX THEN BBHAXlaEJ
MnDE'"TRUEl
lNNpR0O(-GRABnNE(Btn.CC)»O,A[n.xTtCCl.H,HH)|
lNNpROOCH»-HH»RTL(l»»AC I 1 ) • XTt C 1 # H. HH) I
MODEI"TRUE AND PEN«CC"1» bBHll^HI MnDri«TRU£i
00022000?
00022100?
00022200?
00022300?
00022000? E^OJ 00022500?
It 01 GTR 00*0.25 AND L NEO 1 THEN GO TO ILLI 00022600?
OUi-Op 00022700?
IF qZ«1 THEN Go TO l3| 00022800?
33
ILLI 00022900?
ENDI X 0F CaCCSqLVc. 00023000?
00023100?
PREAL VECTOR AAiBB[6<»}» 00023200?
creal vector iNTcHt6*)i 00023300?
CINT I»N»R»OeTE«ITi 00023400?
CREAL DETR.OETII 00023500?
LOOP Ilal.l.N 00 00023600?
AAC I ] l»AC III 00023700?
COECoMP0Se(n»AA»OEtR»dEtI»OETE»InTCH)I 00023800?
CACCSOLvE(NtR»A*AA,lNTCH*B»X«BB*lT)i 00023900?
END| * OF CLlNSYS. 00024000?
00024100?
00024200?
SUBROUTINE INVITERATIONCCINT N»PCP0INT W.pCPOINT X.CINT KK)| 00024300?
BEGIN X ThE SUBROUTINE COMPUTES THE EIGENVECTORS AND PKlNClPAL 00024400?
X VECTORS BY SOLVING ( ( W-LAMBDA*! )**K )*X»0, THE SUBROUTINES 00024500?
I "cDECOMPOSE" AND "CACCSQLyE" IN "cLINSyS" ARE USED* 00024600?
PE hEaL Su8ROuTlNE sQRT As r6A(Pe rEAL ARQ AS r6A)I 00024700?
BEGIN 00024a00?
S* CALL SqRt64()| 00024900?
ENDI % QF SORT. 00025000?
CINT I»J»II»DETE.R»L» 00025100?
CKEAL DETRtDETI.MAXVAL.XXl,XX2| 00025200?
pkeal vector w2» b» b2» bb» xxt641 , 00025300?
ckeal vector p.rootcmji 00025400?
PHeAL Tri»TR2» 00025500?
LABEL INVIT»NEXT.FINISI 00025600?
MODEi«PEN LSS N*N| 00025700?
LOOP I>»1'1»N 00 W2C X 1 1 «•*£ III 00025800?
CDEC0MP0SE{N,W2.DETR#DETI#DETE.P)I 00025900?
M0DE1.PEN LSS 21 00026000?
LOOP ll«l*l#N DO 00026100?
BEGIN B2(I]lal,0l BC 1 1 l«B2C III ENDI 00026200?
KKl.ll 00026300?
Rl»ll 00026^00?
INVITI 00026500?
MOOE»»PEN LSS N*NI 00026600?
CACcS0LVE(N»R»m»m2»P»B*X*BB»L)I 00026700?
HODe««PeN LSS 21 00026800?
LOOP !I»1»1*N DO 00026900?
BEGIN 00027000?
RnOT[I)l"ROHSUH(XIIJ*X[Il)| 00027100?
ROOT t I ]I«SQRT(ROOTII ])l 00027200?
ENOi 00027300?
MAXvALI«ROOT[lll 00027400?
LOOP lt«2»l*N DO 00027500?
IT m*XVAL GEO HOOTtI] THEN MAXVAL I »MAXVAL ELSE MAXVal l»ROOT[ Ij I 00027600?
LOOP lialiUN DO 00027700?
IF ROOTtlj.MAXVAL THEN BEGIN IllaH GO TO NEXTl ENU| 00027800?
MEXTf 00027900?
XXl|»GRABONE(XUn»0)l 00028000?
xX2l«-6RAB0NE(XtII]»l)l 00028100?
IF -XX2 NEQ 0 THEN 00028200?
BEGIN 00028300?
LQOp I i * 1 • 1 • N DO 00028400?
BEGIN 00028500?
3^
MODEl«PEN LSS 21 00028600?
TRll»GRABnNE(Xtl].0)*XXl-GRABoNC(Xtl].l)*XX2| 00028700?
TR2l«GRAB0NE(X[n,0)*XX2*QRAB0NE(XCI].l)*XAil 00028800?
M00El«PEN«0l XXCIll.TRli 00028900?
M0DEl«PEN«ll XX(I}|»TR2| 00029000?
END| 00029100?
MQDEI.PEN LSS 21 00029200?
Xxl|«GRABONE(XXtII).0>| 00029300?
LOOP IIP1.1.N 00 XUJlOXCll/XXU 00029400?
END 00029*00?
ELSE LOOP Il.l.l»N 00 X[ IJ I.XC H/XXl | 00029600?
LOOP II-l.l.N DO 00029700?
ErlN -. «, . 00029800?
IP ABSCB2[I].X£I]> GTR P.lO THEN 00029900?
WE9IN 00030000?
IP KK GTR 10 THEN GO TO FlNISl 00030100?
loop ji.i.i.n oo begin Btnt.xtn» 00030200?
B2tI1l"X[I]l 00030300?
Kkl.KK*it END' 00030400?
*K»"KK*1I 00030500?
G0 TO INVITI 00030600?
If?r rn Tn r»m<;. 00030700?
LLSE GO TO FINISI 00030800?
„„,., H°l 00030900?
*** 0O031OOD?
CNOJ UF INVITERATION. 00031100?
UNCLASSIFIED
Security Classification
DOCUMENT CONTROL DATA • R & D
(SacuHty claaallleatlor, at till: body ol abahmet and humming annotation muat ba antarad whan tha ormr.lt raport t. claaalllad)
NAT.NG ACTIVITY fCo^Or.,. author) »*• ""OUT .ECU R. TV C L ASS. FIC A TION
originating activity (Corporata author)
Center for Advanced Computation University of Illinois at Urbana-Champaign Urbana, Illinois 6l801
UNCLASSIFIED
2b. GROUP
3. REPORT TITLE
A NUMERICAL SOLUTION OF THE MATRIX RICCATI EQUATIONS
4. DESCRIPTIVE NOTES (Type ot report and htelualra datam)
Research Report
8. AUTHOR(S) (Flrat mm, middla Initial, la a ( nam*)
Killion Noh
8. REPORT DATE
January 20, 1972
•a. CONTRACT OR GRANT NO.
DAHCOU 72-C-OOOl
b. PROJECT NO.
ARPA Order 1899
fm. TOTAL NO. OF PAGES
.22.
17b. NO. Or REFS
•a. ORIGINATOR'S REPORT NUMBERIS)
CAC Document No. 2k
Sb. OTHER REPORT NO<3» (Any other number* that may ba aemlgnad thla raport)
UIUCDCS-R-72-U98
10. DISTRIBUTION STATEMENT
Copies may "be obtained from the address given in (l) above. Distribution unlimited; approved for public release.
II. SUPPLEMENTARY NOTES
None
12. SPONSORING MILITARY ACTIVITY
U.S. Army Research Office-Durham Durham, North Carolina
13. ABSTRACT
The eigenvector solution of the time -invariant matrix Riccati equation is discussed. The coefficient matrix of the canonical equation is allowed to have multiple eigenvalues, namely, the matrix could be either derogatory or defective. The solution of matrix Riccati equation is then calculated from a part of similarity transformation which should reduce the coefficient matrix to the Jordan canonical form.
DD
FORM
,1473
UNCLASSIFIED
Security Classification
UNCLASSIFIED
Security Classification
KEY WORDS
LINK A
ROLE WT ROLE WT
ROLE WT
Ordinary and Partial Differential Equations
Matrix Algebra
Nonlinear and Functional Equations
UNCLASSIFIED
Security Classification
BIBLIOGRAPHIC DATA SHEET
1. Report No.
UIUCDCS-R-72-1+98
3. Recipient's Accession No.
4. Title and Subtitle
A NUMERICAL SOLUTION OF THE MATRIX RICCATI EQUATIONS
5. Report Date
January 20, 1972
6.
7. Author(s)
Killion Noh
8. Performing Organization Rept. No.
9. Performing Organization Name and Address
Center for Advanced Computation University of Illinois at Urbana-Champaign Urbana, Illinois 6l801
10. Project/Task/Work Unit No.
11. Contract /Grant No.
DAHCOU 72-C-OOOl
12. Sponsoring Organization Name and Address
U.S. Army Research Office-Durham
Duke Station
Durham, North Carolina
13. Type of Report & Period Covered
Research
14.
15. Supplementary Notes
None
16. Abstracts
The eigenvector solution of the time -invariant matrix Riccati equation is discussed. The coefficient matrix of the canonical equation is allowed to have multiple eigenvalues, namely, the matrix could be either derogatory or defective. The solution of matrix Riccati equation is then calculated from a part of similarity transformation which should reduce the coefficient matrix to the Jordan canonical form.
17. Key Words and Document Analysis. 17o. Descriptors
Ordinary and Partial Differential Equations
Matrix Algebra
Nonlinear and Functional Equations
17b. Identifiers/Open-Ended Terms
17c. COSATI Field/Group
18. Availability Statement Copies may be obtained f rom the
address in (9) above. Distribution unlimited.
19. Security Class (This Report)
UNCLASSIFIED
20. Security Class (This
Page UNCLASSIFIED
21. No. of Pages
38
22. Price
NC
FORM NTIS-35 (10-70)
USCOMM-DC 40329-P71