Re: Matrix Characteristic Polynomial - Reloaded. Message #10 Posted by Thomas Klemm on 21 Aug 2013, 8:24 a.m., in response to message #1 by Ángel Martin
For those interested here's a HP-15C version of the algorithm:
Listing
001 - 42,21,11 LBL A
002 - 45,23,11 RCL DIM A ; -- n n
003 - 1 1 ; -- n n 1
004 - 42,23,14 DIM D ; D(n, 1)
005 - 44,16,14 STO MATRIX D ; D = [1, ..., 1]
006 - 44 0 STO 0 ; a0 = 1
007 - 34 x<>y ; -- n 1 n
008 - 26 EEX
009 - 3 3 ; -- n 1 n 1000
010 - 10 / ; -- n 1 0.00n
011 - 40 + ; -- n 1.00n
012 - 44 25 STO I ; I = 1.00n
013 - 45,16,11 RCL MATRIX A
014 - 44,16,12 STO MATRIX B ; B = A
015 - 42,21, 0 LBL 0
016 - 42,26,13 RESULT C
017 - 45,16,11 RCL MATRIX A
018 - 45,16,12 RCL MATRIX B
019 - 20 x ; C = A * B
020 - 45,23,12 RCL DIM B ; -- n n
021 - 1 1 ; -- n n 1
022 - 40 + ; -- n n+1
023 - 42,23,12 DIM B ; B(n, n+1)
024 - 45,16,12 RCL MATRIX B ; -- n n+1 B(n, n+1)
025 - 42,16, 4 MATRIX 4 ; -- n n+1 B(n+1, n)
026 - 1 1 ; -- n n+1 B(n+1, n) 1
027 - 43 33 R^ ; -- n+1 B(n+1, n) 1 n
028 - 42,23,12 DIM B ; B(1, n)
029 - 42,26,15 RESULT E
030 - 45,16,12 RCL MATRIX B ; -- B(1, n)
031 - 45,16,14 RCL MATRIX D ; -- B(1, n) D(n, 1)
032 - 20 x ; E = B x D
033 - 42,16, 9 MATRIX 9 ; p = |E| = Tr(B)
034 - 16 CHS ; -- -p
035 - 45 25 RCL I ; -- -p k.00n
036 - 43 44 INT ; -- -p k
037 - 10 / ; -- -p/k
038 - 44 24 STO (i) ; -- ak = -p/k
039 - 42,26,12 RESULT B
040 - 45,16,11 RCL MATRIX A
041 - 20 x ; -p/k * A
042 - 45,16,13 RCL MATRIX C
043 - 40 + ; B = C -p/k * A
044 - 42, 6,25 ISG I
045 - 22 0 GTO 0
Calculating the Trace of a Matrix
DIM(n, n+1) AT DIM(1, n) DET
|* . .| |* . . .| |* * *| |1|
|. * .| -> |* . . .| -> |. . 0| -> |* * *| x |1| -> |*| -> *
|. . *| |* 0 0 0| |. . 0| |1|
|. . 0|
Example:
3
ENTER
DIM A
MATRIX 1
USER
3 STO A
1 STO A
5 STO A
3 STO A
3 STO A
1 STO A
4 STO A
6 STO A
4 STO A
GSB A
RCL 0 -> 1
RCL 1 -> -10
RCL 2 -> 4
RCL 3 -> -40
Thus the characteristic polynomial is: p(x) = x3 - 10x2 + 4x - 40.
Eigenvalues
We can use this short program to solve p(x) = 0.
LBL B
RCL+ 1
*
RCL+ 2
*
RCL+ 3
RTN
CLx
SOLVE B
-> 10.0000
The other solutions are 2i and -2i which can't be found with the built-in solver.
Cheers
Thomas
Edited: 21 Aug 2013, 4:39 p.m.
|