Triple Integrals for the HP-41

Overview
 

 1°) Gauss-Legendre 3-point Formula

  a) Program#1
  b) Program#2

 2°) Gauss-Legendre 4-point Formula
 
 

1°)  Gauss-Legendre 3-point Formula
 

     a) Program#1
 

-"TIN" evaluates:  I = §_a^b §_u(x)^v(x)§_w(x;y)^t(x;y)   f (x;y;z) dx dy dz                       ( § = Integral )
 
 

Data Registers:        •  R00 = f name                         ( Registers R00 thru R07 are to be initialized before executing "TIN" )

         •  R01 = a                                   R08 = I               R09 to R23: temp
         •  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 ).
 
 

 
    ( 226 bytes / SIZE 024 )
 
 

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

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

 
( 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"  >>>>    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 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 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 .
 

     b) Program#2
 

-It's in fact preferable to use just 2 subroutines instead of 4 to calculate u(x) & v(x)  and  w(x,y) & t(x,y)
-We can also use a flag to evaluate also double integrals with the same program:
 

Data Registers:           •  R00 = f name                             ( Registers R00 thru R05 are to be initialized before executing "DTI" )

                                      •  R01 = a                                    •  R04 =  uv name       R06 = I      R07 to R23: temp
                                      •  R02 = b                                    •  R05 = tw name
                                      •  R03 = n = Nb of subintervals

Flags: F02                        CF 02 = Triple Integral
                                          SF 02 = Double Integral

Subroutines:    A subroutine that takes x, y and z from the X-register,the Y-register and the Z-register respectively and returns f(x;y;z) in X-register.
                          A subroutine that takes x in X-register and returns v(x) in Y-register & u(x) in X-register
                          A subroutine that takes x in X-register & y in Y-register and returns t(x,y) in Y-register & w(x,y) in X-register.
 

-Line 21 is a three-byte GTO 01
 
 

 
    ( 216 bytes / SIZE 023 )
 
 

 
Example1:      Evaluate   I = §₁² §_x^{x^2} §_x+y^{x^2.y}  [ Ln(1+x.y.z)] / sqrt ( x² + y² + z² ) dx dy dz
 
 

 
    CF 02   ( triple integral )

    S  ASTO 00

    U  ASTO 04
    V  ASTO 05

    1   STO 01
    2   STO 02

-With n = 1    1  STO 03   XEQ "DTI"   >>>>   I = 1.226398672                     ---Execution time = 48s---
-With n = 2    2  STO 03        R/S          >>>>   I = 1.226795831                     ---Execution time = 5m26s---
-With n = 4    4  STO 03        R/S          >>>>   I = 1.226799708
 

Example2:      Evaluate   I = §₁² §_x^{x^2}   [ Ln(1+x.y)] / sqrt ( x² + y² ) dx dy
 
 

 
   SF 02  ( double integral )

    S  ASTO 00

    U  ASTO 04

    1   STO 01
    2   STO 02

-With n = 1    1  STO 03   XEQ "DTI"   >>>>   I = 0.456387227                            ---Execution time = 15s---
-With n = 2    2  STO 03        R/S          >>>>   I = 0.456373589                            ---Execution time = 49s---
-With n = 4    4  STO 03        R/S          >>>>   I = 0.456373361                            ---Execution time = 184s---
 

2°)  Gauss-Legendre 4-point Formula
 

-A better accuracy may of course be obtained by Gauss-Legendre 4-point formula:
 

     §_-1¹ f(x).dx ~ A.f(-a) + B.f(-b) + B.f(b) + A.f(a)

  with     a = sqrt[(3-sqrt(4.8))/7]      A = 0.5+1/sqrt(43.2)
             b = sqrt[(3+sqrt(4.8))/7]     B = 0.5-1/sqrt(43.2)
 
 

Data Registers:           •  R00 = f name                             ( Registers R00 thru R05 are to be initialized before executing "DTI4" )

                                      •  R01 = a                                    •  R04 =  uv name       R06 = I      R07 to R24: temp
                                      •  R02 = b                                    •  R05 = tw name
                                      •  R03 = n = Nb of subintervals

Flags: F02                        CF 02 = Triple Integral
                                          SF 02 = Double Integral

Subroutines:    A subroutine that takes x, y and z from the X-register,the Y-register and the Z-register respectively and returns f(x;y;z) in X-register.
                          A subroutine that takes x in X-register and returns v(x) in Y-register & u(x) in X-register
                          A subroutine that takes x in X-register & y in Y-register and returns t(x,y) in Y-register & w(x,y) in X-register.
 

-Line 40 is a three-byte GTO 01
 
 

 
     ( 302 bytes / SIZE 025 )
 
 

 
Example1:      Evaluate   I = §₁² §_x^{x^2} §_x+y^{x^2.y}  [ Ln(1+x.y.z)] / sqrt ( x² + y² + z² ) dx dy dz
 
 

 
    CF 02   ( triple integral )

    S  ASTO 00

    U  ASTO 04
    V  ASTO 05

    1   STO 01
    2   STO 02

-With n = 1    1  STO 03   XEQ "DTI"   >>>>   I = 1.226803750                    ---Execution time = 1m51s---
-With n = 2    2  STO 03        R/S          >>>>   I = 1.226799889
-With n = 4    4  STO 03        R/S          >>>>   I = 1.226799707
 

Example2:      Evaluate   I = §₁² §_x^{x^2}   [ Ln(1+x.y)] / sqrt ( x² + y² ) dx dy
 
 

 
   SF 02  ( double integral )

    S  ASTO 00

    U  ASTO 04

    1   STO 01
    2   STO 02

-With n = 1    1  STO 03   XEQ "DTI"   >>>>   I = 0.456373416                          ---Execution time = 26s---
-With n = 2    2  STO 03        R/S          >>>>   I = 0.456373357                          ---Execution time = 90s---
-With n = 4    4  STO 03        R/S          >>>>   I = 0.456373358
 
 

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9
[3]   Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes" in Fortran or in C++" ...............