## EIGEN: Characteristic equation+eigenvalues for order 2,3,4,5

Programs for HP-41

### EIGEN: Characteristic equation+eigenvalues for order 2,3,4,5

EIGEN computes the characteristic equation and eigenvalues of matrices of order 2, 3, 4, and 5. It computes the coefficients of the characteristic polynomial directly, then solves for the roots of the polynomial. Routines from the MATHSTAT (or MATH I) module are used -- EIGEN expects that MATRIX is run first to enter the matrix elements, then both DET (XROM 01,06) and ROOTS (XROM 01,12) are called.

Compared to alternate approaches such as Samuelson's formula (see, e.g., Jean-Marc Baillard's program CP), computing the polynomial coefficients directly for low-order matrices has the advantage of speed at the cost of longer code, at least for order 5 matrices, and limited matrix order. Note that the method used by EIGEN is in general the *opposite* approach of an algorithm that obtains the best numerical accuracy (where the eigenvalues are computed first, then the polynomial coefficients; see Matlab's POLY function). However, the accuracy issues associated with computing polynomial coefficients and the ill-conditioning of computing polynomial roots are not a concern for very low order matrices -- compare Matlab's ROOTS(POLY(1:5)) and ROOTS(POLY(1:21)). (And computations for higher order matrices are best done on a computer, if only for data entry!)

To load EIGEN, first install the MATHSTAT (or MATH I) module and XEQ SIZE 090 to allocate sufficient program and memory registers. To run EIGEN, XEQ MATRIX (MATHSTAT/MATH I module) and enter the matrix elements as prompted. Then XEQ EIGEN. The coefficients of the characteristic polynomial will first be displayed, then their roots, the eigenvalues. R/S through each of the displays.

The code below is generated by Antonio Lagana's i41CX+ iPhone emulator.

Code: Select all
`01 LBL "EIGEN"02 RCL 1403 604 X<=Y?05 GTO 0106 X<>Y07 STO 5908 GTO IND X09 LBL 0210 LBL 2011 RCL 1512 RCL 1813 *14 RCL 1615 RCL 1716 *17 -18 STO 0019 LBL 2120 RCL 1521 RCL 1822 +23 CHS24 STO 0125 126 STO 0227 GTO 0628 LBL 0329 LBL 3130 RCL 1931 RCL 2332 +33 RCL 1534 *35 RCL 1936 RCL 2337 *38 +39 RCL 1640 RCL 1841 *42 -43 RCL 1744 RCL 2145 *46 -47 RCL 2048 RCL 2249 *50 -51 STO 3152 LBL 3253 RCL 1554 RCL 1955 +56 RCL 2357 +58 CHS59 STO 3260 LBL 3061 SF 0262 XROM 01,0663 CHS64 STO 0065 RCL 3166 STO 0167 RCL 3268 STO 0269 170 STO 0371 GTO 0672 LBL 0473 LBL 4174 RCL 2575 RCL 3076 *77 RCL 2678 RCL 2979 *80 -81 RCL 2082 *83 RCL 2684 RCL 2885 *86 RCL 2487 RCL 3088 *89 -90 RCL 2191 *92 +93 RCL 2494 RCL 2995 *96 RCL 2597 RCL 2898 *99 -100 RCL 22101 *102 +103 STO 32104 RCL 25105 RCL 30106 +107 RCL 20108 *109 RCL 25110 RCL 30111 *112 +113 RCL 21114 RCL 24115 *116 -117 RCL 22118 RCL 28119 *120 -121 RCL 26122 RCL 29123 *124 -125 STO 33126 RCL 15127 *128 ST+ 32129 RCL 25130 RCL 30131 +132 RCL 19133 *134 RCL 21135 RCL 23136 *137 -138 RCL 22139 RCL 27140 *141 -142 RCL 16143 *144 ST- 32145 RCL 20146 RCL 30147 +148 RCL 23149 *150 RCL 19151 RCL 24152 *153 -154 RCL 26155 RCL 27156 *157 -158 RCL 17159 *160 ST- 32161 RCL 20162 RCL 25163 +164 RCL 27165 *166 RCL 19167 RCL 28168 *169 -170 RCL 23171 RCL 29172 *173 -174 RCL 18175 *176 ST- 32177 RCL 32178 CHS179 STO 41180 LBL 42181 RCL 20182 RCL 25183 +184 RCL 30185 +186 RCL 15187 *188 RCL 33189 +190 RCL 16191 RCL 19192 *193 -194 RCL 17195 RCL 23196 *197 -198 RCL 18199 RCL 27200 *201 -202 STO 42203 LBL 43204 RCL 15205 RCL 20206 +207 RCL 25208 +209 RCL 30210 +211 CHS212 STO 43213 LBL 40214 SF 02215 XROM 01,06216 STO 00217 RCL 41218 STO 01219 RCL 42220 STO 02221 RCL 43222 STO 03223 1224 STO 04225 GTO 06226 LBL 05227 LBL 51228 RCL 33229 RCL 39230 +231 STO 78232 RCL 27233 *234 RCL 33235 RCL 39236 *237 +238 RCL 28239 RCL 32240 *241 -242 RCL 29243 RCL 37244 *245 -246 RCL 34247 RCL 38248 *249 -250 STO 60251 RCL 21252 *253 RCL 78254 RCL 26255 *256 RCL 28257 RCL 31258 *259 -260 RCL 29261 RCL 36262 *263 -264 STO 62265 RCL 22266 *267 -268 RCL 27269 RCL 39270 +271 STO 79272 RCL 31273 *274 RCL 26275 RCL 32276 *277 -278 RCL 34279 RCL 36280 *281 -282 STO 63283 RCL 23284 *285 -286 RCL 27287 RCL 33288 +289 STO 80290 RCL 36291 *292 RCL 26293 RCL 37294 *295 -296 RCL 31297 RCL 38298 *299 -300 STO 64301 RCL 24302 *303 -304 STO 84305 RCL 33306 RCL 39307 *308 RCL 34309 RCL 38310 *311 -312 STO 72313 RCL 27314 *315 RCL 32316 RCL 39317 *318 RCL 34319 RCL 37320 *321 -322 STO 73323 RCL 28324 *325 -326 RCL 32327 RCL 38328 *329 RCL 33330 RCL 37331 *332 -333 STO 74334 RCL 29335 *336 +337 STO 61338 RCL 84339 +340 RCL 15341 *342 STO 51343 RCL 60344 RCL 20345 *346 RCL 28347 RCL 30348 *349 RCL 29350 RCL 35351 *352 +353 STO 69354 CHS355 RCL 78356 RCL 25357 *358 +359 STO 66360 RCL 22361 *362 -363 RCL 25364 RCL 32365 *366 RCL 34367 RCL 35368 *369 +370 STO 70371 CHS372 RCL 79373 RCL 30374 *375 +376 STO 67377 RCL 23378 *379 -380 RCL 25381 RCL 37382 *383 RCL 30384 RCL 38385 *386 +387 STO 71388 CHS389 RCL 80390 RCL 35391 *392 +393 STO 68394 RCL 24395 *396 -397 RCL 16398 *399 ST- 51400 RCL 62401 RCL 20402 *403 RCL 66404 RCL 21405 *406 -407 RCL 23408 RCL 31409 *410 RCL 24411 RCL 36412 *413 +414 RCL 34415 RCL 38416 *417 +418 RCL 33419 RCL 39420 *421 -422 RCL 25423 *424 +425 RCL 23426 RCL 26427 *428 RCL 29429 RCL 38430 *431 +432 RCL 28433 RCL 39434 *435 -436 RCL 30437 *438 -439 RCL 24440 RCL 26441 *442 RCL 28443 RCL 34444 *445 +446 RCL 29447 RCL 33448 *449 -450 RCL 35451 *452 -453 RCL 17454 *455 ST+ 51456 RCL 63457 RCL 20458 *459 RCL 67460 RCL 21461 *462 -463 RCL 22464 RCL 26465 *466 RCL 24467 RCL 36468 *469 +470 RCL 29471 RCL 37472 *473 +474 RCL 27475 RCL 39476 *477 -478 RCL 30479 *480 +481 RCL 22482 RCL 31483 *484 RCL 34485 RCL 37486 *487 +488 RCL 32489 RCL 39490 *491 -492 RCL 25493 *494 -495 RCL 24496 RCL 31497 *498 RCL 29499 RCL 32500 *501 +502 RCL 27503 RCL 34504 *505 -506 RCL 35507 *508 -509 RCL 18510 *511 ST+ 51512 RCL 64513 RCL 20514 *515 RCL 68516 RCL 21517 *518 -519 RCL 22520 RCL 26521 *522 RCL 23523 RCL 31524 *525 +526 STO 81527 RCL 28528 RCL 32529 *530 +531 RCL 27532 RCL 33533 *534 -535 RCL 35536 *537 +538 RCL 22539 RCL 36540 *541 RCL 32542 RCL 38543 *544 +545 RCL 33546 RCL 37547 *548 -549 RCL 25550 *551 -552 RCL 23553 RCL 36554 *555 RCL 28556 RCL 37557 *558 +559 RCL 27560 RCL 38561 *562 -563 RCL 30564 *565 -566 RCL 19567 *568 ST+ 51569 RCL 61570 RCL 21571 *572 ST+ 51573 RCL 72574 RCL 26575 *576 RCL 31577 RCL 39578 *579 RCL 34580 RCL 36581 *582 -583 STO 75584 RCL 28585 *586 -587 RCL 31588 RCL 38589 *590 RCL 33591 RCL 36592 *593 -594 STO 76595 RCL 29596 *597 +598 RCL 22599 *600 ST- 51601 RCL 73602 RCL 26603 *604 RCL 75605 RCL 27606 *607 -608 RCL 31609 RCL 37610 *611 RCL 32612 RCL 36613 *614 -615 STO 77616 RCL 29617 *618 +619 RCL 23620 *621 ST+ 51622 RCL 74623 RCL 26624 *625 RCL 76626 RCL 27627 *628 -629 RCL 77630 RCL 28631 *632 +633 RCL 24634 *635 ST- 51636 LBL 52637 RCL 78638 RCL 27639 +640 STO 82641 RCL 21642 *643 RCL 81644 -645 RCL 24646 RCL 36647 *648 -649 STO 65650 RCL 60651 +652 RCL 15653 *654 CHS655 STO 52656 RCL 82657 RCL 20658 *659 RCL 22660 RCL 25661 *662 -663 RCL 23664 RCL 30665 *666 -667 RCL 24668 RCL 35669 *670 -671 RCL 16672 *673 ST+ 52674 RCL 78675 RCL 21676 +677 RCL 25678 *679 RCL 20680 RCL 26681 *682 -683 RCL 69684 -685 RCL 17686 *687 ST+ 52688 RCL 79689 RCL 21690 +691 RCL 30692 *693 RCL 20694 RCL 31695 *696 -697 RCL 70698 -699 RCL 18700 *701 ST+ 52702 RCL 80703 RCL 21704 +705 RCL 35706 *707 RCL 20708 RCL 36709 *710 -711 RCL 71712 -713 RCL 19714 *715 ST+ 52716 RCL 60717 RCL 21718 *719 ST- 52720 RCL 62721 RCL 22722 *723 ST+ 52724 RCL 63725 RCL 23726 *727 ST+ 52728 RCL 64729 RCL 24730 *731 ST+ 52732 RCL 61733 ST- 52734 LBL 53735 RCL 82736 RCL 21737 +738 STO 83739 RCL 15740 *741 RCL 65742 +743 RCL 60744 +745 RCL 16746 RCL 20747 *748 -749 RCL 17750 RCL 25751 *752 -753 RCL 18754 RCL 30755 *756 -757 RCL 19758 RCL 35759 *760 -761 STO 53762 LBL 54763 RCL 83764 RCL 15765 +766 CHS767 STO 54768 LBL 50769 SF 02770 XROM 01,06771 CHS772 STO 00773 RCL 51774 STO 01775 RCL 52776 STO 02777 RCL 53778 STO 03779 RCL 54780 STO 04781 1782 STO 05783 LBL 06784 RCL 59785 STO 06786 STO 22787 2788 ST+ 06789 LBL 07790 DSE 06791 GTO 08792 SF 00793 SF 21794 XROM 01,12795 RTN796 LBL 08797 RCL 06798 1799 -800 "a"801 FIX 00802 CF 29803 ARCL X804 >"="805 FIX 04806 SF 29807 ARCL IND X808 PROMPT809 GTO 07810 LBL 00811 LBL 01812 "N\1D2,3,4,5"813 AVIEW814 PSE815 "XEQ MATH:MATRIX"816 >" 1ST"817 AVIEW818 END`
eigen
...

Posts: 7
Joined: Sat Jul 18, 2009 6:24 pm

### Re: EIGEN: Characteristic equation+eigenvalues for order 2,3,4,5

I was curious how this code compared to SANDMATH'-5 packages's compiled M-code "CHRPOL", which I believe has the identical functionality. Using the Wilkinson test matrix,

Code: Select all
`>> W5 = wilkinson(5)W5 =     2     1     0     0     0     1     1     1     0     0     0     1     0     1     0     0     0     1     1     1     0     0     0     1     2`

which has characteristic polynomial and eigenvalues (double precision),

Code: Select all
`>> poly(W5), eig(W5).'ans =   1.000000000000000  -6.000000000000000   8.999999999999996   4.000000000000008 -13.000000000000007   4.000000000000000ans =    -1.114907541476756   0.381966011250105   1.254101688365053   2.618033988749895   2.860805853111704`

This code return the eigenvalues

Code: Select all
`-1.114907542   0.3819660170   1.254101674   2.618034045   2.860805070`

(maximum error 8e-7) after about 37 s to compute the characteristic polynomial, and 115 s for MATH:ROOTS to compute its roots (using the i41CX's ridiculously slow "Default Speed").

SANDMATH-5:CHRPOL returns the correct characteristic polynomial after 7+ minutes. However, CHRPOL doesn't return the correct eigenvalues if you hit R/S after the characteristic polynomial coefficients are shown. Perhaps I have not called CHRPOL correctly. A call to SANDMATH-5:POLYN could be used to compute eigenvalues (after reentering the coefficents), but I stopped waiting to verify results after another 9+ more minutes and cranking up the i41CX's speed to the max.

I believe that the low order cases (N = 2, 3, 4, 5) may be sped up considerably by using specialized code, as done above.
eigen
...

Posts: 7
Joined: Sat Jul 18, 2009 6:24 pm

### Who is online

Users browsing this forum: Google [Bot] and 1 guest