wscott@EN.ECN.PURDUE.EDU (Wayne H Scott) (03/14/90)
Posted for a friend:
From: Philip Peleties dominop@ec.ecn.purdue.edu
Solution of matrix equation
A X + X B = C
via Kronecker Product Method
Objective: Compute the matrix solution X of the matrix equation
A X + X B = C.
Method: Utilization of Kronecker products. Solution is:
vec (X) = INV ((TRANS (B) # I + I # A)) vec (C), where #
stands for Kronecker product and I is an appropriately
dimensioned identiry matrix. For Kronecker products
see: "Kronecker Products and Matrix Calculus in System
Theory", by John W. Brewer, IEEE Transactions on
Circuits and Systems, Vol. CAS-25, No. 9, September
1978. If you are REALLY interested in efficient solution
of A X + X B = C look at "Solution of the Matrix
Equation A X + X B = C", by R. H. Bartels and G. W.
Stewart, Communications of the ACM, September 1972,
Volume 15, Number 9. I warn you though that their
algorithm is long and not so easily implementable on a
programmable calculator.
Input: a) Square matrix A
b) Square matrix B
b) Square matrix C
Output: Square matrix C.
Notes: An example of a possible input looks like:
3: [[ -2 0 ]
[ 0 -1 ]]
2: [[ 3 0 ]
[ 0 4 ]]
1: [[ 1 0 ]
[ 0 1 ]]
then the output looks like:
1: [[ 1 0 ]
[ 0 .33333333333 ]]
Bugs: For dimensions in excess of 5, it takes a looong time to
get an answer ...
Also, if the eigenvalues of A are {l1,...,ln} and those
of B are {m1,...,mn} then TRANS (B) # I + I # A will be
non-invertible if and only if li + mj = 0 for every
i=1..n, j=1..n.
---------------------------------------------------------------------
Program AXXB
<< -> A B C
<< A SIZE LIST->
DROP DROP -> n
<< B TRN n IDN
KRON n IDN A KRON +
C VEC SWAP / A SIZE
DVEC
>>
>>
>>
---------------------------------------------------------------------
Program VEC
<< TRN ARRY-> LIST->
DROP * 1 2 ->LIST
->ARRY
>>
---------------------------------------------------------------------
Program DVEC
<< -> n
<< ARRY-> DROP n
->ARRY TRN
>>
>>
---------------------------------------------------------------------
Program KRON
<< -> A B
<< A SIZE LIST->
DROP DROP -> n
<< n SQ DUP 2
->LIST 0 CON 1 n
FOR i 1 n
FOR j A i j
2 ->LIST GET B * -> aa
<< i 1 - n
* 1 + i n *
FOR k j
1 - n * 1 + j n *
FOR l
aa k i 1 - n * - l j
1 - n * - 2 ->LIST
GET k l 2 ->LIST SWAP
PUT
NEXT
NEXT
>>
NEXT
NEXT
>>
>>
>>
------------------------------------------------------------------------
Philip Peleties dominop@ec.ecn.purdue.edu
_______________________________________________________________________________
Wayne Scott | INTERNET: wscott@en.ecn.purdue.edu
Electrical Engineering | BITNET: wscott%ea.ecn.purdue.edu@purccvm
Purdue University | UUCP: {purdue, pur-ee}!en.ecn.purdue.edu!wscott