Artifact 5a5fc55a98aacf749bcf3f4a22b336c65a692c9517616f42c578ea77cc29e5b0:
- Executable file
r38/packages/matrix/matrix.rlg
— 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: 51621) [annotate] [blame] [check-ins using] [more...]
Tue Feb 10 12:26:21 2004 run on Linux % 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 <<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 ] }} % 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 for test: 230 ms, plus GC time: 10 ms