Tue Apr 15 00:32:23 2008 run on win32
% 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: 259 ms, plus GC time: 8 ms