Artifact a8a2ef7118ffb817af9a25a10c0d81993946b2c0f2ede169e20489e77526430f:
- Executable file
r36/XMPL/MATRIX.TST
— 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: 8787) [annotate] [blame] [check-ins using] [more...]
% 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); 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; % Example of degenerate eigenvalues. u := mat((2,-1,1),(0,1,1),(-1,1,1))$ mateigen(u,eta); % 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); eigw := third first ws; w*eigw - al*eigw; % 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); % 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); % 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)); mateigen(j,x); % 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); % 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); % Tests of nullspace operator. a1 := mat((1,2,3,4),(5,6,7,8)); nullspace a1; b1 := {{1,2,3,4},{5,6,7,8}}; nullspace b1; % 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)))}}; nullspace c1; 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))))); nullspace d1; % 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)); eig := mateigen(m,x); % 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; % For the special choice: let e = -7 + sqrt 48; % we get only 7 eigenvectors. eig := mateigen(m,x); 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; % The same behaviour for this choice of e. clear e; let e = -7 - sqrt 48; % we get only 7 eigenvectors. eig := mateigen(m,x); 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; % 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); 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; ma := mat((1,a),(0,b)); % case 1: let a = 0; mateigen(ma,x); % case 2: clear a; let a = 0, b = 1; mateigen(ma,x); % case 3: clear a,b; mateigen(ma,x); % case 4: let b = 1; mateigen(ma,x); % Example from H.G. Graebe: m1:=mat((-sqrt(3)+1,2 ,3 ), (2 ,-sqrt(3)+3,1 ), (3 ,1 ,-sqrt(3)+2)); nullspace m1; for each n in ws collect m1*n; end;