HP Articles Forum
[Return to the Index ]
[ Previous | Next ]
EIGEN: Characteristic equation+eigenvalues for order 2,3,4,5
Posted by Steven Thomas Smith on 1 Sept 2009, 12:39 p.m.
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. It would be straightforward to edit the register numbers to use the (more efficient) matrix storage routines from the Advantage PAC.
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.
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,
>> 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),
>> poly(W5), eig(W5).' ans = 1.000000000000000 -6.000000000000000 8.999999999999996 4.000000000000008 -13.000000000000007 4.000000000000000 ans = -1.114907541476756 0.381966011250105 1.254101688365053 2.618033988749895 2.860805853111704
This code return the eigenvalues
-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 here.
The code below is generated by Antonio Lagana's i41CX+ iPhone emulator.
01 LBL "EIGEN" 02 RCL 14 03 6 04 X<=Y? 05 GTO 01 06 X<>Y 07 STO 59 08 GTO IND X 09 LBL 02 10 LBL 20 11 RCL 15 12 RCL 18 13 * 14 RCL 16 15 RCL 17 16 * 17 - 18 STO 00 19 LBL 21 20 RCL 15 21 RCL 18 22 + 23 CHS 24 STO 01 25 1 26 STO 02 27 GTO 06 28 LBL 03 29 LBL 31 30 RCL 19 31 RCL 23 32 + 33 RCL 15 34 * 35 RCL 19 36 RCL 23 37 * 38 + 39 RCL 16 40 RCL 18 41 * 42 - 43 RCL 17 44 RCL 21 45 * 46 - 47 RCL 20 48 RCL 22 49 * 50 - 51 STO 31 52 LBL 32 53 RCL 15 54 RCL 19 55 + 56 RCL 23 57 + 58 CHS 59 STO 32 60 LBL 30 61 SF 02 62 XROM 01,06 63 CHS 64 STO 00 65 RCL 31 66 STO 01 67 RCL 32 68 STO 02 69 1 70 STO 03 71 GTO 06 72 LBL 04 73 LBL 41 74 RCL 25 75 RCL 30 76 * 77 RCL 26 78 RCL 29 79 * 80 - 81 RCL 20 82 * 83 RCL 26 84 RCL 28 85 * 86 RCL 24 87 RCL 30 88 * 89 - 90 RCL 21 91 * 92 + 93 RCL 24 94 RCL 29 95 * 96 RCL 25 97 RCL 28 98 * 99 - 100 RCL 22 101 * 102 + 103 STO 32 104 RCL 25 105 RCL 30 106 + 107 RCL 20 108 * 109 RCL 25 110 RCL 30 111 * 112 + 113 RCL 21 114 RCL 24 115 * 116 - 117 RCL 22 118 RCL 28 119 * 120 - 121 RCL 26 122 RCL 29 123 * 124 - 125 STO 33 126 RCL 15 127 * 128 ST+ 32 129 RCL 25 130 RCL 30 131 + 132 RCL 19 133 * 134 RCL 21 135 RCL 23 136 * 137 - 138 RCL 22 139 RCL 27 140 * 141 - 142 RCL 16 143 * 144 ST- 32 145 RCL 20 146 RCL 30 147 + 148 RCL 23 149 * 150 RCL 19 151 RCL 24 152 * 153 - 154 RCL 26 155 RCL 27 156 * 157 - 158 RCL 17 159 * 160 ST- 32 161 RCL 20 162 RCL 25 163 + 164 RCL 27 165 * 166 RCL 19 167 RCL 28 168 * 169 - 170 RCL 23 171 RCL 29 172 * 173 - 174 RCL 18 175 * 176 ST- 32 177 RCL 32 178 CHS 179 STO 41 180 LBL 42 181 RCL 20 182 RCL 25 183 + 184 RCL 30 185 + 186 RCL 15 187 * 188 RCL 33 189 + 190 RCL 16 191 RCL 19 192 * 193 - 194 RCL 17 195 RCL 23 196 * 197 - 198 RCL 18 199 RCL 27 200 * 201 - 202 STO 42 203 LBL 43 204 RCL 15 205 RCL 20 206 + 207 RCL 25 208 + 209 RCL 30 210 + 211 CHS 212 STO 43 213 LBL 40 214 SF 02 215 XROM 01,06 216 STO 00 217 RCL 41 218 STO 01 219 RCL 42 220 STO 02 221 RCL 43 222 STO 03 223 1 224 STO 04 225 GTO 06 226 LBL 05 227 LBL 51 228 RCL 33 229 RCL 39 230 + 231 STO 78 232 RCL 27 233 * 234 RCL 33 235 RCL 39 236 * 237 + 238 RCL 28 239 RCL 32 240 * 241 - 242 RCL 29 243 RCL 37 244 * 245 - 246 RCL 34 247 RCL 38 248 * 249 - 250 STO 60 251 RCL 21 252 * 253 RCL 78 254 RCL 26 255 * 256 RCL 28 257 RCL 31 258 * 259 - 260 RCL 29 261 RCL 36 262 * 263 - 264 STO 62 265 RCL 22 266 * 267 - 268 RCL 27 269 RCL 39 270 + 271 STO 79 272 RCL 31 273 * 274 RCL 26 275 RCL 32 276 * 277 - 278 RCL 34 279 RCL 36 280 * 281 - 282 STO 63 283 RCL 23 284 * 285 - 286 RCL 27 287 RCL 33 288 + 289 STO 80 290 RCL 36 291 * 292 RCL 26 293 RCL 37 294 * 295 - 296 RCL 31 297 RCL 38 298 * 299 - 300 STO 64 301 RCL 24 302 * 303 - 304 STO 84 305 RCL 33 306 RCL 39 307 * 308 RCL 34 309 RCL 38 310 * 311 - 312 STO 72 313 RCL 27 314 * 315 RCL 32 316 RCL 39 317 * 318 RCL 34 319 RCL 37 320 * 321 - 322 STO 73 323 RCL 28 324 * 325 - 326 RCL 32 327 RCL 38 328 * 329 RCL 33 330 RCL 37 331 * 332 - 333 STO 74 334 RCL 29 335 * 336 + 337 STO 61 338 RCL 84 339 + 340 RCL 15 341 * 342 STO 51 343 RCL 60 344 RCL 20 345 * 346 RCL 28 347 RCL 30 348 * 349 RCL 29 350 RCL 35 351 * 352 + 353 STO 69 354 CHS 355 RCL 78 356 RCL 25 357 * 358 + 359 STO 66 360 RCL 22 361 * 362 - 363 RCL 25 364 RCL 32 365 * 366 RCL 34 367 RCL 35 368 * 369 + 370 STO 70 371 CHS 372 RCL 79 373 RCL 30 374 * 375 + 376 STO 67 377 RCL 23 378 * 379 - 380 RCL 25 381 RCL 37 382 * 383 RCL 30 384 RCL 38 385 * 386 + 387 STO 71 388 CHS 389 RCL 80 390 RCL 35 391 * 392 + 393 STO 68 394 RCL 24 395 * 396 - 397 RCL 16 398 * 399 ST- 51 400 RCL 62 401 RCL 20 402 * 403 RCL 66 404 RCL 21 405 * 406 - 407 RCL 23 408 RCL 31 409 * 410 RCL 24 411 RCL 36 412 * 413 + 414 RCL 34 415 RCL 38 416 * 417 + 418 RCL 33 419 RCL 39 420 * 421 - 422 RCL 25 423 * 424 + 425 RCL 23 426 RCL 26 427 * 428 RCL 29 429 RCL 38 430 * 431 + 432 RCL 28 433 RCL 39 434 * 435 - 436 RCL 30 437 * 438 - 439 RCL 24 440 RCL 26 441 * 442 RCL 28 443 RCL 34 444 * 445 + 446 RCL 29 447 RCL 33 448 * 449 - 450 RCL 35 451 * 452 - 453 RCL 17 454 * 455 ST+ 51 456 RCL 63 457 RCL 20 458 * 459 RCL 67 460 RCL 21 461 * 462 - 463 RCL 22 464 RCL 26 465 * 466 RCL 24 467 RCL 36 468 * 469 + 470 RCL 29 471 RCL 37 472 * 473 + 474 RCL 27 475 RCL 39 476 * 477 - 478 RCL 30 479 * 480 + 481 RCL 22 482 RCL 31 483 * 484 RCL 34 485 RCL 37 486 * 487 + 488 RCL 32 489 RCL 39 490 * 491 - 492 RCL 25 493 * 494 - 495 RCL 24 496 RCL 31 497 * 498 RCL 29 499 RCL 32 500 * 501 + 502 RCL 27 503 RCL 34 504 * 505 - 506 RCL 35 507 * 508 - 509 RCL 18 510 * 511 ST+ 51 512 RCL 64 513 RCL 20 514 * 515 RCL 68 516 RCL 21 517 * 518 - 519 RCL 22 520 RCL 26 521 * 522 RCL 23 523 RCL 31 524 * 525 + 526 STO 81 527 RCL 28 528 RCL 32 529 * 530 + 531 RCL 27 532 RCL 33 533 * 534 - 535 RCL 35 536 * 537 + 538 RCL 22 539 RCL 36 540 * 541 RCL 32 542 RCL 38 543 * 544 + 545 RCL 33 546 RCL 37 547 * 548 - 549 RCL 25 550 * 551 - 552 RCL 23 553 RCL 36 554 * 555 RCL 28 556 RCL 37 557 * 558 + 559 RCL 27 560 RCL 38 561 * 562 - 563 RCL 30 564 * 565 - 566 RCL 19 567 * 568 ST+ 51 569 RCL 61 570 RCL 21 571 * 572 ST+ 51 573 RCL 72 574 RCL 26 575 * 576 RCL 31 577 RCL 39 578 * 579 RCL 34 580 RCL 36 581 * 582 - 583 STO 75 584 RCL 28 585 * 586 - 587 RCL 31 588 RCL 38 589 * 590 RCL 33 591 RCL 36 592 * 593 - 594 STO 76 595 RCL 29 596 * 597 + 598 RCL 22 599 * 600 ST- 51 601 RCL 73 602 RCL 26 603 * 604 RCL 75 605 RCL 27 606 * 607 - 608 RCL 31 609 RCL 37 610 * 611 RCL 32 612 RCL 36 613 * 614 - 615 STO 77 616 RCL 29 617 * 618 + 619 RCL 23 620 * 621 ST+ 51 622 RCL 74 623 RCL 26 624 * 625 RCL 76 626 RCL 27 627 * 628 - 629 RCL 77 630 RCL 28 631 * 632 + 633 RCL 24 634 * 635 ST- 51 636 LBL 52 637 RCL 78 638 RCL 27 639 + 640 STO 82 641 RCL 21 642 * 643 RCL 81 644 - 645 RCL 24 646 RCL 36 647 * 648 - 649 STO 65 650 RCL 60 651 + 652 RCL 15 653 * 654 CHS 655 STO 52 656 RCL 82 657 RCL 20 658 * 659 RCL 22 660 RCL 25 661 * 662 - 663 RCL 23 664 RCL 30 665 * 666 - 667 RCL 24 668 RCL 35 669 * 670 - 671 RCL 16 672 * 673 ST+ 52 674 RCL 78 675 RCL 21 676 + 677 RCL 25 678 * 679 RCL 20 680 RCL 26 681 * 682 - 683 RCL 69 684 - 685 RCL 17 686 * 687 ST+ 52 688 RCL 79 689 RCL 21 690 + 691 RCL 30 692 * 693 RCL 20 694 RCL 31 695 * 696 - 697 RCL 70 698 - 699 RCL 18 700 * 701 ST+ 52 702 RCL 80 703 RCL 21 704 + 705 RCL 35 706 * 707 RCL 20 708 RCL 36 709 * 710 - 711 RCL 71 712 - 713 RCL 19 714 * 715 ST+ 52 716 RCL 60 717 RCL 21 718 * 719 ST- 52 720 RCL 62 721 RCL 22 722 * 723 ST+ 52 724 RCL 63 725 RCL 23 726 * 727 ST+ 52 728 RCL 64 729 RCL 24 730 * 731 ST+ 52 732 RCL 61 733 ST- 52 734 LBL 53 735 RCL 82 736 RCL 21 737 + 738 STO 83 739 RCL 15 740 * 741 RCL 65 742 + 743 RCL 60 744 + 745 RCL 16 746 RCL 20 747 * 748 - 749 RCL 17 750 RCL 25 751 * 752 - 753 RCL 18 754 RCL 30 755 * 756 - 757 RCL 19 758 RCL 35 759 * 760 - 761 STO 53 762 LBL 54 763 RCL 83 764 RCL 15 765 + 766 CHS 767 STO 54 768 LBL 50 769 SF 02 770 XROM 01,06 771 CHS 772 STO 00 773 RCL 51 774 STO 01 775 RCL 52 776 STO 02 777 RCL 53 778 STO 03 779 RCL 54 780 STO 04 781 1 782 STO 05 783 LBL 06 784 RCL 59 785 STO 06 786 STO 22 787 2 788 ST+ 06 789 LBL 07 790 DSE 06 791 GTO 08 792 SF 00 793 SF 21 794 XROM 01,12 795 RTN 796 LBL 08 797 RCL 06 798 1 799 - 800 "a" 801 FIX 00 802 CF 29 803 ARCL X 804 >"=" 805 FIX 04 806 SF 29 807 ARCL IND X 808 PROMPT 809 GTO 07 810 LBL 00 811 LBL 01 812 "N\1D2,3,4,5" 813 AVIEW 814 PSE 815 "XEQ MATH:MATRIX" 816 >" 1ST" 817 AVIEW 818 END