@@ -1,2200 +1,2200 @@ -REDUCE 3.6, 15-Jul-95, patched to 6 Mar 96 ... - - -% Miscellaneous matrix tests. - -% 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); - - -{{et - 2, - - 1, - - - [ 0 ] - [ ] - [ 0 ] - [ ] - [arbcomplex(1)] - [ ] - [arbcomplex(1)] - [ ] - [arbcomplex(1)] - - }, - - 4 3 2 - {et - 6*et + 13*et + 5*et - 5, - - 1, - - - [ 5*arbcomplex(2)*(et - 2) ] - [ ---------------------------- ] - [ 3 2 ] - [ 2*et - 10*et + 23*et + 5 ] - [ ] - [ 2 ] - [ arbcomplex(2)*et*( - et + 6*et - 8) ] - [--------------------------------------] - [ 3 2 ] - [ 2*et - 10*et + 23*et + 5 ] - [ ] - [ arbcomplex(2)*et*( - 3*et + 7) ] - [ -------------------------------- ] - [ 3 2 ] - [ 2*et - 10*et + 23*et + 5 ] - [ ] - [ 3 2 ] - [ arbcomplex(2)*(et - 4*et + 10) ] - [ ---------------------------------- ] - [ 3 2 ] - [ 2*et - 10*et + 23*et + 5 ] - [ ] - [ 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 ] -[ ] -[ 0 ] -[ ] -[arbcomplex(1)*( - et + 2)] -[ ] -[arbcomplex(1)*( - et + 2)] -[ ] -[arbcomplex(1)*( - et + 2)] - - - -% 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) ] -[ --------------------------- ] -[ 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 *(bx + by + bz + ex + ey + ez ) + bx *ex + 2*bx*by*ex*ey - - 2 2 2 2 - + 2*bx*bz*ex*ez + by *ey + 2*by*bz*ey*ez + bz *ez , - - 1, - - - 2 - mat(((om *arbcomplex(7)*ez + om*arbcomplex(7)*(bx*ey - by*ex) - - 3 2 2 2 - + arbcomplex(7)*bz*(bx*ex + by*ey + bz*ez))/(om + om*(bz + ex + ey ) - - )), - - 2 - (( - om *arbcomplex(7)*by + om*arbcomplex(7)*(bx*bz - ex*ez) - - 3 2 2 2 - - arbcomplex(7)*ey*(bx*ex + by*ey + bz*ez))/(om + om*(bz + ex + ey ) - - )), - - 2 - ((om *arbcomplex(7)*bx + om*arbcomplex(7)*(by*bz - ey*ez) - - 3 2 2 2 - + arbcomplex(7)*ex*(bx*ex + by*ey + bz*ez))/(om + om*(bz + ex + ey ) - - )), - - (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 + bx + bz + ey , - - 1, - - - [ - arbcomplex(10)*ey ] - [ ---------------------- ] - [ bx ] - [ ] - [ - arbcomplex(10)*bz ] - [ ---------------------- ] - [ bx ] - [ ] - [ 2 2 2 ] - [ arbcomplex(10)*(bx + bz + ey ) ] - [----------------------------------] - [ 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 *(p1 + p2 + p3 - s*z - z ) -c1 := {{----------------------------------, - 2 2 - p1 + p3 - - 0, - - 2 2 2 2 - p1*p3*(p1 + p2 + p3 - s*z - z ) - ------------------------------------, - 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*(p1 + p2 + p3 - s*z - z ) - {------------------------------------, - 2 2 - p1 + p3 - - 0, - - 2 2 2 2 2 - p3 *(p1 + p2 + p3 - s*z - z ) - ----------------------------------, - 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 ) - - - p1*p2*p3*z - ---------------, - 2 2 - p1 + p3 - - 0, - - 2 - p1 *p2*z - -----------}, - 2 2 - p1 + p3 - - {p1*(s + z), - - 0, - - p3*(s + z), - - p1*p2*(s + 2*z) - -----------------, - z - - 2 2 2 2 - - p1 *s - 2*p1 *z - p3 *s - 2*p3 *z - --------------------------------------, - z - - p2*p3*(s + 2*z) - -----------------, - z - - p3*z, - - 0, - - - p1*z}, - - - 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 - - p2*p3 *z - -------------, - 2 2 - p1 + p3 - - 0, - - p1*p2*p3*z - ------------}, - 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 - - - p1*p2*p3*z - ---------------, - 2 2 - p1 + p3 - - p3*z, - - 2 - - p2*p3 *z - -------------, - 2 2 - p1 + p3 - - 2 2 2 2 - - p3 *z*(p1 + p2 + p3 ) - -------------------------------, - 2 2 2 2 - p1 *s + p1 *z + p3 *s + p3 *z - - 0, - - 2 2 2 - p1*p3*z*(p1 + p2 + p3 ) - -------------------------------}, - 2 2 2 2 - p1 *s + p1 *z + p3 *s + p3 *z - - {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 - p1 *p2*z - -----------, - 2 2 - p1 + p3 - - - p1*z, - - p1*p2*p3*z - ------------, - 2 2 - p1 + p3 - - 2 2 2 - p1*p3*z*(p1 + p2 + p3 ) - -------------------------------, - 2 2 2 2 - p1 *s + p1 *z + p3 *s + p3 *z - - 0, - - 2 2 2 2 - - p1 *z*(p1 + p2 + p3 ) - -------------------------------}} - 2 2 2 2 - p1 *s + p1 *z + p3 *s + p3 *z - - -nullspace c1; - - -{{p3,0, - p1,0,0,0,0,0,0}, - - {0,1,0,0,0,0,0,0,0}, - - {0,0,0,p3,0, - p1,0,0,0}, - - 2 2 - {0,0,0,0,p2*p3,p1 + p3 ,0,0,0}, - - {0,0,0,0,0,0,p1,0,p3}, - - {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 *(p1 + p2 + p3 - s*z - z ) -d1 := mat((----------------------------------,0, - 2 2 - p1 + p3 - - 2 2 2 2 2 - p1*p3*(p1 + p2 + p3 - s*z - z ) - p1 *p2*(s + z) - ------------------------------------,-------------------,p1*(s + z), - 2 2 2 2 - p1 + p3 p1 + p3 - - 2 2 2 - - p1*p2*p3*(s + z) - p1*p3*(p1 + p2 + p3 ) - ---------------------,----------------------------,0, - 2 2 2 2 - p1 + p3 p1 + p3 - - 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*(p1 + p2 + p3 - s*z - z ) - (------------------------------------,0, - 2 2 - p1 + p3 - - 2 2 2 2 2 - p3 *(p1 + p2 + p3 - s*z - z ) - p1*p2*p3*(s + z) - ----------------------------------,---------------------,p3*(s + z), - 2 2 2 2 - p1 + p3 p1 + p3 - - 2 2 2 2 2 - - p2*p3 *(s + z) - p3 *(p1 + p2 + p3 ) - -------------------,--------------------------,0, - 2 2 2 2 - p1 + p3 p1 + p3 - - 2 2 2 - p1*p3*(p1 + p2 + p3 ) - -------------------------), - 2 2 - p1 + p3 - - 2 2 2 - p1 *p2*(s + z) - p1*p2*p3*(s + z) p1 *p2 *( - s - 2*z) - (----------------,0,---------------------,----------------------, - 2 2 2 2 2 2 - p1 + p3 p1 + p3 z*(p1 + p3 ) - - 2 - p1*p2*(s + 2*z) p1*p2 *p3*( - s - 2*z) - p1*p2*p3*z - -----------------,------------------------,---------------,0, - z 2 2 2 2 - z*(p1 + p3 ) p1 + p3 - - 2 - p1 *p2*z - -----------), - 2 2 - p1 + p3 - - p1*p2*(s + 2*z) - (p1*(s + z),0,p3*(s + z),-----------------, - z - - 2 2 2 2 - - p1 *s - 2*p1 *z - p3 *s - 2*p3 *z p2*p3*(s + 2*z) - --------------------------------------,-----------------,p3*z,0, - z z - - - p1*z), - - 2 2 - - p1*p2*p3*(s + z) - p2*p3 *(s + z) p1*p2 *p3*( - s - 2*z) - (---------------------,0,-------------------,------------------------, - 2 2 2 2 2 2 - p1 + p3 p1 + p3 z*(p1 + p3 ) - - 2 2 2 - p2*p3*(s + 2*z) p2 *p3 *( - s - 2*z) - p2*p3 *z p1*p2*p3*z - -----------------,----------------------,-------------,0,------------ - z 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 2 2 2 2 - - p1*p2*p3*z - p2*p3 *z - p3 *z*(p1 + p2 + p3 ) - ---------------,p3*z,-------------,-------------------------------,0, - 2 2 2 2 2 2 2 2 - p1 + p3 p1 + p3 p1 *s + p1 *z + p3 *s + p3 *z - - 2 2 2 - p1*p3*z*(p1 + p2 + p3 ) - -------------------------------), - 2 2 2 2 - p1 *s + p1 *z + p3 *s + p3 *z - - (0,0,0,0,0,0,0,0,0), - - 2 2 2 2 2 2 2 2 - p1 *(p1 + p2 + p3 ) p1*p3*(p1 + p2 + p3 ) p1 *p2*z - (-----------------------,0,-------------------------,-----------, - 2 2 2 2 2 2 - p1 + p3 p1 + p3 p1 + p3 - - 2 2 2 - p1*p2*p3*z p1*p3*z*(p1 + p2 + p3 ) - - p1*z,------------,-------------------------------,0, - 2 2 2 2 2 2 - p1 + p3 p1 *s + p1 *z + p3 *s + p3 *z - - 2 2 2 2 - - p1 *z*(p1 + p2 + p3 ) - -------------------------------)) - 2 2 2 2 - p1 *s + p1 *z + p3 *s + p3 *z - - - -nullspace d1; - - -{ - - [0] - [ ] - [1] - [ ] - [0] - [ ] - [0] - [ ] - [0] - [ ] - [0] - [ ] - [0] - [ ] - [0] - [ ] - [0] - - , - - - [ 0 ] - [ ] - [ 0 ] - [ ] - [ 0 ] - [ ] - [ p3 ] - [ ] - [ 0 ] - [ ] - [ - p1] - [ ] - [ 0 ] - [ ] - [ 0 ] - [ ] - [ 0 ] - - , - - - [ 0 ] - [ ] - [ 0 ] - [ ] - [ 0 ] - [ ] - [ 0 ] - [ ] - [ p2*p3 ] - [ ] - [ 2 2] - [p1 + p3 ] - [ ] - [ 0 ] - [ ] - [ 0 ] - [ ] - [ 0 ] - - , - - - [0 ] - [ ] - [0 ] - [ ] - [0 ] - [ ] - [0 ] - [ ] - [0 ] - [ ] - [0 ] - [ ] - [p1] - [ ] - [0 ] - [ ] - [p3] - - , - - - [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 6 - mat(((6*arbcomplex(24)*(e *x + 23*e *x - 7*e + 179*e *x - 119*e - - 5 5 4 4 3 3 - + 565*e *x - 581*e + 768*e *x - 890*e + 565*e *x - 581*e - - 2 2 3 8 7 - + 179*e *x - 119*e + 23*e*x - 7*e + x))/(e *(e *x + 30*e *x - - 7 6 6 5 5 - - 7*e + 333*e *x - 168*e + 1692*e *x - 1365*e - - 4 4 3 3 2 - + 4023*e *x - 4368*e + 4470*e *x - 5145*e + 2663*e *x - - 2 - - 2520*e + 576*e*x - 251*e + 36*x))), - - 9 8 8 7 7 - ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e - - 6 6 5 5 4 - + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x - - 4 3 3 2 2 - - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x - - 2 8 7 7 6 6 - - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e - - 5 5 4 4 3 - + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x - - 3 2 2 - - 5145*e + 2663*e *x - 2520*e + 576*e*x - 251*e + 36*x))) - - , - - 9 8 8 7 7 - ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e - - 6 6 5 5 4 - + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x - - 4 3 3 2 2 - - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x - - 2 8 7 7 6 6 - - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e - - 5 5 4 4 3 - + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x - - 3 2 2 - - 5145*e + 2663*e *x - 2520*e + 576*e*x - 251*e + 36*x))) - - , - - 9 8 8 7 7 - ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e - - 6 6 5 5 4 - + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x - - 4 3 3 2 2 - - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x - - 2 8 7 7 6 6 - - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e - - 5 5 4 4 3 - + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x - - 3 2 2 - - 5145*e + 2663*e *x - 2520*e + 576*e*x - 251*e + 36*x))) - - , - - 9 8 8 7 7 - ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e - - 6 6 5 5 4 - + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x - - 4 3 3 2 2 - - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x - - 2 8 7 7 6 6 - - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e - - 5 5 4 4 3 - + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x - - 3 2 2 - - 5145*e + 2663*e *x - 2520*e + 576*e*x - 251*e + 36*x))) - - , - - 9 8 8 7 7 - ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e - - 6 6 5 5 4 - + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x - - 4 3 3 2 2 - - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x - - 2 8 7 7 6 6 - - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e - - 5 5 4 4 3 - + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x - - 3 2 2 - - 5145*e + 2663*e *x - 2520*e + 576*e*x - 251*e + 36*x))) - - , - - 9 8 8 7 7 - ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e - - 6 6 5 5 4 - + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x - - 4 3 3 2 2 - - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x - - 2 8 7 7 6 6 - - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e - - 5 5 4 4 3 - + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x - - 3 2 2 - - 5145*e + 2663*e *x - 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 - <> - 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 - <> - 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 - <> - 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 - <> - 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 ] - - }} - - -% Example from H.G. Graebe: - -m1:=mat((-sqrt(3)+1,2 ,3 ), - (2 ,-sqrt(3)+3,1 ), - (3 ,1 ,-sqrt(3)+2)); - - - [ - sqrt(3) + 1 2 3 ] - [ ] -m1 := [ 2 - sqrt(3) + 3 1 ] - [ ] - [ 3 1 - sqrt(3) + 2] - - - -nullspace m1; - - -{ - - [ 3*sqrt(3) - 7 ] - [ ] - [ sqrt(3) + 5 ] - [ ] - [ - 4*sqrt(3) + 2] - - } - - -for each n in ws collect m1*n; - - -{ - - [0] - [ ] - [0] - [ ] - [0] - - } - - -end; -(TIME: matrix 5070 5619) +REDUCE 3.6, 15-Jul-95, patched to 6 Mar 96 ... + + +% Miscellaneous matrix tests. + +% 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); + + +{{et - 2, + + 1, + + + [ 0 ] + [ ] + [ 0 ] + [ ] + [arbcomplex(1)] + [ ] + [arbcomplex(1)] + [ ] + [arbcomplex(1)] + + }, + + 4 3 2 + {et - 6*et + 13*et + 5*et - 5, + + 1, + + + [ 5*arbcomplex(2)*(et - 2) ] + [ ---------------------------- ] + [ 3 2 ] + [ 2*et - 10*et + 23*et + 5 ] + [ ] + [ 2 ] + [ arbcomplex(2)*et*( - et + 6*et - 8) ] + [--------------------------------------] + [ 3 2 ] + [ 2*et - 10*et + 23*et + 5 ] + [ ] + [ arbcomplex(2)*et*( - 3*et + 7) ] + [ -------------------------------- ] + [ 3 2 ] + [ 2*et - 10*et + 23*et + 5 ] + [ ] + [ 3 2 ] + [ arbcomplex(2)*(et - 4*et + 10) ] + [ ---------------------------------- ] + [ 3 2 ] + [ 2*et - 10*et + 23*et + 5 ] + [ ] + [ 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 ] +[ ] +[ 0 ] +[ ] +[arbcomplex(1)*( - et + 2)] +[ ] +[arbcomplex(1)*( - et + 2)] +[ ] +[arbcomplex(1)*( - et + 2)] + + + +% 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) ] +[ --------------------------- ] +[ 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 *(bx + by + bz + ex + ey + ez ) + bx *ex + 2*bx*by*ex*ey + + 2 2 2 2 + + 2*bx*bz*ex*ez + by *ey + 2*by*bz*ey*ez + bz *ez , + + 1, + + + 2 + mat(((om *arbcomplex(7)*ez + om*arbcomplex(7)*(bx*ey - by*ex) + + 3 2 2 2 + + arbcomplex(7)*bz*(bx*ex + by*ey + bz*ez))/(om + om*(bz + ex + ey ) + + )), + + 2 + (( - om *arbcomplex(7)*by + om*arbcomplex(7)*(bx*bz - ex*ez) + + 3 2 2 2 + - arbcomplex(7)*ey*(bx*ex + by*ey + bz*ez))/(om + om*(bz + ex + ey ) + + )), + + 2 + ((om *arbcomplex(7)*bx + om*arbcomplex(7)*(by*bz - ey*ez) + + 3 2 2 2 + + arbcomplex(7)*ex*(bx*ex + by*ey + bz*ez))/(om + om*(bz + ex + ey ) + + )), + + (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 + bx + bz + ey , + + 1, + + + [ - arbcomplex(10)*ey ] + [ ---------------------- ] + [ bx ] + [ ] + [ - arbcomplex(10)*bz ] + [ ---------------------- ] + [ bx ] + [ ] + [ 2 2 2 ] + [ arbcomplex(10)*(bx + bz + ey ) ] + [----------------------------------] + [ 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 *(p1 + p2 + p3 - s*z - z ) +c1 := {{----------------------------------, + 2 2 + p1 + p3 + + 0, + + 2 2 2 2 + p1*p3*(p1 + p2 + p3 - s*z - z ) + ------------------------------------, + 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*(p1 + p2 + p3 - s*z - z ) + {------------------------------------, + 2 2 + p1 + p3 + + 0, + + 2 2 2 2 2 + p3 *(p1 + p2 + p3 - s*z - z ) + ----------------------------------, + 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 ) + + - p1*p2*p3*z + ---------------, + 2 2 + p1 + p3 + + 0, + + 2 + p1 *p2*z + -----------}, + 2 2 + p1 + p3 + + {p1*(s + z), + + 0, + + p3*(s + z), + + p1*p2*(s + 2*z) + -----------------, + z + + 2 2 2 2 + - p1 *s - 2*p1 *z - p3 *s - 2*p3 *z + --------------------------------------, + z + + p2*p3*(s + 2*z) + -----------------, + z + + p3*z, + + 0, + + - p1*z}, + + - 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 + - p2*p3 *z + -------------, + 2 2 + p1 + p3 + + 0, + + p1*p2*p3*z + ------------}, + 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 + + - p1*p2*p3*z + ---------------, + 2 2 + p1 + p3 + + p3*z, + + 2 + - p2*p3 *z + -------------, + 2 2 + p1 + p3 + + 2 2 2 2 + - p3 *z*(p1 + p2 + p3 ) + -------------------------------, + 2 2 2 2 + p1 *s + p1 *z + p3 *s + p3 *z + + 0, + + 2 2 2 + p1*p3*z*(p1 + p2 + p3 ) + -------------------------------}, + 2 2 2 2 + p1 *s + p1 *z + p3 *s + p3 *z + + {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 + p1 *p2*z + -----------, + 2 2 + p1 + p3 + + - p1*z, + + p1*p2*p3*z + ------------, + 2 2 + p1 + p3 + + 2 2 2 + p1*p3*z*(p1 + p2 + p3 ) + -------------------------------, + 2 2 2 2 + p1 *s + p1 *z + p3 *s + p3 *z + + 0, + + 2 2 2 2 + - p1 *z*(p1 + p2 + p3 ) + -------------------------------}} + 2 2 2 2 + p1 *s + p1 *z + p3 *s + p3 *z + + +nullspace c1; + + +{{p3,0, - p1,0,0,0,0,0,0}, + + {0,1,0,0,0,0,0,0,0}, + + {0,0,0,p3,0, - p1,0,0,0}, + + 2 2 + {0,0,0,0,p2*p3,p1 + p3 ,0,0,0}, + + {0,0,0,0,0,0,p1,0,p3}, + + {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 *(p1 + p2 + p3 - s*z - z ) +d1 := mat((----------------------------------,0, + 2 2 + p1 + p3 + + 2 2 2 2 2 + p1*p3*(p1 + p2 + p3 - s*z - z ) - p1 *p2*(s + z) + ------------------------------------,-------------------,p1*(s + z), + 2 2 2 2 + p1 + p3 p1 + p3 + + 2 2 2 + - p1*p2*p3*(s + z) - p1*p3*(p1 + p2 + p3 ) + ---------------------,----------------------------,0, + 2 2 2 2 + p1 + p3 p1 + p3 + + 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*(p1 + p2 + p3 - s*z - z ) + (------------------------------------,0, + 2 2 + p1 + p3 + + 2 2 2 2 2 + p3 *(p1 + p2 + p3 - s*z - z ) - p1*p2*p3*(s + z) + ----------------------------------,---------------------,p3*(s + z), + 2 2 2 2 + p1 + p3 p1 + p3 + + 2 2 2 2 2 + - p2*p3 *(s + z) - p3 *(p1 + p2 + p3 ) + -------------------,--------------------------,0, + 2 2 2 2 + p1 + p3 p1 + p3 + + 2 2 2 + p1*p3*(p1 + p2 + p3 ) + -------------------------), + 2 2 + p1 + p3 + + 2 2 2 + p1 *p2*(s + z) - p1*p2*p3*(s + z) p1 *p2 *( - s - 2*z) + (----------------,0,---------------------,----------------------, + 2 2 2 2 2 2 + p1 + p3 p1 + p3 z*(p1 + p3 ) + + 2 + p1*p2*(s + 2*z) p1*p2 *p3*( - s - 2*z) - p1*p2*p3*z + -----------------,------------------------,---------------,0, + z 2 2 2 2 + z*(p1 + p3 ) p1 + p3 + + 2 + p1 *p2*z + -----------), + 2 2 + p1 + p3 + + p1*p2*(s + 2*z) + (p1*(s + z),0,p3*(s + z),-----------------, + z + + 2 2 2 2 + - p1 *s - 2*p1 *z - p3 *s - 2*p3 *z p2*p3*(s + 2*z) + --------------------------------------,-----------------,p3*z,0, + z z + + - p1*z), + + 2 2 + - p1*p2*p3*(s + z) - p2*p3 *(s + z) p1*p2 *p3*( - s - 2*z) + (---------------------,0,-------------------,------------------------, + 2 2 2 2 2 2 + p1 + p3 p1 + p3 z*(p1 + p3 ) + + 2 2 2 + p2*p3*(s + 2*z) p2 *p3 *( - s - 2*z) - p2*p3 *z p1*p2*p3*z + -----------------,----------------------,-------------,0,------------ + z 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 2 2 2 2 + - p1*p2*p3*z - p2*p3 *z - p3 *z*(p1 + p2 + p3 ) + ---------------,p3*z,-------------,-------------------------------,0, + 2 2 2 2 2 2 2 2 + p1 + p3 p1 + p3 p1 *s + p1 *z + p3 *s + p3 *z + + 2 2 2 + p1*p3*z*(p1 + p2 + p3 ) + -------------------------------), + 2 2 2 2 + p1 *s + p1 *z + p3 *s + p3 *z + + (0,0,0,0,0,0,0,0,0), + + 2 2 2 2 2 2 2 2 + p1 *(p1 + p2 + p3 ) p1*p3*(p1 + p2 + p3 ) p1 *p2*z + (-----------------------,0,-------------------------,-----------, + 2 2 2 2 2 2 + p1 + p3 p1 + p3 p1 + p3 + + 2 2 2 + p1*p2*p3*z p1*p3*z*(p1 + p2 + p3 ) + - p1*z,------------,-------------------------------,0, + 2 2 2 2 2 2 + p1 + p3 p1 *s + p1 *z + p3 *s + p3 *z + + 2 2 2 2 + - p1 *z*(p1 + p2 + p3 ) + -------------------------------)) + 2 2 2 2 + p1 *s + p1 *z + p3 *s + p3 *z + + + +nullspace d1; + + +{ + + [0] + [ ] + [1] + [ ] + [0] + [ ] + [0] + [ ] + [0] + [ ] + [0] + [ ] + [0] + [ ] + [0] + [ ] + [0] + + , + + + [ 0 ] + [ ] + [ 0 ] + [ ] + [ 0 ] + [ ] + [ p3 ] + [ ] + [ 0 ] + [ ] + [ - p1] + [ ] + [ 0 ] + [ ] + [ 0 ] + [ ] + [ 0 ] + + , + + + [ 0 ] + [ ] + [ 0 ] + [ ] + [ 0 ] + [ ] + [ 0 ] + [ ] + [ p2*p3 ] + [ ] + [ 2 2] + [p1 + p3 ] + [ ] + [ 0 ] + [ ] + [ 0 ] + [ ] + [ 0 ] + + , + + + [0 ] + [ ] + [0 ] + [ ] + [0 ] + [ ] + [0 ] + [ ] + [0 ] + [ ] + [0 ] + [ ] + [p1] + [ ] + [0 ] + [ ] + [p3] + + , + + + [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 6 + mat(((6*arbcomplex(24)*(e *x + 23*e *x - 7*e + 179*e *x - 119*e + + 5 5 4 4 3 3 + + 565*e *x - 581*e + 768*e *x - 890*e + 565*e *x - 581*e + + 2 2 3 8 7 + + 179*e *x - 119*e + 23*e*x - 7*e + x))/(e *(e *x + 30*e *x + + 7 6 6 5 5 + - 7*e + 333*e *x - 168*e + 1692*e *x - 1365*e + + 4 4 3 3 2 + + 4023*e *x - 4368*e + 4470*e *x - 5145*e + 2663*e *x + + 2 + - 2520*e + 576*e*x - 251*e + 36*x))), + + 9 8 8 7 7 + ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e + + 6 6 5 5 4 + + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x + + 4 3 3 2 2 + - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x + + 2 8 7 7 6 6 + - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e + + 5 5 4 4 3 + + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x + + 3 2 2 + - 5145*e + 2663*e *x - 2520*e + 576*e*x - 251*e + 36*x))) + + , + + 9 8 8 7 7 + ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e + + 6 6 5 5 4 + + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x + + 4 3 3 2 2 + - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x + + 2 8 7 7 6 6 + - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e + + 5 5 4 4 3 + + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x + + 3 2 2 + - 5145*e + 2663*e *x - 2520*e + 576*e*x - 251*e + 36*x))) + + , + + 9 8 8 7 7 + ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e + + 6 6 5 5 4 + + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x + + 4 3 3 2 2 + - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x + + 2 8 7 7 6 6 + - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e + + 5 5 4 4 3 + + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x + + 3 2 2 + - 5145*e + 2663*e *x - 2520*e + 576*e*x - 251*e + 36*x))) + + , + + 9 8 8 7 7 + ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e + + 6 6 5 5 4 + + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x + + 4 3 3 2 2 + - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x + + 2 8 7 7 6 6 + - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e + + 5 5 4 4 3 + + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x + + 3 2 2 + - 5145*e + 2663*e *x - 2520*e + 576*e*x - 251*e + 36*x))) + + , + + 9 8 8 7 7 + ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e + + 6 6 5 5 4 + + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x + + 4 3 3 2 2 + - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x + + 2 8 7 7 6 6 + - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e + + 5 5 4 4 3 + + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x + + 3 2 2 + - 5145*e + 2663*e *x - 2520*e + 576*e*x - 251*e + 36*x))) + + , + + 9 8 8 7 7 + ((arbcomplex(24)*(e *x + 29*e *x - 7*e + 310*e *x - 161*e + + 6 6 5 5 4 + + 1520*e *x - 1246*e + 3577*e *x - 3836*e + 4283*e *x + + 4 3 3 2 2 + - 4795*e + 2988*e *x - 3065*e + 978*e *x - 672*e + 132*e*x + + 2 8 7 7 6 6 + - 42*e + 6*x))/(e *(e *x + 30*e *x - 7*e + 333*e *x - 168*e + + 5 5 4 4 3 + + 1692*e *x - 1365*e + 4023*e *x - 4368*e + 4470*e *x + + 3 2 2 + - 5145*e + 2663*e *x - 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 + <> + 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 + <> + 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 + <> + 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 + <> + 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 ] + + }} + + +% Example from H.G. Graebe: + +m1:=mat((-sqrt(3)+1,2 ,3 ), + (2 ,-sqrt(3)+3,1 ), + (3 ,1 ,-sqrt(3)+2)); + + + [ - sqrt(3) + 1 2 3 ] + [ ] +m1 := [ 2 - sqrt(3) + 3 1 ] + [ ] + [ 3 1 - sqrt(3) + 2] + + + +nullspace m1; + + +{ + + [ 3*sqrt(3) - 7 ] + [ ] + [ sqrt(3) + 5 ] + [ ] + [ - 4*sqrt(3) + 2] + + } + + +for each n in ws collect m1*n; + + +{ + + [0] + [ ] + [0] + [ ] + [0] + + } + + +end; +(TIME: matrix 5070 5619)