The Museum of HP Calculators

Triple Integrals for the HP-41

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Introduction

The purpose is to evaluate: ò_a^bò_u(x)^v(x) ò_w(x;y)^t(x;y) f (x;y;z) dx dy dz

The 3 point Gauss-Legendre formula is used in the three directions of x-axis, y-axis and z-axis.

Data Registers: Registers R00 thru R23 are used and some of them must be initialized before executing "TIN", namely:

         R00 = the name of the subroutine which calculates f(x;y;z)
         R01 = a
         R02 = b
         R03 = n = the number of subintervals into which the intervals of integration are to be divided.
         R04 = the name of the subroutine which calculates u(x).
         R05 = ------------------------------------------- v(x).
         R06 = ------------------------------------------- w(x;y)
         R07 = ------------------------------------------- t(x;y)

Flags: none.

Subroutines: 5 subroutines must be keyed into program memory:

      -The first one must take x, y and z from the X-register,the Y-register and the Z-register respectively and calculate f(x;y;z).
      -The second takes x from the X-register and calculates u(x).
      -The third takes x from the X-register and calculates v(x).
      -Another one takes x and y from the X and Y registers respectively and calculates w(x;y)
      -The last one -------------------------------------------------------------------- t(x;y)

N.B: -When the program stops, the result is in X-register and in R08.
-line 43 is a synthetic three-byte GTO ( but it can be replaced by a two-byte GTO
because this line is executed only once ).

001  LBL "TIN"
002  1.6
003  STO 09
004  .15
005  SQRT
006  STO 10
007  RCL 02
008  RCL 01
009  STO 12
010  -
011  RCL 03
012  STO 21
013  /
014  STO 11
015  2
016  /
017  ST+ 12
018  CLX
019  STO 08
020  LBL 01
021  RCL 12
022  RCL 11
023  RCL 10
024  *
025  -
026  XEQ 02
027  ST+ 08
028  RCL 12
029  XEQ 02
030  RCL 09
031  *
032  ST+ 08
033  RCL 12
034  RCL 11
035  ST+ 12
036  RCL 10
037  *
038  +
039  XEQ 02
040  ST+ 08
041  DSE 21
042  GTO 01
043  GTO 06
044  LBL 02
045  STO 15
046  XEQ IND 04
047  STO 14
048  RCL 15
049  XEQ IND 05
050  RCL 14
051  -
052  RCL 03
053  STO 22
054  /
055  STO 13
056  2
057  /
058  ST+ 14
059  CLX
060  STO 19
061  LBL 03
062  RCL 14
063  RCL 13
064  RCL 10
065  *
066  -
067  XEQ 04
068  ST+ 19
069  RCL 14
070  XEQ 04
071  RCL 09
072  *
073  ST+ 19
074  RCL 14
075  RCL 13
076  ST+ 14
077  RCL 10
078  *
079  +
080  XEQ 04
081  ST+ 19
082  DSE 22
083  GTO 03
084  RCL 19
085  RCL 13
086  *
087  RTN
088  LBL 04
089  STO 17
090  RCL 15
091  XEQ IND 06
092  STO 18
093  RCL 17
094  RCL 15
095  XEQ IND 07
096  RCL 18
097  -
098  RCL 03
099  STO 23
100  /
101  STO 16
102  2
103  /
104  ST+ 18
105  CLX
106  STO 20
107  LBL 05
108  RCL 18
109  RCL 16
110  RCL 10
111  *
112  -
113  RCL 17
114  RCL 15
115  XEQ IND 00
116  ST+ 20
117  RCL 18
118  RCL 17
119  RCL 15
120  XEQ IND 00
121  RCL 09
122  *
123  ST+ 20
124  RCL 18
125  RCL 16
126  ST+ 18
127  RCL 10
128  *
129  +
130  RCL 17
131  RCL 15
132  XEQ IND 00
133  ST+ 20
134  DSE 23
135  GTO 05
136  RCL 20
137  RCL 16
138  *
139  RTN
140  LBL 06
141  RCL 08
142  RCL 11
143  *
144  RCL 09
145  2
146  +
147  3
148  Y^X
149  /
150  STO 08
151  END

( 226 bytes / SIZE 024 )

Example: Evaluate I = ò₁²ò_x^{x^2}ò_x+y^xy xyz / Ö ( x² + y² + z² ) dx dy dz

The 5 required subroutines are, for instance:
( with global labels of 6 characters maximum )

01 LBL "FF"
02 ENTER^
03 X^2
04 R^
05 ST* Z
06 X^2
07 +
08 R^
09 ST* Z
10 X^2
11 +
12 SQRT
13 /
14 RTN
15 LBL "U"
16 RTN
17 LBL "V"
18 X^2
19 RTN
20 LBL "W"
21 +
22 RTN
23 LBL "T"
24 *
25 RTN

( All the calculations fit into the stack, but if the functions are very complicated,
one may of course use data registers R24 ...etc... )

Then we initialize registers R00 thru R07:

    "FF"   ASTO 00
      1        STO 01
      2        STO 02
     "U"    ASTO 04
     "V"    ASTO 05
    "W"    ASTO 06
     "T"    ASTO 07

With n = 1      1 STO 03    XEQ "TIN" gives    I₁ = 0.765014888
         n = 2      2 STO 03         R/S            "        I₂ = 0.770640690    ( in 4mn 29s )
         n = 4      4 STO 03         R/S            "        I₄ = 0.770731245
         n = 8      8 STO 03         R/S            "        I₈ = 0.770732669    ( in 4h02mn )                   ( the exact value is 0.7707326901 to ten places )

Unfortunately, with triple integrals, when n is multiplied by 2, execution time is roughly multiplied by 8 !

-An HP-48GX ( in 10 FIX format ) gives 0.7707326900   in 26 hours with the ò function !
                     In 7 FIX format, the result   0.7707326903   is obtained in 3h49mn.
                     In 6 FIX format,   "      "      0.7707326923    "       "       " 1h38mn.
Thus, "TIN" is not so good as the ò function of the HP-48,
but the HP-41 still bears comparison ...

-Execution time can be reduce by keying the 5 subroutines inside the "TIN" program itself ( with local labels ).
The program will run about 17% faster.

-Extrapolation to the limit may also be used to improve accuracy:
There are some theoretical reasons to think that the error is roughly inversely proportional to n⁶: I = I_n + k/n⁶
Applying this formula with n = 4 and n = 8 leads to a system that yields: I = 0.7707326916
the error is only 1.5 10^-9 .

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall