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

Programs for HP-41

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

Postby eigen » Thu Jul 23, 2009 4:08 am

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 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
eigen
...
...
 
Posts: 7
Joined: Sat Jul 18, 2009 6:24 pm

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

Postby eigen » Tue Sep 01, 2009 4:18 pm

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.000000000000000
ans =
    -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


Return to HP-41 Software

Who is online

Users browsing this forum: No registered users and 1 guest

cron