reduce-historical/r37/packages/matrix/matrix.rlg

Sun Jan  3 23:46:20 MET 1999
REDUCE 3.7, 15-Jan-99 ...

1: 1: 
2: 2: 2: 2: 2: 2: 2: 2: 2: 
3: 3: % 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;

4: 4: 4: 4: 4: 4: 4: 4: 4: 

Time for test: 1930 ms, plus GC time: 80 ms

5: 5: 
Quitting
Sun Jan  3 23:46:28 MET 1999