The Museum of HP Calculators

HP Forum Archive 21

 Matrix Characteristic Polynomial - Reloaded.Message #1 Posted by Ángel Martin on 19 Aug 2013, 9:02 a.m. Just to keep you all updated, the SandMatrix module is close to release. While preparing the manual I’ve been tweaking and improving some of the routines, taking advantage of its extended function set. One of these has been CHRPOL, used to obtain the characteristic polynomial of a matrix. The new version code is about 1/3rd. of the original one, and does the job in about half the time – if not faster. Such is the improvement derived from the “direct approach”. An image is worth 1,000 words, so here’s the program listing for you.- Stay tuned, the real think will be packed… Magic in 60 program steps. Uses the Faddevv-Leverrier method - you can almost read the algorithm from the ALPHA register, directly using the abstraction layer of the Matrix functions. ```1 LBL "CHRPOL" MNAME in Alpha 2 DIM? n,00n 3 I#J? not square? 4 -ADV MATRX show error 5 ASTO 01 MNAME 6 -CCD MATRIX shows 'RUNNING…" 7 "|-,P" requires auxiliary array 8 MAT= P = A 9 ASWAP swap alpha around comma 10 DIM? n,00n 11 INT n 12 E 13 + n+1 14 MDET independent term 15 STO IND Y stored in Rn+1 16 ASWAP swap again 17 MAT= avoids LU issues 18 DIM? 19 "#" another auxiliary array 20 MATDIM 21 FRC 0,00n 22 2 23 + 2,00n 24 STO 00 25 CF 21 not halting VIEW 26 LBL 00 27 VIEW 00 shows current index 28 "#" 29 MIDN [#] = [I] 30 "P" 31 MTRACE tr(P) 32 RCL 00 33 INT k+1 34 E 35 - k 36 / 37 CHS 38 STO IND 00 pk = -tr(P) / k 39 "X,#,#" 40 MAT* [#] = pk [I] 41 "P,#,#" 42 MAT+ [#] = [P] + p[ I] 43 CLA 44 ARCL 01 45 "|-,#,P" 46 M*M B= A (B - p I) 47 ISG 00 increment counter 48 GTO 00 next coefficient 49 DIM? n,00n 50 FRC 0,00n 51 E 52 STO 01 it's monic (!) 53 E3/E+ 1.001 54 + 1.00(n+1) - cnt'l word 55 "#" 56 PURFL remove auxiliary arrays 57 "P" 58 PURFL 59 MNAME? bring MNAME back 60 END done ``` The polynomial control word is left in X upon termination - that is, registers containing the n+1 coefficients in the form bbb.eee -from Rbbb to Reee. Cheers, ÁM PS. Here’s an article that may be useful to jog your memory: http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=944 Edited: 19 Aug 2013, 9:10 a.m.

 Re: Matrix Characteristic Polynomial - Reloaded.Message #2 Posted by Massimo Gnerucci (Italy) on 19 Aug 2013, 10:01 a.m.,in response to message #1 by Ángel Martin

 Re: Matrix Characteristic Polynomial - Reloaded.Message #3 Posted by Ángel Martin on 19 Aug 2013, 1:23 p.m.,in response to message #2 by Massimo Gnerucci (Italy) LOL !!! You gotta love this forum :) A suitable ovation for what may well be "the last program ever to be written using the Advantage-based tool-set". Cheers, 'AM

 Re: Matrix Characteristic Polynomial - Reloaded.Message #4 Posted by Massimo Gnerucci (Italy) on 19 Aug 2013, 1:55 p.m.,in response to message #3 by Ángel Martin Quote: LOL !!! You gotta love this forum :) I really do!Besides you deserve our applause for your continued work of love on Mcode routines.Thank you Ángel.Massimo

 Re: Matrix Characteristic Polynomial - Reloaded.Message #5 Posted by Namir on 19 Aug 2013, 2:24 p.m.,in response to message #2 by Massimo Gnerucci (Italy) Hey!! What are MY PARENTS doing in the first row of that picture!!!??????????!!!!!???????????!!!!??? :-) Namir

 Re: Matrix Characteristic Polynomial - Reloaded.Message #6 Posted by Jim Horn on 19 Aug 2013, 3:32 p.m.,in response to message #5 by Namir And why are they sitting with my kids?

 Re: Matrix Characteristic Polynomial - Reloaded.Message #7 Posted by Monte Dalrymple on 19 Aug 2013, 4:54 p.m.,in response to message #1 by Ángel Martin 'Angel, do you have a rough ETA? I'm putting together a new Flash image for the 41CL.Thanks, Monte

 Re: Matrix Characteristic Polynomial - Reloaded.Message #8 Posted by Ángel Martin on 20 Aug 2013, 9:29 a.m.,in response to message #7 by Monte Dalrymple I'm hoping for end of this week, just a couple of wrinkles left. I'll send you the ROM images as soon as they're ready, won't wait for the manual (which is another thick one).

 Re: Matrix Characteristic Polynomial - Reloaded.Message #9 Posted by Robert Prosperi on 19 Aug 2013, 8:18 p.m.,in response to message #1 by Ángel Martin Angel: Quote: Uses the Faddevv-Leverrier method Wow, a blast from my very distant past! In college, I had a Mech Eng & Numerical Analysis Professor that seemed to be able to prove virtually anything he chose, using Feddevv-Leverrier. Always mentioned that Faddevv was married to Faddevva (sp?) and that they worked on these matrix equations together. Such Marital bliss... Were it not for your post I may have believed I imagined it all. Thanks for the mention. It's sometimes odd how these posts will hit you. --bob prosperi

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.

 Re: Matrix Characteristic Polynomial - Reloaded.Message #11 Posted by Ángel Martin on 22 Aug 2013, 2:14 p.m.,in response to message #10 by Thomas Klemm Nice tricks Thomas! Would you also have one to create an identity matrix - the quickest possible way? Also if you're going to end up using the SOLVEr, might as well have started there, solving for x in the equation: Det( A - x I) = 0 Cheers, 'AM

 Re: Matrix Characteristic Polynomial - Reloaded.Message #12 Posted by Thomas Klemm on 22 Aug 2013, 4:48 p.m.,in response to message #11 by Ángel Martin If we know that A is invertible: ```RESULT E RCL MATRIX A ENTER / ``` Otherwise something along this: ``` DIM(n+1, n) AT DIM(n, n) |1 1 1| |1 0 0 0| |1 0 0| |1 1 1| -> |0 0 0| -> |1 0 0 0| -> |0 1 0| |0 0 0| |1 0 0 0| |0 0 1| |0 0 0| ``` Kind regards Thomas

 Re: Matrix Characteristic Polynomial - Reloaded.Message #13 Posted by Thomas Klemm on 22 Aug 2013, 5:33 p.m.,in response to message #11 by Ángel Martin Quote: Det( A - x I) = 0 That's actually a nice idea: ```LBL E RESULT B RCL MATRIX E * RCL MATRIX A - MATRIX 9 RTN 0 SOLVE E -> 10 ``` Cheers Thomas

Go back to the main exhibit hall