hp41programs

Dblinteg

Triple Integrals for the HP-41


Overview
 

The purpose is to evaluate:   òab ò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 = the name of the subroutine which calculates v(x).
         R06 = the name of the subroutine which calculates w(x;y)
         R07 = the name of the subroutine which calculates t(x;y)

Flags:  /

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 takes x and y from the X and Y registers respectively and calculates 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 ).
 
 

  01  LBL "TIN"  
  02  1.6
  03  STO 09
  04  .15
  05  SQRT
  06  STO 10
  07  RCL 02
  08  RCL 01
  09  STO 12
  10  -
  11  RCL 03
  12  STO 21
  13  /
  14  STO 11
  15  2
  16  /
  17  ST+ 12
  18  CLX
  19  STO 08
  20  LBL 01
  21  RCL 12
  22  RCL 11
  23  RCL 10
  24  *
  25  -
  26  XEQ 02
  27  ST+ 08
  28  RCL 12
  29  XEQ 02
  30  RCL 09
  31  *
  32  ST+ 08
  33  RCL 12
  34  RCL 11
  35  ST+ 12
  36  RCL 10
  37  *
  38  +
  39  XEQ 02
  40  ST+ 08
  41  DSE 21
  42  GTO 01
  43  GTO 06
  44  LBL 02
  45  STO 15
  46  XEQ IND 04
  47  STO 14
  48  RCL 15
  49  XEQ IND 05
  50  RCL 14
  51  -
  52  RCL 03
  53  STO 22
  54  /
  55  STO 13
  56  2
  57  /
  58  ST+ 14
  59  CLX
  60  STO 19
  61  LBL 03
  62  RCL 14
  63  RCL 13
  64  RCL 10
  65  *
  66  -
  67  XEQ 04
  68  ST+ 19
  69  RCL 14
  70  XEQ 04
  71  RCL 09
  72  *
  73  ST+ 19
  74  RCL 14
  75  RCL 13
  76  ST+ 14
  77  RCL 10
  78  *
  79  +
  80  XEQ 04
  81  ST+ 19
  82  DSE 22
  83  GTO 03
  84  RCL 19
  85  RCL 13
  86  *
  87  RTN
  88  LBL 04
  89  STO 17
  90  RCL 15
  91  XEQ IND 06
  92  STO 18
  93  RCL 17
  94  RCL 15
  95  XEQ IND 07
  96  RCL 18
  97  -
  98  RCL 03
  99  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 )
 
 

      STACK        INPUTS      OUTPUTS
           X             /       Integral

 
Example:      Evaluate   I = ò12 òxx^2 òx+yxy  xyz / Ö ( x2 + y2 + z2 ) 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

         n = 1      1  STO 03    XEQ "TIN"  >>>>    I1 =  0.765014888
         n = 2      2  STO 03         R/S         >>>>    I2 =  0.770640690    ( in 4mn 29s )
         n = 4      4  STO 03         R/S         >>>>    I4 =  0.770731245
         n = 8      8  STO 03         R/S         >>>>    I8 =  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 built-in  ò function !

       In 7 FIX format, the result   0.7707326903   is obtained in 3h49mn.
       In 6 FIX format, the result   0.7707326923   is obtained in 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 n6:     I = In + k/n6
 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 .