Matrix Characteristic Polynomial - Reloaded.



#14

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.


#15


#16

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


#17

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

#18

Hey!! What are MY PARENTS doing in the first row of that picture!!!??????????!!!!!???????????!!!!???

:-)

Namir


#19

And why are they sitting with my kids?

#20

'Angel, do you have a rough ETA? I'm putting together a new Flash image for the 41CL.

Thanks, Monte


#21

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).

#22

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

#23

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.


#24

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


#25

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

#26

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


Possibly Related Threads...
Thread Author Replies Views Last Post
  AFTER HP-Prime update, Shift+Matrix CRASHES Joseph Ec 3 457 12-06-2013, 11:06 AM
Last Post: Joseph Ec
  HP Prime Matrix TERRIBLE bug and question uklo 19 1,288 11-25-2013, 12:10 PM
Last Post: Mic
  Polynomial program Michael Carey 12 766 11-20-2013, 08:43 PM
Last Post: CR Haeger
  HP Prime: editing a matrix Alberto Candel 6 457 11-20-2013, 06:26 PM
Last Post: Helge Gabert
  Absolute Value and Matrix BruceTTT 5 480 11-11-2013, 11:52 PM
Last Post: Walter B
  HP Prime polynomial long division bluesun08 13 875 10-30-2013, 03:29 AM
Last Post: parisse
  [HP-Prime xCAS] Review Polynomial Tools + BUGs + Request CompSystems 0 212 09-05-2013, 12:53 PM
Last Post: CompSystems
  WP-34S Matrix operations with routine-local registers? Tom Grydeland 1 271 09-04-2013, 10:46 AM
Last Post: Marcus von Cube, Germany
  Matrix Richard Berler 3 321 08-18-2013, 06:24 PM
Last Post: Paul Dale
  Advantage/CCD Matrix Challenge Ángel Martin 1 249 08-09-2013, 06:22 PM
Last Post: Thomas Klemm

Forum Jump: