Artifact 16e9aeb960b2653473e37d3b22e7450ac0a4d51f7e28e4edde2911f017cf760f:
- File
r34.1/xlog/matrix.log
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 50119) [annotate] [blame] [check-ins using] [more...]
Sat May 30 16:11:48 PDT 1992 REDUCE 3.4.1, 15-Jul-92 ... 1: 1: 2: 2: 3: 3: Time: 17 ms 4: 4: % Tests of eigenfunction/eigenvalue code. v := mat((1,1,-1,1,0),(1,2,-1,0,1),(-1,2,3,-1,0), (1,-2,1,2,-1),(2,1,-1,3,0))$ mateigen(v,et); 4 3 2 {{ET - 6*ET + 13*ET + 5*ET - 5, 1, [ 5*ARBCOMPLEX(1)*(ET - 2) ] [ ---------------------------- ] [ 3 2 ] [ 2*ET - 10*ET + 23*ET + 5 ] [ ] [ 2 ] [ ARBCOMPLEX(1)*ET*( - ET + 6*ET - 8) ] [--------------------------------------] [ 3 2 ] [ 2*ET - 10*ET + 23*ET + 5 ] [ ] [ ARBCOMPLEX(1)*ET*( - 3*ET + 7) ] [ -------------------------------- ] [ 3 2 ] [ 2*ET - 10*ET + 23*ET + 5 ] [ ] [ 3 2 ] [ ARBCOMPLEX(1)*(ET - 4*ET + 10) ] [ ---------------------------------- ] [ 3 2 ] [ 2*ET - 10*ET + 23*ET + 5 ] [ ] [ ARBCOMPLEX(1) ] }, {ET - 2, 1, [ 0 ] [ ] [ 0 ] [ ] [ARBCOMPLEX(2)] [ ] [ARBCOMPLEX(2)] [ ] [ARBCOMPLEX(2)] }} eigv := third first ws$ % Now check if the equation for the eigenvectors is fulfilled. Note % that also the last component is zero due to the eigenvalue equation. v*eigv-et*eigv; [ 0 ] [ ] [ 4 3 2 ] [ ARBCOMPLEX(1)*(ET - 6*ET + 13*ET + 5*ET - 5) ] [ ------------------------------------------------- ] [ 3 2 ] [ 2*ET - 10*ET + 23*ET + 5 ] [ ] [ 0 ] [ ] [ 4 3 2 ] [ ARBCOMPLEX(1)*( - ET + 6*ET - 13*ET - 5*ET + 5) ] [ ---------------------------------------------------- ] [ 3 2 ] [ 2*ET - 10*ET + 23*ET + 5 ] [ ] [ 4 3 2 ] [ 2*ARBCOMPLEX(1)*( - ET + 6*ET - 13*ET - 5*ET + 5) ] [------------------------------------------------------] [ 3 2 ] [ 2*ET - 10*ET + 23*ET + 5 ] % Example of degenerate eigenvalues. u := mat((2,-1,1),(0,1,1),(-1,1,1))$ mateigen(u,eta); {{ETA - 1,2, [ARBCOMPLEX(3)] [ ] [ARBCOMPLEX(3)] [ ] [ 0 ] }, {ETA - 2,1, [ 0 ] [ ] [ARBCOMPLEX(4)] [ ] [ARBCOMPLEX(4)] }} % Example of a fourfold degenerate eigenvalue with two corresponding % eigenvectors. w := mat((1,-1,1,-1),(-3,3,-5,4),(8,-4,3,-4), (15,-10,11,-11))$ mateigen(w,al); {{AL + 1, 4, [ ARBCOMPLEX(5) ] [ --------------- ] [ 5 ] [ ] [ - 5*ARBCOMPLEX(6) + 7*ARBCOMPLEX(5) ] [--------------------------------------] [ 5 ] [ ] [ ARBCOMPLEX(5) ] [ ] [ ARBCOMPLEX(6) ] }} eigw := third first ws; [ ARBCOMPLEX(5) ] [ --------------- ] [ 5 ] [ ] [ - 5*ARBCOMPLEX(6) + 7*ARBCOMPLEX(5) ] EIGW := [--------------------------------------] [ 5 ] [ ] [ ARBCOMPLEX(5) ] [ ] [ ARBCOMPLEX(6) ] w*eigw - al*eigw; - ARBCOMPLEX(5)*(AL + 1) MAT((---------------------------), 5 ((5*ARBCOMPLEX(6)*AL + 5*ARBCOMPLEX(6) - 7*ARBCOMPLEX(5)*AL - 7*ARBCOMPLEX(5))/5), ( - ARBCOMPLEX(5)*(AL + 1)), ( - ARBCOMPLEX(6)*(AL + 1))) % Calculate the eigenvectors and eigenvalue equation. f := mat((0,ex,ey,ez),(-ex,0,bz,-by),(-ey,-bz,0,bx), (-ez,by,-bx,0))$ factor om; mateigen(f,om); 4 2 2 2 2 2 2 2 2 2 {{OM + OM *(EX + EY + EZ + BZ + BY + BX ) + EX *BX 2 2 + 2*EX*EY*BY*BX + 2*EX*EZ*BZ*BX + EY *BY + 2*EY*EZ*BZ*BY 2 2 + EZ *BZ , 1, 2 MAT(((OM *ARBCOMPLEX(7)*EZ + OM*ARBCOMPLEX(7)*( - EX*BY + EY*BX) 3 + ARBCOMPLEX(7)*BZ*(EX*BX + EY*BY + EZ*BZ))/(OM 2 2 2 + OM*(EX + EY + BZ ))), 2 (( - OM *ARBCOMPLEX(7)*BY + OM*ARBCOMPLEX(7)*( - EX*EZ + BZ*BX) 3 - ARBCOMPLEX(7)*EY*(EX*BX + EY*BY + EZ*BZ))/(OM 2 2 2 + OM*(EX + EY + BZ ))), 2 ((OM *ARBCOMPLEX(7)*BX + OM*ARBCOMPLEX(7)*( - EY*EZ + BZ*BY) 3 + ARBCOMPLEX(7)*EX*(EX*BX + EY*BY + EZ*BZ))/(OM 2 2 2 + OM*(EX + EY + BZ ))), (ARBCOMPLEX(7))) }} % Specialize to perpendicular electric and magnetic field. let ez=0,ex=0,by=0; % Note that we find two eigenvectors to the double eigenvalue 0 % (as it must be). mateigen(f,om); {{OM, 2, [ ARBCOMPLEX(9)*BX - ARBCOMPLEX(8)*BZ ] [-------------------------------------] [ EY ] [ ] [ ARBCOMPLEX(8) ] [ ] [ 0 ] [ ] [ ARBCOMPLEX(9) ] }, 2 2 2 2 {OM + EY + BZ + BX , 1, [ - ARBCOMPLEX(10)*EY ] [ ---------------------- ] [ BX ] [ ] [ - ARBCOMPLEX(10)*BZ ] [ ---------------------- ] [ BX ] [ ] [ 2 2 2 ] [ ARBCOMPLEX(10)*(EY + BZ + BX ) ] [----------------------------------] [ OM*BX ] [ ] [ ARBCOMPLEX(10) ] }} % The following has 1 as a double eigenvalue. The corresponding % eigenvector must involve two arbitrary constants. j := mat((9/8,1/4,-sqrt(3)/8), (1/4,3/2,-sqrt(3)/4), (-sqrt(3)/8,-sqrt(3)/4,11/8)); [ 9 1 - SQRT(3) ] [ --- --- ------------] [ 8 4 8 ] [ ] [ 1 3 - SQRT(3) ] J := [ --- --- ------------] [ 4 2 4 ] [ ] [ - SQRT(3) - SQRT(3) 11 ] [------------ ------------ ---- ] [ 8 4 8 ] mateigen(j,x); {{X - 1, 2, [SQRT(3)*ARBCOMPLEX(12) - 2*ARBCOMPLEX(11)] [ ] [ ARBCOMPLEX(11) ] [ ] [ ARBCOMPLEX(12) ] }, {X - 2, 1, [ - SQRT(3)*ARBCOMPLEX(13) ] [ --------------------------- ] [ 3 ] [ ] [ - 2*SQRT(3)*ARBCOMPLEX(13) ] [-----------------------------] [ 3 ] [ ] [ ARBCOMPLEX(13) ] }} % The following is a good consistency check. sym := mat( (0, 1/2, 1/(2*sqrt(2)), 0, 0), (1/2, 0, 1/(2*sqrt(2)), 0, 0), (1/(2*sqrt(2)), 1/(2*sqrt(2)), 0, 1/2, 1/2), (0, 0, 1/2, 0, 0), (0, 0, 1/2, 0, 0))$ ans := mateigen(sym,eta); ANS := {{ETA, 1, [ 0 ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ - ARBCOMPLEX(14)] [ ] [ ARBCOMPLEX(14) ] }, {ETA - 1, 1, [ 2*ARBCOMPLEX(15) ] [------------------] [ SQRT(2) ] [ ] [ 2*ARBCOMPLEX(15) ] [------------------] [ SQRT(2) ] [ ] [ 2*ARBCOMPLEX(15) ] [ ] [ ARBCOMPLEX(15) ] [ ] [ ARBCOMPLEX(15) ] }, {2*ETA + 1, 1, [ - ARBCOMPLEX(16)] [ ] [ ARBCOMPLEX(16) ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ 0 ] }, 2 {4*ETA + 2*ETA - 1, 1, [ - ARBCOMPLEX(17) ] [ ------------------- ] [ 2*SQRT(2)*ETA ] [ ] [ - ARBCOMPLEX(17) ] [ ------------------- ] [ 2*SQRT(2)*ETA ] [ ] [ ARBCOMPLEX(17)*( - 2*ETA + 1) ] [-------------------------------] [ 2*ETA ] [ ] [ ARBCOMPLEX(17) ] [ ] [ ARBCOMPLEX(17) ] }} % Check of correctness for this example. for each j in ans do for each k in solve(first j,eta) do write sub(k,sym*third j - eta*third j); [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] % Tests of nullspace operator. a1 := mat((1,2,3,4),(5,6,7,8)); [1 2 3 4] A1 := [ ] [5 6 7 8] nullspace a1; { [ 1 ] [ ] [ 0 ] [ ] [ - 3] [ ] [ 2 ] , [ 0 ] [ ] [ 1 ] [ ] [ - 2] [ ] [ 1 ] } b1 := {{1,2,3,4},{5,6,7,8}}; B1 := {{1,2,3,4},{5,6,7,8}} nullspace b1; {{1,0,-3,2},{0,1,-2,1}} % Example taken from a bug report for another CA system. c1 := {{(p1**2*(p1**2 + p2**2 + p3**2 - s*z - z**2))/(p1**2 + p3**2), 0, (p1*p3*(p1**2 + p2**2 + p3**2 - s*z - z**2))/(p1**2 + p3**2), -((p1**2*p2*(s + z))/(p1**2 + p3**2)), p1*(s + z), -((p1*p2*p3*(s + z))/(p1**2 + p3**2)), -((p1*p3*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)), 0, (p1**2*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)}, {0, 0, 0, 0, 0, 0, 0, 0, 0}, {(p1*p3*(p1**2 + p2**2 + p3**2 - s*z - z**2))/(p1**2 + p3**2), 0, (p3**2*(p1**2 + p2**2 + p3**2 - s*z - z**2))/(p1**2 + p3**2), -((p1*p2*p3*(s + z))/(p1**2 + p3**2)), p3*(s + z), -((p2*p3**2*(s + z))/(p1**2 + p3**2)), -((p3**2*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)), 0, (p1*p3*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)}, {-((p1**2*p2*(s + z))/(p1**2 + p3**2)), 0, -((p1*p2*p3*(s + z))/(p1**2 + p3**2)), -((p1**2*p2**2*(s + 2*z))/((p1**2 + p3**2)*z)), (p1*p2*(s + 2*z))/z, -((p1*p2**2*p3*(s + 2*z))/((p1**2 + p3**2)*z)), -((p1*p2*p3*z)/(p1**2 + p3**2)), 0, (p1**2*p2*z)/(p1**2 + p3**2)}, {p1*(s + z), 0, p3*(s + z), (p1*p2*(s + 2*z))/z, -(((p1**2+p3**2)*(s+ 2*z))/z), (p2*p3*(s + 2*z))/z, p3*z,0, -(p1*z)}, {-((p1*p2*p3*(s + z))/(p1**2 + p3**2)), 0, -((p2*p3**2*(s + z))/(p1**2 + p3**2)), -((p1*p2**2*p3*(s + 2*z))/((p1**2 + p3**2)*z)), (p2*p3*(s + 2*z))/z, -((p2**2*p3**2*(s + 2*z))/((p1**2 + p3**2)*z)), -((p2*p3**2*z)/(p1**2 + p3**2)), 0, (p1*p2*p3*z)/(p1**2 + p3**2)}, {-((p1*p3*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)), 0, -((p3**2*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)), -((p1*p2*p3*z)/(p1**2 + p3**2)),p3*z,-((p2*p3**2*z)/(p1**2 + p3**2)), -((p3**2*(p1**2 + p2**2 + p3**2)*z)/((p1**2 + p3**2)*(s + z))), 0, (p1*p3*(p1**2 + p2**2 + p3**2)*z)/((p1**2 + p3**2)*(s + z))}, {0, 0, 0, 0, 0, 0, 0, 0, 0}, {(p1**2*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2), 0, (p1*p3*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2), (p1**2*p2*z)/(p1**2 + p3**2), -(p1*z), (p1*p2*p3*z)/(p1**2 + p3**2), (p1*p3*(p1**2 + p2**2 + p3**2)*z)/((p1**2 + p3**2)*(s + z)), 0, -((p1**2*(p1**2 + p2**2 + p3**2)*z)/((p1**2 + p3**2)*(s + z)))}}; 2 2 2 2 2 P1 *( - S*Z - Z + P1 + P2 + P3 ) C1 := {{-------------------------------------, 2 2 P1 + P3 0, 2 2 2 2 P1*P3*( - S*Z - Z + P1 + P2 + P3 ) ---------------------------------------, 2 2 P1 + P3 2 - P1 *P2*(S + Z) -------------------, 2 2 P1 + P3 P1*(S + Z), - P1*P2*P3*(S + Z) ---------------------, 2 2 P1 + P3 2 2 2 - P1*P3*(P1 + P2 + P3 ) ----------------------------, 2 2 P1 + P3 0, 2 2 2 2 P1 *(P1 + P2 + P3 ) -----------------------}, 2 2 P1 + P3 {0,0,0,0,0,0,0,0,0}, 2 2 2 2 P1*P3*( - S*Z - Z + P1 + P2 + P3 ) {---------------------------------------, 2 2 P1 + P3 0, 2 2 2 2 2 P3 *( - S*Z - Z + P1 + P2 + P3 ) -------------------------------------, 2 2 P1 + P3 - P1*P2*P3*(S + Z) ---------------------, 2 2 P1 + P3 P3*(S + Z), 2 - P2*P3 *(S + Z) -------------------, 2 2 P1 + P3 2 2 2 2 - P3 *(P1 + P2 + P3 ) --------------------------, 2 2 P1 + P3 0, 2 2 2 P1*P3*(P1 + P2 + P3 ) -------------------------}, 2 2 P1 + P3 2 - P1 *P2*(S + Z) {-------------------, 2 2 P1 + P3 0, - P1*P2*P3*(S + Z) ---------------------, 2 2 P1 + P3 2 2 P1 *P2 *( - S - 2*Z) ----------------------, 2 2 Z*(P1 + P3 ) P1*P2*(S + 2*Z) -----------------, Z 2 P1*P2 *P3*( - S - 2*Z) ------------------------, 2 2 Z*(P1 + P3 ) - Z*P1*P2*P3 ---------------, 2 2 P1 + P3 0, 2 Z*P1 *P2 -----------}, 2 2 P1 + P3 {P1*(S + Z), 0, P3*(S + Z), P1*P2*(S + 2*Z) -----------------, Z 2 2 2 2 - S*P1 - S*P3 - 2*Z*P1 - 2*Z*P3 --------------------------------------, Z P2*P3*(S + 2*Z) -----------------, Z Z*P3, 0, - Z*P1}, - P1*P2*P3*(S + Z) {---------------------, 2 2 P1 + P3 0, 2 - P2*P3 *(S + Z) -------------------, 2 2 P1 + P3 2 P1*P2 *P3*( - S - 2*Z) ------------------------, 2 2 Z*(P1 + P3 ) P2*P3*(S + 2*Z) -----------------, Z 2 2 P2 *P3 *( - S - 2*Z) ----------------------, 2 2 Z*(P1 + P3 ) 2 - Z*P2*P3 -------------, 2 2 P1 + P3 0, Z*P1*P2*P3 ------------}, 2 2 P1 + P3 2 2 2 - P1*P3*(P1 + P2 + P3 ) {----------------------------, 2 2 P1 + P3 0, 2 2 2 2 - P3 *(P1 + P2 + P3 ) --------------------------, 2 2 P1 + P3 - Z*P1*P2*P3 ---------------, 2 2 P1 + P3 Z*P3, 2 - Z*P2*P3 -------------, 2 2 P1 + P3 2 2 2 2 - Z*P3 *(P1 + P2 + P3 ) -------------------------------, 2 2 2 2 S*P1 + S*P3 + Z*P1 + Z*P3 0, 2 2 2 Z*P1*P3*(P1 + P2 + P3 ) -------------------------------}, 2 2 2 2 S*P1 + S*P3 + Z*P1 + Z*P3 {0,0,0,0,0,0,0,0,0}, 2 2 2 2 P1 *(P1 + P2 + P3 ) {-----------------------, 2 2 P1 + P3 0, 2 2 2 P1*P3*(P1 + P2 + P3 ) -------------------------, 2 2 P1 + P3 2 Z*P1 *P2 -----------, 2 2 P1 + P3 - Z*P1, Z*P1*P2*P3 ------------, 2 2 P1 + P3 2 2 2 Z*P1*P3*(P1 + P2 + P3 ) -------------------------------, 2 2 2 2 S*P1 + S*P3 + Z*P1 + Z*P3 0, 2 2 2 2 - Z*P1 *(P1 + P2 + P3 ) -------------------------------}} 2 2 2 2 S*P1 + S*P3 + Z*P1 + Z*P3 nullspace c1; - P1 {{1,0,-------,0,0,0,0,0,0}, P3 {0,1,0,0,0,0,0,0,0}, - P1 {0,0,0,1,0,-------,0,0,0}, P3 2 2 P1 + P3 {0,0,0,0,1,-----------,0,0,0}, P2*P3 P3 {0,0,0,0,0,0,1,0,----}, P1 {0,0,0,0,0,0,0,1,0}} d1 := mat (((p1**2*(p1**2 + p2**2 + p3**2 - s*z - z**2))/(p1**2 + p3**2), 0, (p1*p3*(p1**2 + p2**2 + p3**2 - s*z - z**2))/(p1**2 + p3**2), -((p1**2*p2*(s + z))/(p1**2 + p3**2)), p1*(s + z), -((p1*p2*p3*(s + z))/(p1**2 + p3**2)), -((p1*p3*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)), 0, (p1**2*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)), (0, 0, 0, 0, 0, 0, 0, 0, 0), ((p1*p3*(p1**2 + p2**2 + p3**2 - s*z - z**2))/(p1**2 + p3**2), 0, (p3**2*(p1**2 + p2**2 + p3**2 - s*z - z**2))/(p1**2 + p3**2), -((p1*p2*p3*(s + z))/(p1**2 + p3**2)), p3*(s + z), -((p2*p3**2*(s + z))/(p1**2 + p3**2)), -((p3**2*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)), 0, (p1*p3*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)), ( ((p1**2*p2*(s + z))/(p1**2 + p3**2)), 0, -((p1*p2*p3*(s + z))/(p1**2 + p3**2)), -((p1**2*p2**2*(s + 2*z))/((p1**2 + p3**2)*z)), (p1*p2*(s + 2*z))/z, -((p1*p2**2*p3*(s + 2*z))/((p1**2 + p3**2)*z)), -((p1*p2*p3*z)/(p1**2 + p3**2)), 0, (p1**2*p2*z)/(p1**2 + p3**2)), (p1*(s + z), 0, p3*(s + z), (p1*p2*(s + 2*z))/z, -(((p1**2 + p3**2)*(s + 2*z))/z),(p2*p3*(s + 2*z))/z,p3*z,0,-(p1*z)), (-((p1*p2*p3*(s + z))/(p1**2 + p3**2)), 0, -((p2*p3**2*(s + z))/(p1**2 + p3**2)), -((p1*p2**2*p3*(s + 2*z))/((p1**2 + p3**2)*z)), (p2*p3*(s + 2*z))/z, -((p2**2*p3**2*(s + 2*z))/((p1**2 + p3**2)*z)), -((p2*p3**2*z)/(p1**2 + p3**2)), 0, (p1*p2*p3*z)/(p1**2 + p3**2)), (-((p1*p3*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)), 0, -((p3**2*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2)), -((p1*p2*p3*z)/(p1**2 + p3**2)),p3*z,-((p2*p3**2*z)/(p1**2 + p3**2)), -((p3**2*(p1**2 + p2**2 + p3**2)*z)/((p1**2 + p3**2)*(s + z))), 0, (p1*p3*(p1**2 + p2**2 + p3**2)*z)/((p1**2 + p3**2)*(s + z))), (0, 0, 0, 0, 0, 0, 0, 0, 0), ((p1**2*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2), 0, (p1*p3*(p1**2 + p2**2 + p3**2))/(p1**2 + p3**2), (p1**2*p2*z)/(p1**2 + p3**2), -(p1*z), (p1*p2*p3*z)/(p1**2 + p3**2), (p1*p3*(p1**2 + p2**2 + p3**2)*z)/((p1**2 + p3**2)*(s + z)), 0, -((p1**2*(p1**2 + p2**2 + p3**2)*z)/((p1**2 + p3**2)*(s + z))))); 2 2 2 2 2 P1 *( - S*Z - Z + P1 + P2 + P3 ) D1 := MAT((-------------------------------------,0, 2 2 P1 + P3 2 2 2 2 P1*P3*( - S*Z - Z + P1 + P2 + P3 ) ---------------------------------------, 2 2 P1 + P3 2 - P1 *P2*(S + Z) - P1*P2*P3*(S + Z) -------------------,P1*(S + Z),---------------------, 2 2 2 2 P1 + P3 P1 + P3 2 2 2 2 2 2 2 - P1*P3*(P1 + P2 + P3 ) P1 *(P1 + P2 + P3 ) ----------------------------,0,-----------------------), 2 2 2 2 P1 + P3 P1 + P3 (0,0,0,0,0,0,0,0,0), 2 2 2 2 P1*P3*( - S*Z - Z + P1 + P2 + P3 ) (---------------------------------------,0, 2 2 P1 + P3 2 2 2 2 2 P3 *( - S*Z - Z + P1 + P2 + P3 ) -------------------------------------, 2 2 P1 + P3 2 - P1*P2*P3*(S + Z) - P2*P3 *(S + Z) ---------------------,P3*(S + Z),-------------------, 2 2 2 2 P1 + P3 P1 + P3 2 2 2 2 2 2 2 - P3 *(P1 + P2 + P3 ) P1*P3*(P1 + P2 + P3 ) --------------------------,0,-------------------------), 2 2 2 2 P1 + P3 P1 + P3 2 P1 *P2*(S + Z) - P1*P2*P3*(S + Z) (----------------,0,---------------------, 2 2 2 2 P1 + P3 P1 + P3 2 2 P1 *P2 *( - S - 2*Z) P1*P2*(S + 2*Z) ----------------------,-----------------, 2 2 Z Z*(P1 + P3 ) 2 2 P1*P2 *P3*( - S - 2*Z) - Z*P1*P2*P3 Z*P1 *P2 ------------------------,---------------,0,-----------), 2 2 2 2 2 2 Z*(P1 + P3 ) P1 + P3 P1 + P3 P1*P2*(S + 2*Z) (P1*(S + Z),0,P3*(S + Z),-----------------, Z 2 2 2 2 - S*P1 - S*P3 - 2*Z*P1 - 2*Z*P3 P2*P3*(S + 2*Z) --------------------------------------,-----------------, Z Z Z*P3,0, - Z*P1), 2 - P1*P2*P3*(S + Z) - P2*P3 *(S + Z) (---------------------,0,-------------------, 2 2 2 2 P1 + P3 P1 + P3 2 P1*P2 *P3*( - S - 2*Z) P2*P3*(S + 2*Z) ------------------------,-----------------, 2 2 Z Z*(P1 + P3 ) 2 2 2 P2 *P3 *( - S - 2*Z) - Z*P2*P3 Z*P1*P2*P3 ----------------------,-------------,0,------------), 2 2 2 2 2 2 Z*(P1 + P3 ) P1 + P3 P1 + P3 2 2 2 2 2 2 2 - P1*P3*(P1 + P2 + P3 ) - P3 *(P1 + P2 + P3 ) (----------------------------,0,--------------------------, 2 2 2 2 P1 + P3 P1 + P3 2 - Z*P1*P2*P3 - Z*P2*P3 ---------------,Z*P3,-------------, 2 2 2 2 P1 + P3 P1 + P3 2 2 2 2 - Z*P3 *(P1 + P2 + P3 ) -------------------------------,0, 2 2 2 2 S*P1 + S*P3 + Z*P1 + Z*P3 2 2 2 Z*P1*P3*(P1 + P2 + P3 ) -------------------------------), 2 2 2 2 S*P1 + S*P3 + Z*P1 + Z*P3 (0,0,0,0,0,0,0,0,0), 2 2 2 2 2 2 2 P1 *(P1 + P2 + P3 ) P1*P3*(P1 + P2 + P3 ) (-----------------------,0,-------------------------, 2 2 2 2 P1 + P3 P1 + P3 2 Z*P1 *P2 Z*P1*P2*P3 -----------, - Z*P1,------------, 2 2 2 2 P1 + P3 P1 + P3 2 2 2 Z*P1*P3*(P1 + P2 + P3 ) -------------------------------,0, 2 2 2 2 S*P1 + S*P3 + Z*P1 + Z*P3 2 2 2 2 - Z*P1 *(P1 + P2 + P3 ) -------------------------------)) 2 2 2 2 S*P1 + S*P3 + Z*P1 + Z*P3 nullspace d1; { [0] [ ] [1] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] , [ 0 ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ 1 ] [ ] [ 0 ] [ ] [ - P1 ] [-------] [ P3 ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ 0 ] , [ 0 ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ 1 ] [ ] [ 2 2 ] [ P1 + P3 ] [-----------] [ P2*P3 ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ 0 ] , [ 0 ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ 0 ] [ ] [ 1 ] [ ] [ 0 ] [ ] [ P3 ] [----] [ P1 ] , [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [1] [ ] [0] } % The following example, by Kenton Yee, was discussed extensively by % the sci.math.symbolic newsgroup. m := mat((e^(-1), e^(-1), e^(-1), e^(-1), e^(-1), e^(-1), e^(-1), 0), (1, 1, 1, 1, 1, 1, 0, 1),(1, 1, 1, 1, 1, 0, 1, 1), (1, 1, 1, 1, 0, 1, 1, 1),(1, 1, 1, 0, 1, 1, 1, 1), (1, 1, 0, 1, 1, 1, 1, 1),(1, 0, 1, 1, 1, 1, 1, 1), (0, e, e, e, e, e, e, e)); [ 1 1 1 1 1 1 1 ] [--- --- --- --- --- --- --- 0] [ E E E E E E E ] [ ] [ 1 1 1 1 1 1 0 1] [ ] [ 1 1 1 1 1 0 1 1] [ ] M := [ 1 1 1 1 0 1 1 1] [ ] [ 1 1 1 0 1 1 1 1] [ ] [ 1 1 0 1 1 1 1 1] [ ] [ 1 0 1 1 1 1 1 1] [ ] [ 0 E E E E E E E] eig := mateigen(m,x); EIG := {{X - 1, 3, [ 0 ] [ ] [ - ARBCOMPLEX(20)] [ ] [ - ARBCOMPLEX(19)] [ ] [ - ARBCOMPLEX(18)] [ ] [ ARBCOMPLEX(18) ] [ ] [ ARBCOMPLEX(19) ] [ ] [ ARBCOMPLEX(20) ] [ ] [ 0 ] }, {X + 1, 3, ARBCOMPLEX(23) MAT((----------------), E (ARBCOMPLEX(22)), (ARBCOMPLEX(21)), (( - ARBCOMPLEX(23)*E - ARBCOMPLEX(23) - 2*ARBCOMPLEX(22)*E - 2*ARBCOMPLEX(21)*E)/(2*E)), (( - ARBCOMPLEX(23)*E - ARBCOMPLEX(23) - 2*ARBCOMPLEX(22)*E - 2*ARBCOMPLEX(21)*E)/(2*E)), (ARBCOMPLEX(21)), (ARBCOMPLEX(22)), (ARBCOMPLEX(23))) }, 2 2 { - E *X + E*X - 6*E*X + 7*E - X, 1, 8 7 7 6 MAT(((6*ARBCOMPLEX(24)*(E *X + 23*E *X - 7*E + 179*E *X 6 5 5 4 4 - 119*E + 565*E *X - 581*E + 768*E *X - 890*E 3 3 2 2 + 565*E *X - 581*E + 179*E *X - 119*E + 23*E*X 3 8 7 7 6 - 7*E + X))/(E *(E *X + 30*E *X - 7*E + 333*E *X 6 5 5 4 - 168*E + 1692*E *X - 1365*E + 4023*E *X 4 3 3 2 - 4368*E + 4470*E *X - 5145*E + 2663*E *X 2 - 2520*E + 576*E*X - 251*E + 36*X))), 9 8 8 7 ((ARBCOMPLEX(24)*(E *X + 29*E *X - 7*E + 310*E *X 7 6 6 5 - 161*E + 1520*E *X - 1246*E + 3577*E *X 5 4 4 3 - 3836*E + 4283*E *X - 4795*E + 2988*E *X 3 2 2 - 3065*E + 978*E *X - 672*E + 132*E*X - 42*E 2 8 7 7 6 + 6*X))/(E *(E *X + 30*E *X - 7*E + 333*E *X 6 5 5 4 - 168*E + 1692*E *X - 1365*E + 4023*E *X 4 3 3 2 - 4368*E + 4470*E *X - 5145*E + 2663*E *X 2 - 2520*E + 576*E*X - 251*E + 36*X))), 9 8 8 7 ((ARBCOMPLEX(24)*(E *X + 29*E *X - 7*E + 310*E *X 7 6 6 5 - 161*E + 1520*E *X - 1246*E + 3577*E *X 5 4 4 3 - 3836*E + 4283*E *X - 4795*E + 2988*E *X 3 2 2 - 3065*E + 978*E *X - 672*E + 132*E*X - 42*E 2 8 7 7 6 + 6*X))/(E *(E *X + 30*E *X - 7*E + 333*E *X 6 5 5 4 - 168*E + 1692*E *X - 1365*E + 4023*E *X 4 3 3 2 - 4368*E + 4470*E *X - 5145*E + 2663*E *X 2 - 2520*E + 576*E*X - 251*E + 36*X))), 9 8 8 7 ((ARBCOMPLEX(24)*(E *X + 29*E *X - 7*E + 310*E *X 7 6 6 5 - 161*E + 1520*E *X - 1246*E + 3577*E *X 5 4 4 3 - 3836*E + 4283*E *X - 4795*E + 2988*E *X 3 2 2 - 3065*E + 978*E *X - 672*E + 132*E*X - 42*E 2 8 7 7 6 + 6*X))/(E *(E *X + 30*E *X - 7*E + 333*E *X 6 5 5 4 - 168*E + 1692*E *X - 1365*E + 4023*E *X 4 3 3 2 - 4368*E + 4470*E *X - 5145*E + 2663*E *X 2 - 2520*E + 576*E*X - 251*E + 36*X))), 9 8 8 7 ((ARBCOMPLEX(24)*(E *X + 29*E *X - 7*E + 310*E *X 7 6 6 5 - 161*E + 1520*E *X - 1246*E + 3577*E *X 5 4 4 3 - 3836*E + 4283*E *X - 4795*E + 2988*E *X 3 2 2 - 3065*E + 978*E *X - 672*E + 132*E*X - 42*E 2 8 7 7 6 + 6*X))/(E *(E *X + 30*E *X - 7*E + 333*E *X 6 5 5 4 - 168*E + 1692*E *X - 1365*E + 4023*E *X 4 3 3 2 - 4368*E + 4470*E *X - 5145*E + 2663*E *X 2 - 2520*E + 576*E*X - 251*E + 36*X))), 9 8 8 7 ((ARBCOMPLEX(24)*(E *X + 29*E *X - 7*E + 310*E *X 7 6 6 5 - 161*E + 1520*E *X - 1246*E + 3577*E *X 5 4 4 3 - 3836*E + 4283*E *X - 4795*E + 2988*E *X 3 2 2 - 3065*E + 978*E *X - 672*E + 132*E*X - 42*E 2 8 7 7 6 + 6*X))/(E *(E *X + 30*E *X - 7*E + 333*E *X 6 5 5 4 - 168*E + 1692*E *X - 1365*E + 4023*E *X 4 3 3 2 - 4368*E + 4470*E *X - 5145*E + 2663*E *X 2 - 2520*E + 576*E*X - 251*E + 36*X))), 9 8 8 7 ((ARBCOMPLEX(24)*(E *X + 29*E *X - 7*E + 310*E *X 7 6 6 5 - 161*E + 1520*E *X - 1246*E + 3577*E *X 5 4 4 3 - 3836*E + 4283*E *X - 4795*E + 2988*E *X 3 2 2 - 3065*E + 978*E *X - 672*E + 132*E*X - 42*E 2 8 7 7 6 + 6*X))/(E *(E *X + 30*E *X - 7*E + 333*E *X 6 5 5 4 - 168*E + 1692*E *X - 1365*E + 4023*E *X 4 3 3 2 - 4368*E + 4470*E *X - 5145*E + 2663*E *X 2 - 2520*E + 576*E*X - 251*E + 36*X))), (ARBCOMPLEX(24))) }} % Now check the eigenvectors and calculate the eigenvalues in the % respective eigenspaces: factor expt; for each eispace in eig do begin scalar eivaleq,eival,eivec; eival := solve(first eispace,x); for each soln in eival do <<eival := rhs soln; eivec := third eispace; eivec := sub(soln,eivec); write "eigenvalue = ", eival; write "check of eigen equation: ", m*eivec - eival*eivec>> end; eigenvalue = 1 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] eigenvalue = -1 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] 4 3 2 2 SQRT(E + 12*E + 10*E + 12*E + 1) + E + 6*E + 1 eigenvalue = ---------------------------------------------------- 2*E check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] 4 3 2 2 - SQRT(E + 12*E + 10*E + 12*E + 1) + E + 6*E + 1 eigenvalue = ------------------------------------------------------- 2*E check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] % For the special choice: let e = -7 + sqrt 48; % we get only 7 eigenvectors. eig := mateigen(m,x); EIG := {{X + 1, 4, ARBCOMPLEX(27) MAT((----------------), 4*SQRT(3) - 7 (ARBCOMPLEX(26)), (ARBCOMPLEX(25)), ((2*SQRT(3)*( - ARBCOMPLEX(27) - 2*ARBCOMPLEX(26) - 2*ARBCOMPLEX(25)) + 3*ARBCOMPLEX(27) + 7*ARBCOMPLEX(26) + 7*ARBCOMPLEX(25))/(4*SQRT(3) - 7 )), ((2*SQRT(3)*( - ARBCOMPLEX(27) - 2*ARBCOMPLEX(26) - 2*ARBCOMPLEX(25)) + 3*ARBCOMPLEX(27) + 7*ARBCOMPLEX(26) + 7*ARBCOMPLEX(25))/(4*SQRT(3) - 7 )), (ARBCOMPLEX(25)), (ARBCOMPLEX(26)), (ARBCOMPLEX(27))) }, {X - 1, 3, [ 0 ] [ ] [ - ARBCOMPLEX(30)] [ ] [ - ARBCOMPLEX(29)] [ ] [ - ARBCOMPLEX(28)] [ ] [ ARBCOMPLEX(28) ] [ ] [ ARBCOMPLEX(29) ] [ ] [ ARBCOMPLEX(30) ] [ ] [ 0 ] }, {X + 7, 1, [ ARBCOMPLEX(31) ] [ ----------------- ] [ 56*SQRT(3) - 97 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(31) + 24*ARBCOMPLEX(31) ] [--------------------------------------------------] [ 168*SQRT(3) - 291 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(31) + 24*ARBCOMPLEX(31) ] [--------------------------------------------------] [ 168*SQRT(3) - 291 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(31) + 24*ARBCOMPLEX(31) ] [--------------------------------------------------] [ 168*SQRT(3) - 291 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(31) + 24*ARBCOMPLEX(31) ] [--------------------------------------------------] [ 168*SQRT(3) - 291 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(31) + 24*ARBCOMPLEX(31) ] [--------------------------------------------------] [ 168*SQRT(3) - 291 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(31) + 24*ARBCOMPLEX(31) ] [--------------------------------------------------] [ 168*SQRT(3) - 291 ] [ ] [ ARBCOMPLEX(31) ] }} for each eispace in eig do begin scalar eivaleq,eival,eivec; eival := solve(first eispace,x); for each soln in eival do <<eival := rhs soln; eivec := third eispace; eivec := sub(soln,eivec); write "eigenvalue = ", eival; write "check of eigen equation: ", m*eivec - eival*eivec>> end; eigenvalue = -1 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] eigenvalue = 1 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] eigenvalue = -7 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] % The same behaviour for this choice of e. clear e; let e = -7 - sqrt 48; % we get only 7 eigenvectors. eig := mateigen(m,x); EIG := {{X + 1, 4, - ARBCOMPLEX(34) MAT((-------------------), 4*SQRT(3) + 7 (ARBCOMPLEX(33)), (ARBCOMPLEX(32)), ((2*SQRT(3)*( - ARBCOMPLEX(34) - 2*ARBCOMPLEX(33) - 2*ARBCOMPLEX(32)) - 3*ARBCOMPLEX(34) - 7*ARBCOMPLEX(33) - 7*ARBCOMPLEX(32))/(4*SQRT(3) + 7 )), ((2*SQRT(3)*( - ARBCOMPLEX(34) - 2*ARBCOMPLEX(33) - 2*ARBCOMPLEX(32)) - 3*ARBCOMPLEX(34) - 7*ARBCOMPLEX(33) - 7*ARBCOMPLEX(32))/(4*SQRT(3) + 7 )), (ARBCOMPLEX(32)), (ARBCOMPLEX(33)), (ARBCOMPLEX(34))) }, {X - 1, 3, [ 0 ] [ ] [ - ARBCOMPLEX(37)] [ ] [ - ARBCOMPLEX(36)] [ ] [ - ARBCOMPLEX(35)] [ ] [ ARBCOMPLEX(35) ] [ ] [ ARBCOMPLEX(36) ] [ ] [ ARBCOMPLEX(37) ] [ ] [ 0 ] }, {X + 7, 1, [ - ARBCOMPLEX(38) ] [ ------------------- ] [ 56*SQRT(3) + 97 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(38) - 24*ARBCOMPLEX(38) ] [--------------------------------------------------] [ 168*SQRT(3) + 291 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(38) - 24*ARBCOMPLEX(38) ] [--------------------------------------------------] [ 168*SQRT(3) + 291 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(38) - 24*ARBCOMPLEX(38) ] [--------------------------------------------------] [ 168*SQRT(3) + 291 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(38) - 24*ARBCOMPLEX(38) ] [--------------------------------------------------] [ 168*SQRT(3) + 291 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(38) - 24*ARBCOMPLEX(38) ] [--------------------------------------------------] [ 168*SQRT(3) + 291 ] [ ] [ - 14*SQRT(3)*ARBCOMPLEX(38) - 24*ARBCOMPLEX(38) ] [--------------------------------------------------] [ 168*SQRT(3) + 291 ] [ ] [ ARBCOMPLEX(38) ] }} for each eispace in eig do begin scalar eivaleq,eival,eivec; eival := solve(first eispace,x); for each soln in eival do <<eival := rhs soln; eivec := third eispace; eivec := sub(soln,eivec); write "eigenvalue = ", eival; write "check of eigen equation: ", m*eivec - eival*eivec>> end; eigenvalue = -1 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] eigenvalue = 1 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] eigenvalue = -7 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] % For this choice of values clear e; let e = 1; % the eigenvalue 1 becomes 4-fold degenerate. However, we get a complete % span of 8 eigenvectors. eig := mateigen(m,x); EIG := {{X - 1, 4, [ - ARBCOMPLEX(42)] [ ] [ - ARBCOMPLEX(41)] [ ] [ - ARBCOMPLEX(40)] [ ] [ - ARBCOMPLEX(39)] [ ] [ ARBCOMPLEX(39) ] [ ] [ ARBCOMPLEX(40) ] [ ] [ ARBCOMPLEX(41) ] [ ] [ ARBCOMPLEX(42) ] }, {X + 1, 3, [ ARBCOMPLEX(45) ] [ ] [ ARBCOMPLEX(44) ] [ ] [ ARBCOMPLEX(43) ] [ ] [ - (ARBCOMPLEX(45) + ARBCOMPLEX(44) + ARBCOMPLEX(43))] [ ] [ - (ARBCOMPLEX(45) + ARBCOMPLEX(44) + ARBCOMPLEX(43))] [ ] [ ARBCOMPLEX(43) ] [ ] [ ARBCOMPLEX(44) ] [ ] [ ARBCOMPLEX(45) ] }, {X - 7, 1, [ARBCOMPLEX(46)] [ ] [ARBCOMPLEX(46)] [ ] [ARBCOMPLEX(46)] [ ] [ARBCOMPLEX(46)] [ ] [ARBCOMPLEX(46)] [ ] [ARBCOMPLEX(46)] [ ] [ARBCOMPLEX(46)] [ ] [ARBCOMPLEX(46)] }} for each eispace in eig do begin scalar eivaleq,eival,eivec; eival := solve(first eispace,x); for each soln in eival do <<eival := rhs soln; eivec := third eispace; eivec := sub(soln,eivec); write "eigenvalue = ", eival; write "check of eigen equation: ", m*eivec - eival*eivec>> end; eigenvalue = 1 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] eigenvalue = -1 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] eigenvalue = 7 check of eigen equation: [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] [ ] [0] ma := mat((1,a),(0,b)); [1 A] MA := [ ] [0 B] % case 1: let a = 0; mateigen(ma,x); {{X - 1,1, [ARBCOMPLEX(47)] [ ] [ 0 ] }, { - B + X,1, [ 0 ] [ ] [ARBCOMPLEX(48)] }} % case 2: clear a; let a = 0, b = 1; mateigen(ma,x); {{X - 1,2, [ARBCOMPLEX(49)] [ ] [ARBCOMPLEX(50)] }} % case 3: clear a,b; mateigen(ma,x); {{ - B + X, 1, [ ARBCOMPLEX(51)*A ] [------------------] [ B - 1 ] [ ] [ ARBCOMPLEX(51) ] }, {X - 1,1, [ARBCOMPLEX(52)] [ ] [ 0 ] }} % case 4: let b = 1; mateigen(ma,x); {{X - 1,2, [ARBCOMPLEX(53)] [ ] [ 0 ] }} end; 5: 5: Time: 11713 ms plus GC time: 697 ms 6: 6: Quitting Sat May 30 16:12:03 PDT 1992