[66] in mathematical software users group
MAPLE question- inveritng complex symbolic matirix
daemon@ATHENA.MIT.EDU (Saurav Dev Bhatta)
Sat Jan 23 16:40:46 1993
Date: Sat, 23 Jan 93 16:40:03 EST
From: bhattas@allegro.mit.edu (Saurav Dev Bhatta)
To: msug@MIT.EDU
Hi:
I have been trying to invert a 9 by 9 complex symbolic matrix using MAPLE. The
problem I have been facing is that Maple crashes while doing this.
My problem is as follows:
PROBLEM:
I want to invert a 9 by 9 matrix R. Rather than try to invert R by
taking the whole matrix at once, it is usually easier to break it into
blocks and try to invert it as follows:
|B C| |X Y|
R = | | Rinv = R^(-1)=| |
|D E| |Z U|
where:
X = (B - C E^-1 D)^-1
U = (E - D B^-1 C)^-1
Y = - B^-1 C U
Z = - E^-1 D X
I have written a program called "sec_oder" to first create the matrix
R (R has a very special form; hence sec_oder needs to call another
procedure called "expectation" when creating the matrix R). From what
I understand, MAPLE ran into problems when computing the inversion
(B - C E^-1 D)^-1, which, in the program, is done by the statement
X := inverse(temp2). I have included the listings of the two
procedures below.
QUESTION:
If there is some way to make this procedure work, or if there is some
other more effiecient way to do this inversion, I would be very glad
to hear of it (R is a complex hermetian 9 by 9 matrix). Thanks a lot
for your time.
Saurav Dev Bhattas
(bhattas@allegro.mit.edu)
***************************************************************
Address:
Home: 88 Beacon Street, Apt 54 Off: Room 36-315
Somerville, MA 02143 MIT
Ph: 617-547-6426 Ph: 617-253-7322
***************************************************************
#**************************************************************
#need to do with(linalg) first;
sec_odr := proc(z)
local i,ii,j,k,l,m,ww0,ww,tildw,conjtildw,R,temp, B,C,D,E,X,U,Y,Z,Binv,Einv,temp1,temp2,tempinv;
ww := arr(1..8);
ww0 := z[i,j];
ww[1] := z[i,j-1];
ww[3] := z[i-1,j];
ww[5] := z[i,j+1];
ww[7] := z[i+1,j];
ww[2] := z[i-1,j-1];
ww[4] := z[i-1,j+1];
ww[6] := z[i+1,j+1];
ww[8] := z[i+1,j-1];
tildw := matrix(1,1,[ww0]);
# create a 9 element long vector called tildw whose elements are
# ww0,..,ww[8].
for ii by 1 to 8 do
tildw := extend(tildw,1,0,ww[ii]);
od;
conjtildw := map(conjugate,tildw);
# create a 9 by 9 matrix called temp by doing mulstiplying
# tildw by its conjugate transpose
temp := multiply(tildw,transpose(conjtildw));
# create the matrix R by applying an operator called expectation
# on the matrix temp
R := map(expectation,temp);
# extract the blocks B, C, D and E from R.
B := submatrix(R,1..5,1..5);
C := submatrix(R,1..5,6..9);
D := submatrix(R,6..9,1..5);
E := submatrix(R,6..9,6..9);
Einv := inverse(E);
# print(`Einv=`,Einv);
temp1 := evalm(C &* Einv &* D);
# print(`CEinvD= `,temp1);
temp2 := evalm(B-temp1);
# print(`E-CEinvD= `,temp2);
X := inverse(temp2);
Binv := inverse(B);
temp1 := evalm(D &* Binv &* C);
temp2 := evalm(E - temp1);
U := inverse(temp2);
Y := evalm(-Binv &* C &* U);
Z := evalm(-Einv &* D &* X);
writeto(foo);
print(`X= `,X);
print('Y= `,Y);
print(`Z= `,Z);
print(`U= `,U);
writeto(terminal);
end;
#************************************************
`type/expectationargs` := [indexed &* function(indexed)];
expectation := proc(expr)
local i, j, k, l, rh, rv,rhoh,rhov,rhoh_re,rhoh_im,rhov_re,rhov_im;
if not type([expr], expectationargs) then
ERROR(
`invalid arguments; must be in form of indexed * function(indexed)`)
fi;
i := op(1,op(1,expr));
j := op(2,op(1,expr));
k := op(1,op(1,op(2,expr)));
l := op(2,op(1,op(2,expr)));
if (j-l) >= 0 then
rh := rhoh^(j-l);
else
rh := conjugate(rhoh^(l-j));
fi;
if (i-k) >= 0 then
rv := rhov^(i-k);
else
rv := conjugate(rhov^(k-i));
fi;
ans := rh*rv;
end;
