Overview
1°) Characteristic Polynomial & Minimal Polynomial
a) Program#1
b) Program#2
c) Minimal Polynomial
2°) The Power Method: Dominant Eigenvalue
a) Program#1
b) Program#2
3°) Deflation
4°) The Jacobi's Method
a) Symmetric ( and some non-symmetric
) Matrices
b) Symmetric Matrices only
5°) Complex Matrices
-The following programs - except the last one - only deal with real n x n matrices A
1°) "CP" calculates the characteristic polynomial of A (
n < 11 )
"CP2" ------------------------------------------
( n < 17 )
"PMIN" ----------
minimal ----------------
( n < 13 )
2°) "EGV" & "EGV2" give the dominant real eigenvalue
of A and a corresponding eigenvector by the power method
( n < 17 )
3°) "DFL" uses a deflation method to produce the real
eigenvalues ki and the eigenvectors, provided that | k1
| > | k2 | > ....... > | kn | > 0
( n < 17 )
4°) "JCB" uses the Jacobi's method to produce all the n
real eigenvalues and eigenvectors of a symmetric matrix of order
n ( n < 13 )
-However, it also works
if the matrix is not symmetric provided all the eigenvalues
are real and distinct. ( n < 13 )
"JCB2" runs faster and uses
less data registers but it works for symmetric matrices only ( n
< 15 )
5°) "JCBZ" uses a generalization of the Jacobi's method
to compute all the eigenvalues and eigenvectors of a Hermitian matrix.
( n < 9 )
-It also works for a general
complex matrix provided all the eigenvalues are distinct.
( n < 9 )
-All these programs - except "PMIN" - use the REGMOVE function
of the X-functions module.
1°) Characteristic Polynomial & Minimal
Polynomial
a) Program#1
-If there exist a number k and a non-zero vector V such that A.V
= k.V k is called an eigenvalue and V an eigenvector
of the matrix A.
-The characteristic polynomial p(x) = xn + c1.xn-1
+ c2.xn-2 + ......... + cn-2.x2
+ cn-1.x + cn is defined by p(x) = det ( x.I - A
)
where I is the n x n unit matrix.
-Thus, its roots are actually the n eigenvalues of A.
-"CP" uses the following formulae: we define a sequence of matrices Mk by:
M0 = I and
Mk+1 = Mk.A - trace(Mk.A) I / (k+1)
k = 0 , 1 , .... ( the trace of a matrix
is the sum of its diagonal elements )
and we have:
ck = - trace(Mk-1.A)/k
-The elements of A are to be stored into contiguous registers from R01
to Rn2 , n being stored in R00
-When the program stops, R01 = c1 , R02 = c2
, ............... , Rnn = cn
( c1 = -trace(A) and cn = (-1)n
det(A) )
Data Registers:
• R00 = n
• R01 = a11 ......................
• Rn2-n+1 = a1n
...................................................................
( the n2 elements of A )
INPUTS
• Rnn = an1 ......................
• Rn2 = ann
---------
R00 = n
R01 = c1 Rn+1 =
a11 Rn2+1
= a1n
.............
..........................................
OUTPUTS
Rnn = cn
R2n = an1
Rn2+n = ann
Flag: F10
Subroutine: "M*" ( cf "Matrix Operations
for the HP-41" )
| 01 LBL "CP"
02 RCL 00 03 ENTER^ 04 X^2 05 E3 06 / 07 ISG X 08 + 09 E3 10 / 11 1 12 + 13 REGMOVE 14 RCL 00 15 .1 16 % 17 STO IND 00 18 ISG X |
19 *
20 + 21 REGMOVE 22 RCL 00 23 .1 24 % 25 + 26 RCL 00 27 X^2 28 .1 29 % 30 ST+ X 31 + 32 + 33 ISG X 34 RCL 00 35 E5 36 / |
37 +
38 CF 10 39 LBL 01 40 ENTER^ 41 ISG IND 00 42 E-5 43 + 44 STO M 45 X<>Y 46 0 47 LBL 02 48 RCL IND Z 49 - 50 ISG Z 51 GTO 02 52 RCL IND 00 53 INT 54 ST/ Y |
55 RDN
56 STO IND L 57 ISG L 58 FS? 30 59 GTO 04 60 LBL 03 61 ST+ IND M 62 ISG M 63 GTO 03 64 X<>Y 65 STO Y 66 RCL 00 67 X^2 68 .1 69 % 70 + 71 FS? 10 72 ST+ X |
73 -
74 ENTER^ 75 INT 76 RCL 00 77 X^2 78 FC? 10 79 ST+ X 80 + 81 XEQ "M*" 82 FC?C 10 83 SF 10 84 GTO 01 85 LBL 04 86 CF 10 87 CLA 88 END |
( 138 bytes / SIZE 3n2+n+1 )
| STACK | INPUT | OUTPUT |
| X | / | cn |
Example: Find the characteristic polynomial
of the matrix
1 2 4
A = 4 3 5
7 4 7
-First we store these 9 elements in registers R01 thru R09
R01 R04 R07
1 2 4
R02 R05 R08
= 4 3 5 respectively
and n = 3 STO 00
R03 R06 R09
7 4 7
-Then XEQ "CP" >>>> c3 = 5 = R03
( in 31 seconds )
and R01 = -11 , R02 = -25 ,
R03 = 5 Thus p(x) = x3 - 11 x2
- 25 x + 5
Notes:
-A is saved in registers Rn+1 thru Rn2+n and n is saved in
R00
-The elements of A may be stored either in column order or in row order
since A and tA ( the transpose of A ) have the same characteristic
polynomial.
-Execution time = 1mn25s for n = 4
--------------- = 3mn18s for n = 5
-Any root-finding program yields the 3 eigenvalues
k1 = 12.90692994 , k2
= -2.092097589 , k3 = 0.185167648
-Then, subtracting ( for instance ) k1 to the diagonal elements
of A leads to a singular system with an infinite set of solutions ( the
eigenvectors )
-One solution is ( 0.451320606 ; 0.686921424 ;
1 ) which is an eigenvector corresponding to the first
eigenvalue.
-This method could be used to obtain all the eigenvalues and eigenvectors
of a n x n matrix A.
b) Program#2 ( HP-41CX
)
-"CP" uses 3 blocks of n2 registers to store the matrices
A , Mk and the products Mk.A so that n cannot exceed
10.
-"CP2" creates a temporary data file to store the matrix A and uses
another method to perform the multiplication of 2 square-matrices:
-Thus, it can produce the characteristic polynomial of a 16x16 matrix!
Data Registers:
• R00 = n
• R01 = a11 ......................
• Rn2-n+1 = a1n
..................................................................
( the n2 elements of A )
INPUTS
• Rnn = an1 ......................
• Rn2 = ann
---------
R00 = n
OUTPUTS R01 = c1
...................... Rnn = cn
Flags: /
Subroutine: "TRN2" Transpose of a square
matrix ( cf "Matrix Operations for the HP-41" )
-Line 110 is a three-byte GTO 01
-Do not key in lines 112-113 if you want to save the matrix A in this
data file.
| 01 LBL "CP2"
02 RCL 00 03 ENTER^ 04 X^2 05 "T" 06 CRFLD 07 E3 08 / 09 ISG X 10 SAVERX 11 + 12 E3 13 / 14 1 15 + 16 REGMOVE 17 RCL 00 18 .1 19 % 20 STO IND 00 21 1.001 22 + 23 * |
24 ISG X
25 RCL 00 26 E5 27 / 28 + 29 STO M 30 XEQ "TRN2" 31 LBL 01 32 ISG IND 00 33 RCL M 34 E-5 35 + 36 ENTER^ 37 ENTER^ 38 CLX 39 LBL 02 40 RCL IND Y 41 - 42 ISG Y 43 GTO 02 44 RCL IND 00 45 INT 46 ST/ Y |
47 X<>Y
48 STO IND L 49 ISG L 50 FS? 30 51 GTO 07 52 LBL 03 53 ST+ IND T 54 ISG T 55 GTO 03 56 R^ 57 INT 58 STO Y 59 RCL 00 60 ST- Z 61 SIGN 62 - 63 E3 64 / 65 + 66 STO N 67 RCL 00 68 .1 69 % |
70 ST+ X
71 + 72 ISG X 73 STO O 74 LBL 04 75 CLX 76 SEEKPT 77 LBL 05 78 RCL O 79 0 80 LBL 06 81 RCL IND Y 82 GETX 83 * 84 + 85 ISG Y 86 GTO 06 87 STO IND N 88 ISG N 89 GTO 05 90 RCL O 91 INT 92 RCL 00 |
93 ST- N
94 .1 95 % 96 + 97 ST+ O 98 X<> L 99 + 100 E3 101 / 102 RCL N 103 INT 104 + 105 REGMOVE 106 RCL N 107 RCL O 108 X#Y? 109 GTO 04 110 GTO 01 111 LBL 07 112 "T" 113 PURFL 114 CLA 115 END |
( 192 bytes / SIZE n2+2n+1
)
| STACK | INPUT | OUTPUT |
| X | / | cn |
Example1: Find the characteristic polynomial
of the matrix
1 2 4
A = 4 3 5
7 4 7
-First we store these 9 elements in registers R01 thru R09
R01 R04 R07
1 2 4
R02 R05 R08
= 4 3 5 respectively
and n = 3 STO 00
R03 R06 R09
7 4 7
-Then XEQ "CP2" >>>> c3 = 5 = R03
( in 37 seconds )
and R01 = -11 , R02 = -25 ,
R03 = 5 whence p(x) = x3
- 11 x2 - 25 x + 5
Example2: The 10x10 matrix A = [ aij ] is defined by aij = i j mod 13
-Store these 100 elements into registers R01 thru R100 ( you can
use for instance the "MATRIX" program listed in "Creating a Matrix for
the HP-41"
and the subroutine LBL "U" X<>Y Y^X
13 MOD RTN then "U" ASTO 00 1.10010
XEQ "MATRIX" )
-Then 10 STO 00 XEQ "CP2" and you should find ( in about 54mn45s ):
p(x) = x10 - 43 x9 - 968 x8 - 2462 x7 + 40796 x6 - 488852 x5 - 10916340 x4 + 15630136 x3 + 441980832 x2 - 1282786560 x + 155105280
-There is no roundoff error: all the coefficients are exact.
-Obviously, this is not a very fast routine for large n-values!
-For n = 16 , the execution time is probably greater than 5 hours,
so a good emulator like V41 in "turbo mode" is quite useful...
c) The Minimal Polynomial
-The characteristic polynomial P verifies P(A) = 0 ( Cayley-Hamilton
theorem )
-The minimal polynomial p is the unitary least degree
polynomial that satisfies p(A) = 0.
-It is a divisor of the characteristic polynomial.
-In the following program, we start with a random vector V and we compute
AV , A2V , ........... , AnV
-Then, "PMIN" seeks the smallest k such that V , AV , A2V
, ........... , AkV are linearly dependent.
-A relation a0 V + a1 AV + ........... +
ak-1 Ak-1V + AkV = 0 is found, whence
p(x)
Data Registers: • R00 = n ( Registers R00 thru Rn2 are to be initialized before executing "PMIN" )
• R01 = a11 ......................
• Rn2-n+1 = a1n
...................................................................
( the n2 elements of A )
• Rnn = an1 ...................... • Rn2 = ann Rn2+1 to R2n2+n: temp
>>>> When the program stops, R01 thru Ree = the coefficients of p(x)
Flags: /
Subroutine: "M*" "MCO"
( cf "Matrix Operations for the HP-41" )
"LS3" ( cf "Linear and Non-Linear Systems for the
HP-41" )
-Lines 11-12 may be replaced by the M-code function RAND ( cf
"A few M-Code Routines fot the HP-41" )
| 01 LBL "PMIN"
02 RCL 00 03 X^2 04 ST+ X 05 .1 06 % 07 + 08 RCL 00 09 ST+ Y 10 LBL 00 11 R-D 12 FRC 13 STO IND Y 14 DSE Y 15 GTO 00 16 CLX 17 RCL 00 18 .00101 19 * 20 + 21 1 22 + |
23 ENTER^
24 LBL 01 25 CLX 26 .01 27 RCL 00 28 ST+ Y 29 * 30 E3 31 / 32 1 33 + 34 X<>Y 35 ENTER^ 36 INT 37 RCL 00 38 - 39 XEQ "M*" 40 STO Y 41 INT 42 DSE X 43 RCL 00 44 X^2 |
45 -
46 X#0? 47 GTO 01 48 LASTX 49 ST+ X 50 RCL 00 51 + 52 E3 53 / 54 RCL 00 55 X^2 56 + 57 1 58 ST+ Y 59 XEQ "MCO" 60 RCL 00 61 RCL 00 62 1 63 + 64 XEQ "LS3' 65 R^ 66 INT |
67 STO 00
68 1 69 + 70 E5 71 / 72 RCL 00 73 X^2 74 + 75 RCL 00 76 + 77 SIGN 78 ENTER^ 79 LBL 02 80 RCL IND L 81 X#Y? 82 GTO 03 83 RDN 84 ST+ Y 85 DSE L 86 GTO 02 87 ENTER^ 88 LBL 03 |
89 LASTX
90 INT 91 R^ 92 ST+ Y 93 ENTER^ 94 DSE Z 95 LBL 04 96 RCL IND Z 97 CHS 98 STO IND Z 99 RDN 100 DSE Z 101 "" 102 DSE Y 103 GTO 04 104 E3 105 / 106 1 107 STO 01 108 + 109 END |
( 164 bytes / SIZE 2n2+n+1
)
| STACK | INPUT | OUTPUT |
| X | / | 1.eee |
Example: Find the minimal polynomial of the
matrix
[ [
3 -1 1 ]
A = [ -1 3 1
]
[ 1 1 3 ] ]
-Store
3 -1 1
R01 R04 R07
-1 3 1
into R02 R05 R08
respectively and
n = 3 STO 00
1 1 3
R03 R06 R09
XEQ "PMIN" >>>> 1.003 --- Execution time = 36s --- and we get R01 = 1 R02 = -5 R03 = 4 ( with some roundoff-errors )
-So, p(x) = x2 - 5 x + 4
( the characteristic polynomial is P(x) = x3 - 9 x2
+ 24 x - 16 = ( x - 4 ) p(x) )
Notes:
-Since we use a pseudo-random vector V, there is a small risk of getting
an eigenvector, which would lead to a wrong result.
-In practice however, that seems highly improbable.
-If the coefficients of A are integers, the coefficients of p are integers
too.
-Roundoff-errors may be important, especially if the system is ill-conditioned,
and E-7 is sometimes too small to identify tiny pivots ( cf
"LS3" listing, line 61 ).
-If you prefer "LS2" instead of "LS3", use the listing below:
-Lines 13-14 may be replaced by the M-code function RAND ( cf
"A few M-Code Routines fot the HP-41" )
| 01 LBL "PMIN"
02 RCL 00 03 X^2 04 STO Y 05 LASTX 06 + 07 E3 08 / 09 + 10 ISG X 11 RCL 00 12 LBL 00 13 R-D 14 FRC 15 STO IND Y 16 ISG Y 17 GTO 00 18 CLX 19 RCL 00 20 ST- Y 21 E5 22 / 23 + |
24 ENTER^
25 LBL 01 26 CLX 27 .01 28 RCL 00 29 ST+ Y 30 * 31 E3 32 / 33 1 34 + 35 X<>Y 36 ENTER^ 37 INT 38 RCL 00 39 + 40 XEQ "M*" 41 STO Y 42 INT 43 DSE X 44 RCL 00 45 X^2 46 ST+ X |
47 -
48 X#0? 49 GTO 01 50 LASTX 51 RCL 00 52 + 53 E3 54 / 55 RCL 00 56 X^2 57 + 58 1 59 ST+ Y 60 XEQ "MCO" 61 RCL 00 62 RCL 00 63 1 64 + 65 XEQ "LS2' 66 X<> Z 67 INT 68 ENTER^ 69 SQRT |
70 STO 00
71 1 72 + 73 E2 74 / 75 + 76 E3 77 / 78 1 79 + 80 SIGN 81 ENTER^ 82 LBL 02 83 RCL IND L 84 X#Y? 85 GTO 03 86 RDN 87 ST+ Y 88 ISG L 89 GTO 02 90 ENTER^ 91 LBL 03 92 LASTX |
93 INT
94 R^ 95 E3 96 / 97 2 98 + 99 DSE Y 100 LBL 04 101 RCL IND Y 102 CHS 103 STO IND Y 104 RDN 105 DSE Y 106 "" 107 ISG X 108 GTO 04 109 FRC 110 1 111 STO 01 112 + 113 END |
( 165 bytes / SIZE 2n2+n+1
)
| STACK | INPUT | OUTPUT |
| X | / | 1.eee |
-The same example gives the same results - with small differences in
the last decimals.
2°) The Power method: Dominant Eigenvalue
a) Program#1
-Here, we assume that the matrix A has n eigenvalues k1
, k2 , ....... , kn corresponding to n independent
eigenvectors V1 , V2 , ....... , Vn
and that k1 is the dominant eigenvalue | k1
| > | ki | i = 2 , 3 , ....... , n
-Any vector V can be written V = a1
V1 + a2 V2 + ....... + an Vn
whence AV = a1 k1 V1 + a2
k2 V2 + ....... + an kn Vn
and ApV = k1p ( a1
V1 + a2 ( k2/k1 )p
V2 + ....... + an ( kn/k1 )p
Vn ) which tends to k1p
a1 V1 as p tends to infinity.
-Furthermore, if we define Wp+1 = ( A.Wp
) / || A.Wp || with W0 = V ( almost
) arbitrary , then Wp tends to a unit eigenvector U
and tWp AWp tends to
the corresponding eigenvalue k1
Data Registers:
• R00 = n
• R01 = v1
• Rn+1 = a11
• Rn2+1 = a1n
.................
( the n2 elements of A )
where V ( v1 ; ....... ; vn )
is an initial guess of the eigenvector.
INPUTS
• Rnn = vn
• R2n = an1
• Rn2+n = ann
---------
R00 = n
R01 = u1 Rn+1 =
a11 Rn2+1
= a1n
.............
.............................................
and
k in X-register. U ( u1; ........ ; un)
= a unit eigenvector
OUTPUTS
Rnn = un
R2n = an1
Rn2+n = ann
Flags: /
Subroutine: "M*" ( cf "Matrix Operations
for the HP-41" )
| 01 LBL "EGV"
02 LBL 01 03 RCL 00 04 X^2 05 E3 06 / 07 RCL 00 08 ST+ Y 09 1 10 % 11 + 12 E3 13 / 14 1 15 + 16 ST+ Y |
17 LASTX
18 RCL 00 19 ST+ Y 20 X^2 21 + 22 XEQ "M*" 23 LASTX 24 RCL 00 25 E3 26 / 27 ISG X 28 LASTX 29 / 30 + 31 REGSWAP 32 INT |
33 RCL 00
34 STO N 35 + 36 DSE X 37 STO O 38 0 39 LBL 02 40 RCL IND Y 41 RCL IND N 42 X^2 43 ST+ M 44 X<> L 45 * 46 + 47 DSE Y 48 DSE N |
49 GTO 02
50 VIEW X 51 STO N 52 RCL 00 53 RCL 00 54 RCL M 55 SQRT 56 R^ 57 SIGN 58 * 59 LBL 03 60 ST/ IND Y 61 DSE Y 62 GTO 03 63 RDN 64 LBL 04 |
65 RCL IND Y
66 RCL IND O 67 - 68 ABS 69 + 70 DSE O 71 DSE Y 72 GTO 04 73 E-8 74 X<Y? 75 GTO 01 76 RCL N 77 CLA 78 END |
( 122 bytes / SIZE n2+2n+1 )
| STACK | INPUTS | OUTPUTS |
| X | / | k |
Example: Find the dominant eigenvalue k and
the corresponding unit eigenvector U of the matrix
1
2 4
A = 4 3 5
7 4 7
-Here, we have n = 3 therefore 3
STO 00
and if we choose ( 1 ; 1 ; 1 ) as initial guess, we store the
12 numbers
1 1 2 4
R01 R04 R07 R10
1 4 3 5 into
R02 R05 R08 R11 respectively.
1 7 4 7
R03 R06 R09 R12
-Then XEQ "EGV" the HP-41 displays the successive k-approximations
and 104 seconds later, it yields k = 12.90692995
in X-register and U is in registers R01 R02 R03 :
U ( 0.348663346 ; 0.530674468 ; 0.772540278 )
Notes:
-The greater | k1/ k2 | , the faster the convergence.
-If you calculate A-1 ( provided A is non singular ) and
apply "EGV" again , you'll find the reciprocal of the smallest eigenvalue
and the corresponding eigenvector.
in the above example, 1/k3 = 5.400511417
and the eigenvector ( -0.094824736 ; 0.897989404 ; -0.429678136
) which is also an eigenvector of A.
-Subtracting a number N to the diagonal elements of A can make the
other extreme eigenvalue dominant:
for instance choosing N = 13 and applying "EGV" leads
to k2 - 13 = -15.092... ( add 13 to the result
) and the last eigenvector is correctly found.
-Unfortunately in our example, the convergence is very slow
since the new ratio | k'1/ k'2 | = 1.17...
is very close to 1. More than 100 iterations are needed
to achieve a good accuracy! All the digits of k are correct
, but the components of U are only obtained to 7 decimals.
-Thus, in this case, the convergence of k is faster ( I mean less slow
) than the convergence of U:
That's why the test ( line 74 ) compares 2 successive
approximations of U instead of 2 successive approximations of k.
Otherwise, the program could have been both shorter and faster:
see below!
b) Program#2
-The following program uses the same "power method" but the iteration
stops when 2 successive k-values are (approximately ) equal.
-The eigenvector V( x1 , ....... , xn
) is choosen so that its first coordinate = the eigenvalue k.
-So, it will not work if x1 = 0
| 01 LBL "EGV2"
02 LBL 01 03 RCL 00 04 X^2 05 E3 06 / 07 RCL 00 08 ST+ Y 09 1 10 % |
11 +
12 E3 13 / 14 1 15 + 16 ST+ Y 17 LASTX 18 RCL 00 19 ST+ Y 20 X^2 |
21 +
22 XEQ "M*" 23 LASTX 24 RCL 00 25 E3 26 / 27 ISG X 28 LASTX 29 / 30 + |
31 REGSWAP
32 RCL 00 33 RCL IND Y 34 LBL 02 35 ST/ IND Y 36 DSE Y 37 GTO 02 38 RCL 01 39 VIEW X 40 ST- Y |
41 /
42 ABS 43 2 E-9 44 X<Y? 45 GTO 01 46 RCL 01 47 END |
( 77 bytes / SIZE n2+2n+1
)
| STACK | INPUTS | OUTPUTS |
| X | / | k |
-With the same example as "EGV" it yields k = 12.90692994
in X-register in 90 seconds
and V is in registers R01-R02-R03 : V ( 12.90692994 ; 19.64467513 ; 28.59814010 )
-With this less reliable method, always check that the eigenvector V
is also in registers Rn2+n+1 ...... Rn2+2n
( in this example, R13 R14 R15 )
-Of course, the last decimals may differ slightly.
3°) Deflation
-If k1 is an eigenvalue , U1 the corresponding
unit eigenvector of A and U1' the corresponding unit eigenvector
of tA
then A and the new matrix A' = A - k1
( U1 tU1' ) / ( U1.U1'
) have the same eigenvalues and eigenvectors , except that
k1 has been replaced by 0.
-Thus, the 2nd dominant eigenvalue of A is now the 1st
dominant eigenvalue of A'
-"DFL" uses this method to calculate all the eigenvalues and eigenvectors , provided they are real and verify | k1 | > | k2 | > ....... > | kn | > 0
1°) "DFL" calculates k1 and U1 ( with flag F02 clear ) and the program stops ( line 87 )
>>>> k1
is in X-register
>>>> U1
is in registers R01 thru Rnn
2°) If you press R/S , the HP-41 sets flag F02 and
computes k1 once again & U1'
the matrix A is replaced
by A' in registers Rn+1 thru Rn2+n and when the program
stops again:
>>>> k2
is in X-register
>>>> U2
is in registers R01 thru Rnn
... etc ...
Data Registers:
• R00 = n
• R01 = v1
• Rn+1 = a11
• Rn2+1 = a1n
.................
( the n2 elements of A )
where V ( v1 ; ....... ; vn )
is an initial guess of the eigenvector.
INPUTS
• Rnn = vn
• R2n = an1
• Rn2+n = ann
---------
R00 = n
R01 = u1 Rn+1 =
a'11 Rn2+1
= a'1n
.............
A will be replaced with
k in X-register. U ( u1 , ........
, un ) = a unit eigenvector
OUTPUTS
A' = A - k ( U.tU' ) / ( U.U' ) when F02 is set.
Rnn = un
R2n = a'n1
Rn2+n = a'nn
Flag: F02
Subroutine: "M*" ( cf "Matrix Operations
for the HP-41" )
-Lines 82-102 are three-byte GTO 01
-Line 142 is a three-byte GTO 12
| 01 LBL "DFL"
02 LBL 12 03 CF 02 04 LBL 01 05 RCL 00 06 X^2 07 LASTX 08 ST+ Y 09 E3 10 ST/ Z 11 / 12 1 13 ST+ Z 14 % 15 ISG Y 16 ST+ Z 17 FS? 02 18 CLX 19 + 20 ISG Y 21 FS? 02 22 X<>Y 23 RCL 00 24 ENTER^ 25 X^2 26 + 27 1 28 + 29 XEQ "M*" |
30 LASTX
31 RCL 00 32 E3 33 / 34 ISG X 35 LASTX 36 / 37 + 38 REGSWAP 39 INT 40 RCL 00 41 STO N 42 + 43 DSE X 44 STO O 45 0 46 LBL 02 47 RCL IND Y 48 RCL IND N 49 X^2 50 ST+ M 51 X<> L 52 * 53 + 54 DSE Y 55 DSE N 56 GTO 02 57 VIEW X 58 STO N |
59 RCL 00
60 RCL 00 61 RCL M 62 SQRT 63 R^ 64 SIGN 65 * 66 LBL 03 67 ST/ IND Y 68 DSE Y 69 GTO 03 70 RDN 71 LBL 04 72 RCL IND Y 73 RCL IND O 74 - 75 ABS 76 + 77 DSE O 78 DSE Y 79 GTO 04 80 E-8 81 X<Y? 82 GTO 01 83 RCL N 84 FS? 02 85 GTO 05 86 CLA 87 RTN |
88 SF 02
89 RCL 00 90 1 91 + 92 X^2 93 E3 94 / 95 ISG X 96 RCL 00 97 LASTX 98 X^2 99 / 100 + 101 REGMOVE 102 GTO 01 103 LBL 05 104 3 105 RCL 00 106 ST+ Y 107 ST* Y 108 STO M 109 CLX 110 LBL 06 111 RCL IND Y 112 RCL IND M 113 * 114 + 115 DSE Y 116 DSE M |
117 GTO 06
118 ST/ N 119 CLX 120 .1 121 % 122 X<>Y 123 ST+ Y 124 RCL 00 125 ST- Y 126 ST+ Z 127 STO M 128 LBL 07 129 CLX 130 RCL IND Z 131 RCL IND M 132 * 133 RCL N 134 * 135 ST- IND Y 136 DSE Y 137 DSE Z 138 GTO 07 139 CLX 140 RCL 00 141 ST+ Z 142 DSE M 143 GTO 07 144 GTO 12 145 END |
( 229 bytes / SIZE n2+3n+1 )
| STACK | INPUTS | OUTPUTS |
| X | / | k1 , k2 , ... |
Example: Find all the eigenvalues and the
corresponding unit eigenvectors of the same matrix
1 2 4
A = 4 3 5
7 4 7
1°) Once again, n = 3 therefore 3 STO 00
and if we choose ( 1 ; 1 ; 1 ) as initial guess, we store the 12 numbers
1 1 2 4
R01 R04 R07 R10
1 4 3 5 in
R02 R05 R08 R11 respectively.
1 7 4 7
R03 R06 R09 R12
XEQ "DFL" the HP-41 displays the successive k1-approximations ( with F02 clear )
- Eventually, k1 = 12.90692995 is in X-register and U1 is in registers R01 R02 R03 : U1 ( 0.348663346 ; 0.530674468 ; 0.772540278 )
2°) To obtain the second eigenvalue simply press R/S ( you may re-initialize R01 to Rnn but in most cases, it's not necessary )
the HP-41 sets flag F02 and displays the successive
k1-approximations converging to 12.90692995 again
then flag F02 is cleared and the HP-41 displays
the successive k2-approximations and
k2 = -2.092097594 in X-register and U2 in R01R02R03 : U2 ( 0.800454175 ; -0.041651079 ; -0.597945065 )
3°) Once again R/S ( similar displays ) and finally,
k3 = 0.185167649 in X-register and U3 in R01R02R03 : U3 ( -0.094824730 ; 0.897989404 ; -0.429678136 )
Warning:
Watch the successive k-numbers displayed by the HP-41 ( with flag
F02 set and clear ): the 2 limits must be ( approximately ) the
same.
Otherwise change the initial guess and start over again.
( It may happen that an initial vector V leads to a dominant
eigenvalue of A but a non-dominant eigenvalue of tA.
In such a case, A' will be miscalculated and the next
eigenvalues will be wrong! )
-If k1 and k2 are real but k3
... are complex, the first 2 eigenvalues and eigenvectors will
be correctly computed.
4°) The Jacobi's method
a) Symmetric ( and some non-symmetric
) Matrices
*** A non-symmetric real matrix may have complex eigenvalues, but all the eigenvalues of a real symmetric matrix A are necessarily real !
-In the Jacobi's method, the eigenvalues are the elements of a diagonal
matrix P equal to the infinite product ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....
where the Ok are rotation matrices. The eigenvectors
are the columns of O1.O2....Ok-1.Ok....
-"JCB" finds the greatest off-diagonal element ai j
( lines 22 to 46 )
-Then, O1 is determined so that it makes a pair
of off-diagonal elements zero in O1-1.A.O1
-It happens if the rotation angle µ is choosen so that
Tan 2µ = 2.ai j / ( ai i - aj j
)
-The process is repeated until the greatest off-diagonal element is
smaller than E-8 in absolute value ( line 48 )
-The successive greatest | ai j | ( i > j ) are displayed
( line 47 ) when the routine is running.
*** Though it's not really orthodox, we can try to apply the same program to non-symmetric matrices: of course, it cannot work if some eigenvalues are complex.
But if all the eigenvalues are real, the infinite product
T = [ ti j ] = ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....
is an upper triangular matrix,
the diagonal elements of T ( i-e ti i = ki
) are the n eigenvalues and the first column of U = O1.O2....Ok-1.Ok....
is an eigenvector corresponding to k1
-If all the eigenvalues are distinct, we can create an upper triangular matrix W = [ wi j ] defined by
wi j = 0 if i > j
wi i = 1
wi j = Sumk = i+1, i+2 , .... , j
ti k wk j / ( tj j - ti i )
if i < j
-Then, the n columns of U.W are the eigenvectors corresponding to the n eigenvalues t11 = k1 , ............... , tnn = kn
*** When the following program stops:
the n eigenvalues are stored in registers R01 thru
Rnn
the first eigenvector is stored in registers Rn+1 thru R2n
..................................................................................
the nth eigenvector is stored in registers
Rn2+1 thru Rn2+n as shown below.
Data Registers:
• R00 = n
• R01 = a11 .................
• Rn2-n+1 = a1n
..................................................................
( the n2 elements of A in column order )
INPUTS
• Rnn = an1 .................
• Rn2 = ann
---------
DURING
R01 thru Rn2 = the "infinite" product ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....
THE
Rn2+1 thru R2n2 = O1.O2........Ok.....
ITERATIONS R2n2+1 thru
R2n2+n: temp ( for non-symmetric matrices only )
---------
R00 = n
R01 = k1 Rn+1 = V1(1)
R2n+1 = V1(2) ....................
Rn2+1 = V1(n)
.............
OUTPUTS
Rnn = kn R2n =
Vn(1)
R3n = Vn(2) .....................
Rn2+n = Vn(n)
( n eigenvalues ; 1st eigenvector ; 2nd eigenvector ;
.................. ; nth eigenvector )
Flag: F02 Set flag F02
for a symmetric matrix
Clear flag F02 for a non-symmetric matrix
Subroutines: /
-If you don't have an HP-41CX, replace line 09 by ENTER^
SIGN CLX LBL 12 STO IND L ISG L GTO 12
CLX
-Line 47 , the HP-41 displays the successive greatest sub-diagonal
elements | aij | ( i > j )
-Line 50 is a three-byte GTO 05
-Line 132 is a three-byte GTO 01
-Lines 177 to 269 are unuseful if the matrix is symmetric
-Line 178 is a three-byte GTO 11
-Line 268 is a three-byte GTO 07
| 01 LBL "JCB"
02 RCL 00 03 X^2 04 .1 05 % 06 ST+ X 07 + 08 ISG X 09 CLRGX 10 RCL 00 11 1 12 + 13 E5 14 / 15 + 16 SIGN 17 LBL 00 18 STO IND L 19 ISG L 20 GTO 00 21 LBL 01 22 RCL 00 23 E3 24 / 25 STO M 26 2 27 + 28 ENTER^ 29 ENTER^ 30 CLX 31 LBL 02 32 RCL IND Z 33 ABS 34 X>Y? 35 STO Z 36 X>Y? 37 + 38 RDN 39 ISG Z 40 GTO 02 41 ISG M 42 RCL M 43 ST+ T 44 RDN 45 ISG Z 46 GTO 02 47 VIEW X 48 E-8 |
49 X>Y?
50 GTO 05 51 RDN 52 SIGN 53 - 54 INT 55 ENTER^ 56 STO M 57 RCL 00 58 ST/ Z 59 MOD 60 ENTER^ 61 STO N 62 LASTX 63 * 64 + 65 X<>Y 66 INT 67 RCL X 68 RCL 00 69 * 70 STO O 71 + 72 1 73 ST+ M 74 ST+ N 75 ST+ Z 76 + 77 RCL IND M 78 ST+ X 79 RCL IND Y 80 RCL IND T 81 - 82 R-P 83 CLX 84 2 85 / 86 1 87 P-R 88 X<> O 89 X<>Y 90 X<> N 91 RCL 00 92 .1 93 % 94 + 95 * 96 1 |
97 ST+ Z
98 + 99 RCL 00 100 ST- Y 101 X^2 102 ST+ Z 103 .1 104 % 105 + 106 + 107 XEQ 03 108 ST- Y 109 ST- Z 110 X^2 111 ST- Z 112 .1 113 % 114 + 115 - 116 XEQ 03 117 DSE Z 118 ST/ Z 119 / 120 INT 121 RCL 00 122 1 123 % 124 X<>Y 125 X^2 126 + 127 E3 128 / 129 ST+ Z 130 + 131 XEQ 03 132 GTO 01 133 LBL 03 134 STO M 135 LBL 04 136 CLX 137 RCL IND M 138 RCL O 139 * 140 RCL IND Y 141 RCL N 142 * 143 - 144 X<> IND M |
145 RCL N
146 * 147 X<> IND Y 148 RCL O 149 * 150 ST+ IND Y 151 ISG Y 152 CLX 153 ISG M 154 GTO 04 155 X<> M 156 RCL 00 157 RTN 158 LBL 05 159 CLD 160 RCL 00 161 X^2 162 LASTX 163 1 164 + 165 E5 166 / 167 + 168 STO O 169 RCL 00 170 LBL 06 171 RCL IND Y 172 STO IND Y 173 DSE Y 174 RDN 175 DSE Y 176 GTO 06 177 FS? 02 178 GTO 11 179 RCL 00 180 DSE X 181 E3 182 / 183 ISG X 184 STO M 185 LBL 07 186 RCL 00 187 ENTER^ 188 X^2 189 ST+ X 190 + 191 STO Y 192 RCL M |
193 INT
194 - 195 E3 196 / 197 + 198 STO N 199 RCL O 200 RCL M 201 INT 202 - 203 E-5 204 - 205 1 206 STO IND Z 207 CLX 208 LBL 08 209 RCL IND Y 210 RCL IND N 211 * 212 + 213 DSE Y 214 DSE N 215 GTO 08 216 RCL IND O 217 RCL IND Z 218 - 219 / 220 STO IND N 221 ISG M 222 GTO 07 223 RCL 00 224 DSE M 225 LBL 09 226 RCL 00 227 X^2 228 LASTX 229 RCL M 230 INT 231 * 232 + 233 E3 234 / 235 RCL 00 236 X^2 237 + 238 + 239 RCL 00 240 E5 |
241 /
242 + 243 STO P 244 CLX 245 RCL N 246 STO Q 247 CLX 248 LBL 10 249 RCL IND P 250 RCL IND Q 251 * 252 + 253 ISG Q 254 CLX 255 ISG P 256 GTO 10 257 ST+ IND P 258 X<>Y 259 DSE X 260 GTO 09 261 RCL M 262 FRC 263 E-3 264 - 265 STO M 266 DSE O 267 ISG M 268 GTO 07 269 LBL 11 270 RCL 00 271 X^2 272 LASTX 273 1 274 ST+ Z 275 + 276 E3 277 / 278 + 279 RCL 00 280 X^2 281 E6 282 / 283 + 284 REGMOVE 285 RCL 01 286 CLA 287 END |
( 429 bytes / SIZE 2n2+1 if A is symmetric
/ SIZE 2n2+n+1 if A is not symmetric )
| STACK | INPUTS | OUTPUTS |
| X | / | k1 |
Example1: Let's compute all the eigenvalues and
the eigenvectors of the matrix
1
2 4
A = 2 7 3
4 3 9
n = 3 therefore 3 STO 00 and we store the 9 numbers
1 2 4
R01 R04 R07
2 7 3
into R02 R05 R08
respectively.
4 3 9
R03 R06 R09
SF 02 ( A is symmetric ) XEQ "JCB" yields ( in 1m33s )
k1 = 12.81993499 in R01
k2 = 4.910741214 in R02
k3 = -0.730676199 in R03
and the 3 corresponding eigenvectors:
V1 ( 0.351369026 ; 0.521535689 ; 0.777521917
) in ( R04 ; R05 ; R06 )
V2 ( -0.101146468 ; 0.846760701 ; -0.522269766
) in ( R07 ; R08 ; R09 )
V3 ( -0.930757326 ; 0.104865823 ; 0.350276976
) in ( R10 ; R11 ; R12 )
Example2:
1 2 4
A = 4 3 5
7 4 7
n = 3 STO 00 and we store the 9 numbers
1 2 4
R01 R04 R07
4 3 5
into R02 R05 R08
respectively.
7 4 7
R03 R06 R09
CF 02 ( A is not symmetric )
XEQ "JCB" >>>> k1 = 12.90692994 ( in 2mn44s ) and we have
k1 =
12.90692994 in R01
k2 =
0.185167649 in R02
k3 = -2.092097593
in R03
and the 3 corresponding eigenvectors:
V1
( 0.348663347 ; 0.530674467 ; 0.772540277
) in ( R04 ; R05 ; R06 )
V2
( -0.095420104 ; 0.903627529 ; -0.432375917 )
in ( R07 ; R08 ; R09 )
V3
( -0.830063031 ; 0.043191757 ; 0.620063097 )
in ( R10 ; R11 ; R12 )
-The eigenvectors are not unit vectors except the first
one.
-Divide V2 & V3 by || V2
|| & || V3 || respectively if you want unit eigenvectors.
Notes:
-If A is non-symmetric, this program only works if all the eigenvalues
are real and distinct.
-If all the eigenvalues are real but not distinct, they will be correctly
computed but not the eigenvectors ( DATA ERROR line 219 )
-The former version of this program used the subroutine "M*" to multiply
matrices by the successive "rotation" matrices, but it was prohibitive!
-This new "JCB" uses a much simpler ( and faster ) way to perform these
products: in fact, at each step, only 2 columns and 2 rows are modified
in Oi-1.M.Oi
and only 2 columns in M.Oi .
-Execution time = 13 seconds/iteration if n = 3
-Execution time = 16 seconds/iteration if n = 4
-Execution time = 27 seconds/iteration if n = 7
-Unfortunately, the number of required iterations may increase very
much with n, especially if A is not symmetric.
Improvement?
-Actually, the Jacobi's rotations do not zero the element in position
( i , j ) in Ok-1.A.Ok if the matrix
is non-symmetric!
-We could use the formula of paragraph 5 below. Unfortunately, it involves
complex arithmetic if the radicand is negative.
-But we can replace this radicand by 0 in this case!
-If you want to use this hybrid method, replace lines 77 thru 87 by
STO P
RCL M
ST+ X RCL IND
P X^2
SQRT SIGN
X<>Y
-
STO Z RCL IND
Q +
+
P-R
STO Q
RCL IND X *
-
X<0? R-P
+
RCL IND M ST+ X
STO Z
CLX CLX
( 302 lines / 451 bytes )
-Numerical results are mixed: for the matrix of example 2 above, the
problem is solved in 2mn14s instead of 2mn44s
but in other examples, the number of iterations may be greater
and execution time longer!
-So I let you decide whether it's worthwhile or not!
Exercise: Find all the eigenvalues and eigenvectors of the 4x4 matrix
1 2 4 7
2 3 7 1
4 7 2 4
7 1 4 9
Answer: k1 = 16.97583168 , k2 = 6.365547530 , k3 = -3.301311094 , k4 = -5.040068160
V1 ( 0.455772321
; 0.346041152 ; 0.464075961 ;
0.676136537 )
V2 ( -0.142731960
; 0.681492880 ; 0.448494335 ; -0.560399745
)
V3 ( 0.842568185
; -0.247658954 ; 0.050153515 ; -0.475634862
)
V4 ( -0.248953877
; -0.595388965 ; 0.762214511 ; -0.050625961
)
-These results are obtained in 4mn22s.
b) Symmetric Matrices
only
-If the matrix is symmetric, we can save data registers: it is unuseful
to store ai j and aj i
( i # j ) in 2 registers since they are equal!
-So the program is approximately 20% faster than "JCB" and it can be
used with a matrix of order 14
Data Registers: • R00 = n ( All these registers are to be initialized before executing "JCB2" )
• R01 = a11 • Rn+1 = a1,2
• R2n = a2,3 ..................................
• R(n2+n)/2 = an-1,n
• R02 = a22 • Rn+2 = a1,3
• R2n+1 = a2,4 ...................
INPUTS
............... .................
...................
............... .................
• R3n-2 = a2,n
................ • R2n-1 = a1,n
• Rnn = ann
-----------------------------------------------------------------------------------------------------------------------------------------------------------------
When
the program stops:
R00 = n
R01 = k1 Rn+1 = V1(1)
R2n+1 = V1(2) ....................
Rn2+1 = V1(n)
.............
OUTPUTS
Rnn = kn R2n =
Vn(1)
R3n = Vn(2) .....................
Rn2+n = Vn(n)
( n eigenvalues ; 1st eigenvector ; 2nd eigenvector ; .................. ; nth eigenvector )
Flags: /
Subroutines: /
-If you don't have an HP-41CX, replace line 17 by ENTER^
SIGN CLX LBL 12 STO IND L ISG L GTO 12
CLX
-Line 56 is a three-byte GTO 10
-Line 197 is a three-byte GTO 01
| 01 LBL "JCB2"
02 RCL 00 03 X^2 04 STO Y 05 3 06 * 07 RCL 00 08 ST+ Z 09 + 10 2 11 ST/ Z 12 E3 13 * 14 / 15 + 16 ISG X 17 CLRGX 18 RCL 00 19 1 20 + 21 E5 22 / 23 + 24 SIGN 25 LBL 08 26 STO IND L 27 ISG L 28 GTO 08 29 LBL 01 30 RCL 00 31 ENTER^ 32 X^2 33 + 34 2 35 / 36 RCL 00 37 E3 38 / 39 + 40 ENTER^ 41 ENTER^ 42 CLX |
43 LBL 02
44 RCL IND Z 45 ABS 46 X>Y? 47 STO Z 48 X>Y? 49 + 50 RDN 51 DSE Z 52 GTO 02 53 VIEW X 54 E-8 55 X>Y? 56 GTO 10 57 CLX 58 ENTER^ 59 R^ 60 INT 61 STO M 62 RCL 00 63 - 64 X<>Y 65 LBL 03 66 ISG Z 67 CLX 68 RCL 00 69 + 70 R^ 71 - 72 X<Y? 73 GTO 03 74 - 75 RCL 00 76 + 77 RCL IND Y 78 RCL IND Y 79 - 80 RCL IND M 81 ST+ X 82 X<>Y 83 R-P 84 CLX |
85 2
86 / 87 1 88 P-R 89 STO N 90 RDN 91 STO O 92 X^2 93 RCL IND Y 94 * 95 STO P 96 CLX 97 RCL IND Z 98 RCL N 99 X^2 100 * 101 ST+ P 102 CLX 103 X<> IND M 104 ST+ X 105 RCL N 106 * 107 RCL O 108 * 109 ENTER^ 110 X<> P 111 + 112 X<> IND Z 113 RCL O 114 X^2 115 * 116 RCL P 117 - 118 X<> IND Y 119 RCL N 120 X^2 121 * 122 ST+ IND Y 123 CLX 124 RCL 00 125 ENTER^ 126 X^2 |
127 +
128 2 129 / 130 X<>Y 131 STO Q 132 R^ 133 STO P 134 ENTER^ 135 SIGN 136 ST- Z 137 - 138 RCL 00 139 ST* Z 140 * 141 ST+ T 142 RDN 143 + 144 RCL X 145 RCL 00 146 + 147 PI 148 INT 149 10^X 150 / 151 + 152 ISG X 153 ISG Y 154 LBL 09 155 XEQ 07 156 ISG Y 157 CLX 158 ISG X 159 GTO 09 160 RCL 00 161 DSE X 162 X<> P 163 RCL Q 164 LBL 04 165 RCL P 166 ST+ Z 167 + 168 RCL M |
169 X=Y?
170 GTO 00 171 RDN 172 XEQ 07 173 DSE P 174 GTO 04 175 LBL 00 176 RDN 177 LBL 05 178 RCL P 179 ENTER^ 180 SIGN 181 ST+ T 182 - 183 + 184 X<>Y 185 RCL M 186 X=Y? 187 GTO 00 188 X<> Z 189 XEQ 07 190 DSE P 191 GTO 05 192 LBL 00 193 X<> Z 194 LBL 06 195 DSE P 196 X=0? 197 GTO 01 198 ISG X 199 CLX 200 ISG Y 201 CLX 202 XEQ 07 203 GTO 06 204 LBL 07 205 RCL IND Y 206 RCL O 207 * 208 SIGN 209 CLX 210 RCL IND Y |
211 RCL N
212 ST* Y 213 X<> L 214 - 215 X<> IND Y 216 RCL O 217 * 218 X<> IND Z 219 RCL N 220 * 221 ST+ IND Z 222 RDN 223 RTN 224 LBL 10 225 RCL 00 226 ENTER^ 227 X^2 228 + 229 2 230 / 231 RCL 00 232 1 233 ST+ Z 234 + 235 E3 236 / 237 + 238 RCL 00 239 X^2 240 E6 241 / 242 + 243 REGMOVE 244 RCL 01 245 CLA 246 CLD 247 END |
( 358 bytes / SIZE (3n2+n+2)/2
)
| STACK | INPUTS | OUTPUTS |
| X | / | k1 |
Example:
1 2 4
A = 2 7 3
4 3 9
n = 3 therefore 3 STO 00 and we store the 6 numbers
1
R01
7 2
into R02 R04
respectively.
9 4 3
R03 R05 R06
>>> XEQ "JCB2" yields ( in 1m12s )
k1 = 12.81993499 in R01
k2 = 4.910741213 in R02
k3 = -0.730676198 in R03
and the 3 corresponding eigenvectors:
V1
( 0.351369026 ; 0.521535689 ; 0.777521917
) in ( R04 ; R05 ; R06 )
V2
( -0.101146468 ; 0.846760701 ; -0.522269766 )
in ( R07 ; R08 ; R09 )
V3
( -0.930757326 ; 0.104865823 ; 0.350276976 )
in ( R10 ; R11 ; R12 )
Notes:
-Here, the coefficients are stored differently.
-Store in successive registers starting with R01:
-the diagonal elements a11 a22
........ ann
-then a12 a13
a14 a15 ......... a1n
-then a23 a24
a25 .......... a2n
............ and so on ......
-until an-1,n
-Execution time = 10 seconds/iteration if n = 3
-Execution time = 12 seconds/iteration if n = 4
-Execution time = 22 seconds/iteration if n = 7
-For the matrix:
1 2 4 7
A = 2 3 7 1
execution time = 3mn25s ( instead of 4mn22s with "JCB"
)
4 7 2 4
7 1 4 9
-This program is very similar to "QUAD" ( cf "Quadratic Hypersurfaces
for the HP-41" )
5°) Complex Matrices
-This program uses a variant of the Jacobi's algorithm:
*** The eigenvalues of A are the diagonal-elements of an upper triangular
matrix T equal to the infinite product ...Uk-1.Uk-1-1....U2-1.U1-1.A.U1.U2....Uk-1.Uk....
where the Uk are unitary matrices. The eigenvectors
are the columns of U = U1.U2....Uk-1.Uk....
if A is Hermitian ( i-e if A equals its transconjugate )
Actually, T is diagonal if A is Hermitian.
-"JCB" finds the greatest element ai j below the main diagonal
( lines 23 to 60 )
-Then, U1 is determined so that it places a
zero in position ( i , j ) in U1-1.A.U1
U1 has the same elements as the Identity matrix, except that ui i = uj j = x and ui j = y + i.z , uj i = -y + i.z
with y + i.z = C.x , x = ( 1 + | C |2 ) -1/2 and C = 2.ai j / [ ai i - ai j + ( ( ai i - ai j )2 + 4 ai j aj i )1/2 ] ( line 156 R-P avoids a test if the denominator = 0 )
-The process is repeated until the greatest sub-diagonal element is smaller than E-8 in magnitude ( line 62 )
-The successive greatest | ai j | ( i > j ) are displayed ( line 61 ) when the routine is running.
*** If the matrix is non-Hermitian, and if all the eigenvalues are distinct, we define an upper triangular matrix W = [ wi j ] by
wi j = 0 if i > j
wi i = 1
wi j = Sumk = i+1, i+2 , .... , j
ti k wk j / ( tj j - ti i )
if i < j
-Then, the n columns of U.W are the eigenvectors corresponding
to the n eigenvalues t11 = k1 , ...............
, tnn = kn
Data Registers:
• R00 = n
• R01,R02 = a11 .................
• R2n2-2n+1,R2n2-2n+2 = a1n
.................................................................................................
( the n2 elements of A in column order )
INPUTS
• R2n-1,R2n = an1 .................
• R2n2-1,R2n2 = ann
---------
odd registers = real part of the coefficients , even registers =
imaginary part of the coefficients
DURING
R01 thru R2n2 = the "infinite" product ...Uk-1.Uk-1-1....U2-1.U1-1.A.U1.U2....Uk-1.Uk....
THE
R2n2+1 thru R4n2 = U1.U2........Uk.....
ITERATIONS R4n2+1 thru
R4n2+2n: temp ( for non-Hermitian matrices only )
---------
R00 = n
R01,R02 = k1
R2n+1,R2n+2 = V1(1)
.................... R2n2+1,R2n2+2
= V1(n)
.......................
........................
.................... ..................................
OUTPUTS
R2n-1,R2n = kn R4n-1,R4n
= Vn(1)
..................... R2n2+2n-1,R2n2+2n
= Vn(n)
( n eigenvalues ;
1st eigenvector ; ..................................
; nth eigenvector )
Flag: F02 Set flag F02
for a Hermitian matrix
Clear flag F02 for a non-Hermitian matrix
Subroutines: /
-If you don't have an HP-41CX, replace line 10 by ENTER^
SIGN CLX LBL 13 STO IND L ISG L GTO 13
CLX
-Line 64 is a three-byte GTO 06
-Line 237 is a three-byte GTO 01
-Line 332 is a three-byte GTO 11
-Lines 333 thru 471 ( and lines 312 thru 328 ) are unuseful if the
matrix is Hermitian
-There will be a DATA ERROR line 406 if all the eigenvalues are not
distinct
-Line 471 is a three-byte GTO 07
| 01 LBL "JCBZ"
02 RCL 00 03 X^2 04 ST+ X 05 .1 06 % 07 ST+ X 08 + 09 ISG X 10 CLRGX 11 RCL 00 12 1 13 + 14 ST+ X 15 E5 16 / 17 + 18 SIGN 19 LBL 00 20 STO IND L 21 ISG L 22 GTO 00 23 LBL 01 24 RCL 00 25 ST+ X 26 E3 27 / 28 ISG X 29 STO M 30 2 31 + 32 ENTER^ 33 ENTER^ 34 CLX 35 LBL 02 36 RCL IND Z 37 X^2 38 SIGN 39 ST+ T 40 CLX 41 RCL IND T 42 ST* X 43 ST+ L 44 X<> L 45 X>Y? 46 STO Z 47 X>Y? 48 + 49 RDN 50 ISG Z 51 GTO 02 52 2 53 ST+ M 54 CLX 55 RCL M 56 ST+ T 57 RDN 58 ISG Z 59 GTO 02 60 SQRT 61 VIEW X 62 E-8 63 X>Y? 64 GTO 06 65 X<> Z 66 INT 67 STO M 68 2 69 ST- Y 70 / 71 STO Y 72 RCL 00 73 ST/ Z 74 MOD 75 ST+ X 76 RCL X 77 LASTX 78 * 79 + 80 X<>Y 81 INT 82 ST+ X 83 RCL X 84 RCL 00 85 * 86 + |
87 2
88 ST+ Z 89 + 90 STO N 91 X<>Y 92 STO O 93 + 94 RCL M 95 - 96 RCL IND X 97 DSE Y 98 RCL IND Y 99 R-P 100 4 101 * 102 RCL IND M 103 RCL M 104 DSE X 105 RDN 106 RCL IND T 107 R-P 108 X<>Y 109 ST+ T 110 RDN 111 * 112 P-R 113 RCL IND N 114 RCL IND O 115 - 116 DSE N 117 DSE O 118 RCL IND N 119 STO N 120 CLX 121 RCL IND O 122 ST- N 123 RDN 124 STO O 125 RCL N 126 R-P 127 X^2 128 X<>Y 129 ST+ X 130 X<>Y 131 P-R 132 X<>Y 133 ST+ T 134 RDN 135 + 136 R-P 137 SQRT 138 X<>Y 139 2 140 / 141 X<>Y 142 P-R 143 RCL O 144 ST+ Z 145 X<> N 146 + 147 R-P 148 RCL IND M 149 DSE M 150 RCL IND M 151 R-P 152 ST+ X 153 R^ 154 ST- Z 155 X<> T 156 R-P 157 CLX 158 SIGN 159 ST- M 160 P-R 161 STO O 162 RDN 163 P-R 164 STO N 165 X<>Y 166 X<> M 167 2 168 / 169 STO Y 170 RCL 00 171 ST/ Z 172 MOD |
173 X<>Y
174 INT 175 RCL 00 176 ST+ X 177 ST* Y 178 ST* Z 179 + 180 .1 181 % 182 + 183 RCL 00 184 ST+ X 185 - 186 1 187 ST+ Z 188 + 189 RCL 00 190 X^2 191 ST+ X 192 ST+ Z 193 .1 194 % 195 + 196 + 197 XEQ 03 198 RCL 00 199 ST+ X 200 ST- Z 201 - 202 RCL 00 203 X^2 204 ST+ X 205 ST- Z 206 .1 207 % 208 + 209 - 210 XEQ 03 211 RCL 00 212 ST/ Z 213 / 214 INT 215 X<>Y 216 INT 217 RCL 00 218 X^2 219 ST+ X 220 E3 221 / 222 + 223 RCL 00 224 ST+ X 225 1 226 CHS 227 ST* M 228 ST* N 229 ST+ Z 230 ST+ T 231 + 232 E5 233 / 234 ST+ Z 235 + 236 XEQ 03 237 GTO 01 238 LBL 03 239 STO P 240 LBL 04 241 CLX 242 RCL IND P 243 RCL O 244 * 245 RCL IND Y 246 RCL N 247 * 248 + 249 PI 250 SIGN 251 ST+ T 252 CLX 253 RCL IND T 254 RCL M 255 * 256 - 257 STO Q 258 CLX |
259 RCL IND Y
260 RCL O 261 * 262 RCL IND P 263 RCL N 264 * 265 - 266 PI 267 SIGN 268 ST+ P 269 CLX 270 RCL IND P 271 RCL M 272 * 273 - 274 X<> IND Y 275 RCL M 276 * 277 RCL IND P 278 RCL O 279 * 280 + 281 PI 282 SIGN 283 ST+ Z 284 CLX 285 RCL IND Z 286 RCL N 287 * 288 + 289 X<> IND P 290 RCL N 291 * 292 RCL IND Y 293 RCL O 294 * 295 X<>Y 296 - 297 RCL P 298 SIGN 299 ST- L 300 X<> Q 301 X<> IND L 302 RCL M 303 * 304 + 305 STO IND Y 306 ISG Y 307 CLX 308 ISG P 309 GTO 04 310 X<> P 311 RTN 312 LBL 05 313 RCL P 314 SIGN 315 ST- L 316 RCL IND P 317 RCL IND L 318 R-P 319 RCL IND 04 320 DSE 04 321 RCL IND 04 322 R-P 323 X<>Y 324 ST+ T 325 RDN 326 * 327 P-R 328 RTN 329 LBL 06 330 CLD 331 FS? 02 332 GTO 11 333 RCL 00 334 X^2 335 ST+ X 336 LASTX 337 1 338 + 339 ST+ X 340 E5 341 / 342 + 343 STO O 344 RCL 00 |
345 DSE X
346 E3 347 / 348 ISG X 349 STO M 350 LBL 07 351 RCL 00 352 ST+ X 353 ENTER^ 354 X^2 355 + 356 STO Y 357 RCL M 358 INT 359 ST+ X 360 STO 03 361 - 362 E3 363 / 364 + 365 STO 04 366 RCL O 367 RCL 03 368 - 369 2 E-5 370 - 371 STO P 372 CLX 373 STO 03 374 STO N 375 STO IND Y 376 SIGN 377 ST- Y 378 STO IND Y 379 LBL 08 380 XEQ 05 381 ST+ 03 382 X<>Y 383 ST+ N 384 DSE P 385 DSE 04 386 GTO 08 387 RCL IND O 388 RCL IND P 389 - 390 PI 391 SIGN 392 ST- 04 393 ST- O 394 ST- P 395 CLX 396 RCL IND O 397 RCL IND P 398 - 399 R-P 400 RCL N 401 RCL 03 402 R-P 403 R^ 404 ST- Z 405 X<> T 406 / 407 P-R 408 STO IND 04 409 CLX 410 SIGN 411 ST+ 04 412 ST+ O 413 X<>Y 414 STO IND 04 415 ISG M 416 GTO 07 417 RCL 00 418 STO 03 419 DSE M 420 LBL 09 421 RCL 00 422 X^2 423 ENTER^ 424 STO 04 425 ST+ 04 426 LASTX 427 ST+ 04 428 RCL M 429 INT 430 ST- 04 |
431 *
432 + 433 STO N 434 E3 435 / 436 + 437 RCL 03 438 ST+ N 439 + 440 RCL 00 441 E5 442 / 443 + 444 2 445 ST* 04 446 ST* N 447 * 448 STO P 449 LBL 10 450 XEQ 05 451 RCL N 452 DSE X 453 RDN 454 ST+ IND T 455 X<>Y 456 ST+ IND N 457 PI 458 INT 459 ST+ 04 460 ISG P 461 GTO 10 462 DSE 03 463 GTO 09 464 RCL M 465 FRC 466 E-3 467 - 468 STO M 469 DSE O 470 ISG M 471 GTO 07 472 LBL 11 473 RCL 00 474 X^2 475 ST+ X 476 STO M 477 LASTX 478 ST+ X 479 STO Z 480 1 481 + 482 E5 483 / 484 + 485 LBL 12 486 RCL IND X 487 STO IND Z 488 CLX 489 SIGN 490 ST- Z 491 - 492 RCL IND X 493 STO IND Z 494 DSE Y 495 RDN 496 DSE Y 497 GTO 12 498 CLX 499 RCL M 500 R^ 501 1 502 ST+ Z 503 + 504 E3 505 / 506 + 507 RCL M 508 E6 509 / 510 + 511 REGMOVE 512 RCL 02 513 RCL 01 514 CLA 515 END |
( 780 bytes / SIZE 4n2+1 if A
is Hermitian / SIZE 4n2+2n+1 if A is non-Hermitian )
| STACK | INPUTS | OUTPUTS |
| Y | / | y1 |
| X | / | x1 |
where k1 = x1 + y1 = the first eigenvalue.
Example1:
1 4 -7.i 3-4.i
A = 4+7.i 6
1-5.i
A is equal to its transconjugate so A is Hermitian
3+4.i
1+5.i 7
1 0 4 -7 3 -4
R01 R02 R07 R08 R13 R14
Store 4 7 6 0
1 -5 into R03 R04 R09 R10
R15 R16 respectively and
n = 3 STO 00
3 4 1 5 7 0
R05 R06 R11 R12 R17 R18
SF 02 ( the matrix is Hermitian ) XEQ "JCBZ" yields in 5mn07s
k1
= 15.61385271 + 0.i = ( X , Y ) = ( R01 , R02 )
k2
= 5.230678474 + 0.i = ( R03 , R04 )
k3
= -6.844531162 + 0.i = ( R05 , R06 ) and the 3
corresponding eigenvectors
V1( 0.481585119 + 0.108201123
i , 0.393148652 + 0.527305194 i , -0.142958847 + 0.550739892
i ) = ( R07 R08 , R09 R10 , R11 R12 )
V2( -0.436763638 + 0.075843695
i , 0.210271354 - 0.463555612 i , -0.516799072 + 0.526598642
i ) = ( R13 R14 , R15 R16 , R17 R18 )
V3( -0.706095475 + 0.247553486
i , 0.326530385 + 0.449069516 i , 0.363126603 - 0.
i )
= ( R19 R20 , R21 R22 , R23 R24 )
-Due to roundoff errors, registers R02 R04 R06 contain very small numbers
instead of 0
but actually, the eigenvalues of a Hermitian matrix are always
real.
Example2:
1+2.i
2+5.i 4+7.i
A = 4+7.i 3+6.i 3+4.i
A is non-Hermitian
3+4.i
1+7.i 2+4.i
1 2 2 5 4 7
R01 R02 R07 R08 R13 R14
Store 4 7 3 6
3 4 into R03 R04 R09 R10
R15 R16 respectively and
n = 3 STO 00
3 4 1 7 2 4
R05 R06 R11 R12 R17 R18
CF 02 ( the matrix is non-Hermitian ) XEQ "JCBZ" produces in 6mn45s
k1
= 7.656606009 + 15.61073835 i = ( X , Y ) = ( R01 ,
R02 )
k2
= 1.661248138 - 1.507335315 i = ( R03 , R04 )
k3
= -3.317854151 - 2.103403073 i = ( R05 , R06 )
and the 3 corresponding eigenvectors
V1( 0.523491429 + 0.008555748
i , 0.650242685 - 0.046353000 i , 0.546288477 + 0.049882590
i ) = ( R07 R08 , R09 R10 , R11 R12
)
V2( -0.4777850326 - 0.196368365
i , 0.660547554 - 0.265888145 i , -0.356842635 + 0.318267519
i ) = ( R13 R14 , R15 R16 , R17 R18 )
V3( -0.567221101 - 0.574734664
i , 0.141603481 + 0.549379028 i , 0.480533312 - 0.090337847
i ) = ( R19 R20 , R21 R22 , R23 R24 )
Notes:
-If n = 4 , a typical execution time is 19 minutes for a non-Hermitian
matrix.
-If you don't want to get the eigenvectors, but you do want to compute
the eigenvalues:
>>> Delete lines 498 to 511 , line 476 , lines 331 to 472 , lines 312 to 328 , lines 189 to 209 and lines 02 to 22 ( 299 lines / 450 bytes / SIZE 2n2+1 )
-The eigenvalues of the 2 matrices above are now obtained in 3mn59s
& 4mn39s respectively.
-Since the SIZE is now 2n2+1 , you might calculate the eigenvalues
of a 12x12 complex matrix! Not very quickly however:
-For a complex matrix of order 7 , the eigenvalues are computed in
about 2h30mn, but with a good emulator ...
References:
[1] Francis Scheid "Numerical Analysis"
McGraw-Hill ISBN 07-055197-9
[2] J-M Ferrard "Mastering your HP-48"
D3I Diffusion ISBN 2-908791-04-8