HP Articles Forum
[Return to the Index ]
[ Previous | Next ]
Eighth Order Linear Equation Solver for the hp 33s
Posted by Palmer O. Hanson, Jr. on 6 Oct 2006, 12:24 a.m.
I started with a TI-59 linear equation solution program by Henrik Ohlsson which appeared in the 81-2 issue of the Swedish newsletter Programbiten. Henrik's program expanded the capability of the TI-59 from order 8 to order 16.
I translated Henrik's program for the Durabrand 828 as a means of demonstrating the capability of that little low cost machine. I was able to fit a sixth order linear equation solution into the Durabrand. The translation was essentially a straightforward, brute force modification from AOS to EOS.
I translated the Durabrand program into a sixth order linear equation solution for the hp 33s. The resulting hp 33s program was essentially a brute force translation of Henrik's program from the TI-59 for the Durabrand and from the Durabrand for the hp 33s. The resulting sixth order linear equation solution was documented as Article 676. The hp 33s program wasn't pretty but it was functional. Only the 26 data registers were used. Register arithmetic was used throughout (no use of the stack in RPN and no use of parentheses or equal signs in ALG) so the the program as written would run in either RPN or ALG mode. That wasn't necessary but it was something that I had wanted to demonstrate.
Next, I found a way to use the statistics registers to expand the capability of the hp 33s program to solve a seventh order set of linear equations. Only limited changes were required to make the expansion. The changes were posted in a September 17 message "Expansion to seventh order" in the Forum.
Finally, I expanded the solution to solve an eighth order set of linear equations. To do that Henrik's method required 33 data registers. The combination of the data registers and the statistics registers of the hp 33s yields only 32 registers. I found a modification to Henrik's method which eliminated the need for one register. The sixth and seventh order solutions had used data registers Y and Z as scratch pads in order to implement the register arithmetic type solution. I had to free up those registers to yield the 32 registers needed for my modification to Henrik's method. To do that I changed substantial parts of the program to use the RPN stack. The resulting program follows.
I won't be doing an expansion to ninth order. Henrik's method would require five more data registers to do that.
Program Listing:
N0001 LBL N N0002 CLVARS N0003 CLSum N0004 INPUT N N0005 STO A N0006 0 N0007 STO N N0008 8 N0009 STO C O0001 LBL O O0002 0 O0003 STO G O0004 RCL A O0005 STO B O0006 1 O0007 STO+ B O0008 8 O0009 STO F O0010 RCL D O0011 IP O0012 STO E O0013 STO D O0014 1 O0015 STO+ D O0016 RCL E O0017 x=0? O0018 GTO C A0001 LBL A A0002 XEQ S A0003 +/- B0001 LBL B B0002 ENTER B0003 ENTER B0004 RCL F B0005 XEQ T B0006 R Down B0007 RCL (i) B0008 x B0009 RCL C B0010 XEQ T B0011 R Down B0012 STO+ (i) B0013 R Down B0014 R Down B0015 1 B0016 STO+ C B0017 STO+ F B0018 R Down B0019 DSE B B0020 GTO B B0021 RCL A B0022 STO B B0023 STO- C B0024 1 B0025 STO+ B B0026 STO- C B0027 DSE E B0028 GTO A C0001 LBL C C0002 XEQ S C0003 STO E C0004 RCL C C0005 XEQ T C0006 RCL (i) C0007 STO+ E C0008 1 C0009 STO- B D0001 LBL D D0002 XEQ S D0003 ENTER D0004 1 D0005 STO+ C D0006 R Down D0007 RCL C D0008 XEQ T D0009 R Down D0010 RCL (i) D0011 + D0012 0 D0013 STO (i) D0014 1 D0015 STO- C D0016 R Down D0017 R Down D0018 RCL C D0019 XEQ T D0020 R Down D0021 STO (i) D0022 RCL E D0023 STO/ (i) D0024 1 D0025 STO+ C D0026 DSE B D0027 GTO D D0028 RCL D D0029 IP D0030 STO E D0031 DSE E D0032 GTO I D0033 GTO H I0001 LBL I I0002 7 I0003 STO C I0004 RCL A I0005 STO B E0001 LBL E E0002 1 E0003 STO+ C E0004 STO+ G E0005 RCL C E0006 XEQ T E0007 0 E0008 x<> (i) F0001 LBL F F0002 ENTER F0003 ENTER F0004 1 F0005 STO+ C F0006 R Down F0007 RCL F F0008 XEQ T F0009 R Down F0010 RCL (i) F0011 +/- F0012 x F0013 RCL C F0014 XEQ T F0015 R Down F0016 RCL (i) F0017 + F0018 0 F0019 STO (i) F0020 R Down F0021 RCL G F0022 STO- C F0023 R Down F0024 RCL C F0025 XEQ T F0026 R Down F0027 STO (i) F0028 RCL G F0029 STO+ C F0030 R Down F0031 1 F0032 STO+ F F0033 R Down F0034 R Down F0035 DSE B F0036 GTO F F0037 RCL A F0038 STO- F F0039 STO B F0040 DSE E F0041 GTO E F0042 RCL G F0043 STO- C F0044 1 F0045 STO+ C G0001 LBL G G0002 RCL F G0003 XEQ T G0004 RCL (i) G0005 RCL C G0006 XEQ T G0007 R Down G0008 STO (i) G0009 RCL F G0010 XEQ T G0011 0 G0012 STO (i) G0013 1 G0014 STO+ C G0015 STO+ F G0016 DSE B G0017 GTO G H0001 LBL H H0002 DSE A H0003 GTO O H0004 7 H0005 STO i H0006 RCL D H0007 IP H0008 STO+ i H0009 1000 H0010 STO/ i H0011 8 H0012 STO+ i J0001 LBL J J0002 RCL (i) J0003 STOP J0004 ISG i J0005 GTO J J0006 0 J0007 1/x J0008 STOP S0001 LBL S S0002 0.01 S0003 STO+ D S0004 RCL D S0005 FIX 2 S0006 STOP S0007 FIX 9 S0008 RTN T0001 LBL T T0002 STO i T0003 27 T0004 - T0005 x<0? T0006 RTN T0007 R Down T0008 1 T0009 STO+ i T0019 RTNTo use the program press XEQ N and see the prompt N? in the display. Enter the order of the solution desired (n = 2 to 6) and press R/S. 1.01 will appear in the display indicating that the machine is ready to accept the A(1,1) element of the matrix Enter the A(1,1) element and press R/S. 1.02 will appear in the display indicating that the machine is ready to accept the A(1,2) element of the matrix. Procced as before to enter the A(1,2) element and see the prompt for the A(1.3) element, and so on. When the column element of the prompt is n+1 enter the corresponding element of the vector. When the last element of the vector has been entered the machine will stop with 8.000000000 in the Y register display and the first element of the solution, X(1), in the X register display. Press R/S to see X(2), and so on. When all of elements of the solution vector have been displayed the message DIVIDE BY 0 will be displayed. The elements of the solution can be recalled starting with register H followed by subsequent registers as needed by the order of the system.. .
If you run the problem with the matrix as a 8x8 sub-Hilbert and the vector all ones you will get the following results:
Exact hp 33s Relative Errorfor a mean relative error of 7.316E-02. Hendrik's program running on the TI-59 yields a mean relative error of 3.331E-02. The ML-02 program running on the TI-59 yields a mean relative error of 2.212E-03. The matrix solution in the HP-28 yield a mean relative error of 1.801E-03. .-72 -79.85564194 0.10674 2520 2752.76793665 0.09237 -27720 -29976.5053369 0.08140 138600 148685.022829 0.07276 -360360 -384064.400797 0.06578 504504 534783.269832 0.06002 -360360 -380245.814081 0.05518 102960 108218.029350 0.05107
Some notes on the program:
1. With Henrik's program on the TI-59 the calculations after the entry of an element required a significant amount of time so that the user was always waiting on the machine. For example, in a fifth order solution when the vector element is entered in response to the prompt 3.06 the TI-59 will run for nine seconds in fast mode or fifteen seconds in normal mode before the next prompt (4.01) appears. The much higher calculating speed of the hp 33s yields much shorter wait times.
2. Henrik used an INV Dsz command in his TI-59 program. I had to synthesize the equivalent in the hp 33s which yielded the curious code at steps D0021 through D0033..
3. The T routine provides the methodology for easily adding the statistics registers immediately after the data registers, essenially skipping the i register.
4. The STO/ (i) command at step D0023 will cause a divide by zero error for some systems. Consider the following problem:
2 -1 -3 1 1where the solution 5.0, 8.25, 0.5, 0.75 can be readily obtained with the TI-59 ML-02 program, with the H-41 Math Pac and with others. If you attempt to solve the problem with this hp 33s program you will get a "Divide by 0" error after the A(2,3) element is entered. If you start the program over again and enter the first row and then the third row you will again get the"Divide by 0" error after the A(2,3) element is entered. However, if you start the program one more time and enter the first row, followed by the fourth row, the second row and the third row in that order you will see the correct solution. The point of all this is that if you get the "Divide by 0" error after the entry of an element you must try another arrangement of the rows.-2 1 1 3 1
2 -1 0 -1 1
-3 2 -1 0 1