Matrix Characteristic Polynomial - Reloaded.

[Ángel Martin]

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.

[Massimo Gnerucci (Italy)]

[Ángel Martin]

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

[Massimo Gnerucci (Italy)]

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

[Namir]

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

:-)

Namir

[Jim Horn]

And why are they sitting with my kids?

[Monte Dalrymple]

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

Thanks, Monte

[Ángel Martin]

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

[Robert Prosperi]

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

[Thomas Klemm]

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) = x³ - 10x² + 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.

[Ángel Martin]

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

[Thomas Klemm]

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

[Thomas Klemm]

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