The Museum of HP Calculators

Numerical Integration for the HP-41

This program is Copyright © 2005 by Jean-Marc Baillard and is used here by permission.

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.

Overview

 1°) Continuous Data

  a)  Gauss-Legendre 3-point Formula
  b)  Gauss-Legendre 16-point ( & 48-point ) Formula
  c)  Gauss-Chebyshev Formula
  d)  Filon's Formulae
  e)  An Example of Curvilinear Integral

 2°) Discrete Data

  a)  Equally Spaced Abscissas

   a-1) Simpson's Rule
   a-2) Other Newton-Cotes Formulae

  b)  Unequally Spaced Abscissas

   b-1) Trapezoidal Rule  ( new )
   b-2) Connecting Parabolic Segments
   b-3) Connecting Cubic Segments ( new )
   b-4) Natural Cubic Spline ( new )
   b-5) Integration of the Lagrange Polynomial ( new )

  c)  Double Integrals ( Equally Spaced Arguments only )
  d)  Triple Integrals ( Equally Spaced Arguments only )
 

-For continuous data, see also  "Double Integrals for the HP-41"
                                                "Triple Integrals for the HP-41"
                                         and "Multiple Integration for the HP-41"

-Integrals are denoted §      §_a^b f(x).dx
 

1°) Continuous Data
 

        a) Gauss-Legendre 3-Point  Formula
 

-Gauss-Legendre formulas are defined by  §_-1¹ f(x).dx  ~  w₁.x₁ + w₂.x₂ + ....... + w_n.x_n
  where w_i and x_i  are choosen so that the formula produces perfect accuracy when  f(x) is a polynomial of degree < 2n
-The  x_i  are the zeros of the Legendre polynomials and they can't be easily computed for large n-values.

-The 3-point Gauss-Legendre formula is   §_-1¹ f(x).dx  ~  (5/9).f(-(3/5)^1/2) + (8/9).f(0) + (5/9).f((3/5)^1/2)
-Then, a linear change of variable transforms  [-1;1]  into any ( finite ) interval [a;b]

-Actually, the interval [a;b] is divided into  n sub-intervals  [ a+k.(b-a)/n ; a+(k+1).(b-a)/n ]    k = 0 , 1 , ...... , n-1
 and, at least theoretically, the results tend to the exact integral as n tends to infinity
 ( provided f and its derivative have no "strong" singularity ).
 

Data Registers:       •  R00 = Function name                             ( Registers R00 thru R03 are to be initialized before executing "GL3" )
                                  •  R01 = a
                                  •  R02 = b
                                  •  R03 = n = number of subintervals over which the formula is applied.          R04 =  §_a^b f(x).dx        R05 thru R08: temp
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.
 

01  LBL "GL3"
02  RCL 02
03  RCL 01
04  STO 07
05  -
06  RCL 03
07  STO 05
08  /
09  STO 06
10  2
11  /
12  ST+ 07
13  .6
14  SQRT
15  *
16  STO 08
17  CLX
18  STO 04
19  LBL 01
20  RCL 07
21  RCL 08
22  -
23  XEQ IND 00
24  ST+ 04
25  RCL 07
26  XEQ IND 00
27  1.6
28  *
29  ST+ 04
30  RCL 07
31  RCL 08
32  +
33  XEQ IND 00
34  ST+ 04
35  RCL 06
36  ST+ 07
37  DSE 05
38  GTO 01
39  RCL 04
40  *
41  3.6
42  /
43  STO 04
44  END

( 67 bytes / SIZE 009 )
 

Example:     Evaluate  §₁³  e^{-x^2} dx

-Load this short routine:

  01  LBL "T"   ( any global Label , maximum 6 characters )
  02  X^2
  03  CHS
  04  E^X
  05  RTN

 "T"  ASTO 00
  1    STO 01
  3    STO 02       if  n = 2    2  STO 03  XEQ "GL3"  >>>>  0.139390854
                               n = 4    4  STO 03       R/S          >>>>  0.139383255
                               n = 8    8  STO 03       R/S          >>>>  0.139383216      ( execution time = 17 seconds )

-If a or b is infinite, make a change of variable ( like x = Tan u ) to reduce [ a ; b ] to a finite interval.
 

        b) Gauss-Legendre 16-Point Formula
 

-The 16 coefficients w_i and x_i  are to be stored in registers R11 to R26 as listed below.
 

Data Registers:       •  R00 = Function name            ( Registers R00 thru R03 and R11 thru R26  are to be initialized before executing "GL16" )
                                  •  R01 = a
                                  •  R02 = b
                                  •  R03 = n = number of subintervals over which the formula is applied.    R04 = Integral         R05 thru R10: temp

                                  •  R11 = 0.02715245941         •  R19 = 0.1495959888
                                  •  R12 = 0.9894009350           •  R20 = 0.6178762444
                                  •  R13 = 0.06225352394         •  R21 = 0.1691565194
                                  •  R14 = 0.9445750231           •  R22 = 0.4580167777
                                  •  R15 = 0.09515851168         •  R23 = 0.1826034150
                                  •  R16 = 0.8656312024           •  R24 = 0.2816035508
                                  •  R17 = 0.1246289713           •  R25 = 0.1894506105
                                  •  R18 = 0.7554044084           •  R26 = 0.09501250984
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.
 

01  LBL "GL16"
02  CLX
03  STO 04
04  RCL 02
05  RCL 01
06  STO 07
07  -
08  RCL 03
09  STO 05
10  ST+ X
11  /
12  STO 06
13  LBL 01
14  ST+ 07
15  26.01
16  STO 08
17  LBL 02
18  RCL 07
19  RCL IND 08
20  RCL 06
21  *
22  STO 09
23  -
24  XEQ IND 00
25  STO 10
26  RCL 07
27  RCL 09
28  +
29  XEQ IND 00
30  RCL 10
31  +
32  DSE 08
33  RCL IND 08
34  *
35  ST+ 04
36  DSE 08
37  GTO 02
38  RCL 06
39  ST+ 07
40  DSE 05
41  GTO 01
42  ST* 04
43  RCL 04
44  END

( 71 bytes / SIZE 027 )
 

Example:     Evaluate  §₀³  e^{-x^4} dx

  01  LBL "T"   ( any global Label , maximum 6 characters )
  02  X^2
  03  X^2
  04  CHS
  05  E^X
  06  RTN

 "T"  ASTO 00
  0    STO 01
  3    STO 02       if  n = 1    1  STO 03  XEQ "GL16"  >>>>  0.906402825
                               n = 2    2  STO 03        R/S           >>>>  0.906402476    ( in 28 seconds )

Remark:

-If you want to use the 48-point formula,
    replace line 01 by  LBL "GL48"
    ----------- 15 by  58.01            ( instead of 26.01 )
  and store these 48 coefficients in registers R11 to R58
 

    R11 = 0.003153346052           R23 = 0.02742650971             R35 = 0.04761665849          R47 = 0.06070443917
    R12 = 0.9987710073               R24 = 0.9058791367               R36 = 0.6778723796            R48 = 0.3487558863
    R13 = 0.007327553901           R25 = 0.03116722783             R37 = 0.05035903555          R49 = 0.06203942316
    R14 = 0.9935301723               R26 = 0.8765720203               R38 = 0.6288673968            R50 = 0.2873624874
    R15 = 0.01147723458             R27 = 0.03477722256             R39 = 0.05289018949          R51 = 0.06311419229
    R16 = 0.9841245837               R28 = 0.8435882616               R40 = 0.5772247261            R52 = 0.2247637904
    R17 = 0.01557931572             R29 = 0.03824135107             R41 = 0.05519950370          R53 = 0.06392423858
    R18 = 0.9705915925               R30 = 0.8070662040               R42 = 0.5231609747            R54 = 0.1612223561
    R19 = 0.01961616046             R31 = 0.04154508294             R43 = 0.05727729210          R55 = 0.06446616444
    R20 = 0.9529877032               R32 = 0.7671590325               R44 = 0.4669029048            R56 = 0.09700469921
    R21 = 0.02357076084             R33 = 0.04467456086             R45 = 0.05911483970          R57 = 0.06473769681
    R22 = 0.9313866907               R34 = 0.7240341309               R46 = 0.4086864820            R58 = 0.03238017096
 

        c) Gauss-Chebyshev Formula
 

-This program calculates the singular integral:    §_a^b  ((x-a)(x-b))^-1/2  f(x).dx  ~  (pi/n) [ f(x₁) + f(x₂) + ..... + f(x_n) ]
  where   x_i = (a+b)/2 + ((b-a)/2).cos((2i-1)pi/(2n))
-Unlike "GL3" and "GL16", "GC" determines the coefficients w_i and x_i directly since they are much easier to obtain.
 

Data Registers:       •  R00 = Function name        ( Registers R00 thru R03 are to be initialized before executing "GC" )
                                  •  R01 = a
                                  •  R02 = b
                                  •  R03 = n          R04 = Integral         R05 thru R07: temp
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.
 

01  LBL "GC"
02  RCL 03
03  STO 05
04  RCL 02
05  RCL 01
06  -
07  2
08  /
09  STO 06
10  CLX
11  STO 04
12  LBL 01
13  1                       CLX  SIGN  is faster
14  ASIN
15  RCL 05
16  ST+ X
17  DSE X              or  LASTX   -
18  RCL 03
19  /
20  *
21  COS
22  RCL 06
23  ST* Y
24  +
25  RCL 01
26  +
27  XEQ IND 00
28  ST+ 04
29  DSE 05
30  GTO 01
31  PI
32  RCL 03
33  /
34  ST* 04
35  RCL 04
36  END

( 51 bytes / SIZE 007 )
 

Example:     Compute  §₁³  e^x(-x²+4.x-3)^-1/2 dx  =  §₁³  e^x((3-x)(x-1))^-1/2 dx

  01  LBL "T"   ( any global Label , maximum 6 characters )
  02  E^X
  03  RTN

 "T"  ASTO 00
  1    STO 01
  3    STO 02       if  n = 2    2  STO 03  XEQ "GC"  >>>>  29.262
                               n = 4    4  STO 03       R/S         >>>>  29.389695
                               n = 8    8  STO 03       R/S         >>>>  29.38969917
                              n = 16  16  STO 03      R/S         >>>>  29.38969918    ( in 25 seconds )

Note:   This program works in all angular mode,
             but "GC" and the subroutine must be executed in the same angular mode.
 

        d) Filon's Formulae
 

-This method computes integrals of the form:    §_a^b  f(x).cos(k.x).dx   and   §_a^b  f(x).sin(k.x).dx

    §_a^b  f(x).cos(k.x).dx  ~  h.[ A(k.h) (  f(x_2n).sin(k.x_2n)  -  f(x₀).sin(k.x₀) ) + B(k.h) C_2n + C(k.h) C_2n-1 ]                 a = x₀ ;  b = x_2n ;  h = (b-a)/(2n)
    §_a^b  f(x).sin(k.x).dx   ~  h.[ A(k.h) ( -f(x_2n).cos(k.x_2n) + f(x₀).cos(k.x₀) ) + B(k.h) S_2n + C(k.h) S_2n-1 ]

where   C_2n = f(x₀).cos(k.x₀) + f(x₂).cos(k.x₂) + ......... + f(x_2n).cos(k.x_2n)  -  (1/2) ( f(x_2n).cos(k.x_2n) + f(x₀).cos(k.x₀) )
            S_2n = f(x₀).sin(k.x₀) + f(x₂).sin(k.x₂) + ......... + f(x_2n).sin(k.x_2n)  -  (1/2) ( f(x_2n).sin(k.x_2n) + f(x₀).sin(k.x₀) )

            C_2n-1 = f(x₁).cos(k.x₁) + f(x₃).cos(k.x₃) + ......... + f(x_2n-1).cos(k.x_2n-1)
            S_2n-1 = f(x₁).sin(k.x₁) + f(x₃).sin(k.x₃) + ......... + f(x_2n-1).sin(k.x_2n-1)

 and     A(µ) = 1/µ + (sin(2µ))/(2µ²) - 2(sin²µ)/µ³
           B(µ) =  2(1+cos²µ)/µ² - 2(sin(2µ))/µ³
           C(µ) =  4(sinµ)/µ³ - 4(cosµ)/µ²

-If  k = 0 , Simpson's rule is used
 

Data Registers:       •  R00 = Function name        ( Registers R00 thru R03 are to be initialized before executing "FIL" )
                                  •  R01 = a
                                  •  R02 = b         R06 = k             R07 to R12: temp
                                  •  R03 = n
                  When the program stops:   R04 =   §_a^b  f(x).cos(k.x).dx  = X-register
                                                           R05 =   §_a^b  f(x).sin(k.x).dx   = Y-register
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.
 

  01  LBL "FIL"
  02  RAD
  03  STO 06
  04  RCL 02
  05  RCL 01
  06  -
  07  RCL 03
  08  ST+ X
  09  STO 07
  10  /
  11  STO 08
  12  *
  13  STO 09
  14  CLX
  15  STO 04
  16  STO 05
  17  STO 10
  18  1.5
  19  1/X
  20  STO 11
  21  ST+ X
  22  STO 12
  23  RCL 09
  24  X=0?
  25  GTO 01
  26  ST+ X
  27  SIN
  28  STO 12
  29  LASTX
  30  /
  31  RCL 09
  32  SIN
  33  LASTX
  34  /
  35  STO 11
  36  X^2
  37  ST+ X
  38  -
  39  1
  40  +
  41  STO 10
  42  RCL 09
  43  COS
  44  X^2
  45  1
  46  +
  47  RCL 12
  48  RCL 09
  49  /
  50  -
  51  ST+ X
  52  X<> 11
  53  RCL 09
  54  COS
  55  -
  56  4
  57  *
  58  RCL 09
  59  ST/ 10
  60  X^2
  61  ST/ 11
  62  /
  63  STO 12
  64  GTO 01
  65  LBL 00
  66  RCL 06
  67  *
  68  X<>Y
  69  P-R
  70  RCL X
  71  RCL 10
  72  *
  73  ST+ 05
  74  CLX
  75  RCL 11
  76  *
  77  ST+ 04
  78  CLX
  79  RCL 10
  80  *
  81  ST- 04
  82  CLX
  83  RCL 11
  84  *
  85  ST+ 05
  86  RTN
  87  LBL 01
  88  RCL 01
  89  RCL 07
  90  RCL 08
  91  *
  92  +
  93  STO 09
  94  XEQ IND 00
  95  RCL 06
  96  RCL 09
  97  *
  98  X<>Y
  99  P-R
100  RCL 12
101  X<> 11
102  STO 12
103  *
104  ST+ 04
105  CLX
106  RCL 12
107  *
108  ST+ 05
109  DSE 07
110  GTO 01
111  2
112  ST/ 11
113  RCL 01
114  XEQ IND 00
115  RCL 01
116  XEQ 00
117  RCL 02
118  XEQ IND 00
119  CHS
120  RCL 02
121  XEQ 00
122  RCL 08
123  ST* 04
124  ST* 05
125  RCL 05
126  RCL 04
127  END

( 167 bytes / SIZE 013 )
 
 

STACK	INPUT	OUTPUTS
Y	/	§ab  f(x).sin(k.x).dx
X	k	§ab  f(x).cos(k.x).dx

Example:   Evaluate   I = §₁⁶  Ln(x).cos(10x).dx    and   J = §₁⁶  Ln(x).sin(10x).dx

 01  LBL "T"
 02  LN
 03  RTN

 "T"  ASTO 00
  1    STO 01
  6    STO 02

with n = 8   ;    8   STO 03        10      R/S     >>>>    I =  -0.047890755
                                                                     X<>Y   J =   0.175512930

with n = 16  ;  16  STO 03        10      R/S     >>>>    I =  -0.047429223
                                                                     X<>Y   J =   0.174731804

with n = 32  ;  32  STO 03        10      R/S     >>>>    I =  -0.047453034
                                                                     X<>Y   J =   0.174714501

with n = 64  ;  64  STO 03        10      R/S     >>>>    I =  -0.047454443     ( exact results are  I = -0.047454534
                                                                     X<>Y   J =   0.174713854                           and  J =  0.174713817   to 9 decimals )
 

-The complete Filon's formulas involve the 3rd derivative of the function f
  but this term is omitted here ( f'''(x) is difficult to evaluate on an HP-41 ... )
 

        e) An Example of Curvilinear Integral
 

-The following program computes    §_(C)  f(x,y).ds    where  (C) is the cicumference of the circle:  x² + y² = R²     (  ds = ( dx²+dy²)^1/2 )
-The function f is evaluated 2n times.
 

Data Registers:        •  R00:  function name       ( Registers R00 thru R02 are to be initialized before executing "CINT" )
                                   •  R01 = R
                                   •  R02 = n
                         when the program stops:  R03 = §_(C)  f(x,y).dx         R04-R05: temp

Flags: /
Subroutine:  A program which calculates  f(x,y)  assuming x is in X-register and y is in Y-register upon entry.
 

01  LBL "CINT"
02  RCL 02
03  ST+ X
04  STO 04
05  CLX
06  STO 03
07  SIGN
08  CHS
09  ACOS
10  RCL 02
11  /
12  STO 05
13  LBL 01
14  RCL 04
15  RCL 05
16  *
17  RCL 01
18  P-R
19  XEQ IND 00
20  ST+ 03
21  DSE 04
22  GTO 01
23  RCL 03
24  PI
25  *
26  RCL 01
27  *
28  RCL 02
29  /
30  STO 03
31  END

( 45 bytes / SIZE 006 )
 

Example:   Evaluate    §_(C)  Ln(3+x.y).ds    where  (C):  x² + y² = 1

 01  LBL "T"
 02  *
 03  3
 04  +
 05  LN
 06  RTN

 "T"  ASTO 00    R = 1  whence  1  STO 01

-with  n =  4    4   STO 02    XEQ "CINT"  >>>>  6.858533883
-with  n =  8    8   STO 02          R/S           >>>>  6.858689700
-with  n = 16  16  STO 02          R/S           >>>>  6.858689706
 
 

2°) Discrete Data
 

        a) Equally Spaced Abscissas
 

        a-1) Simpson's Rule
 

-We want to evaluate  §_x1^xn  f(x).dx  but we only know   f(x₁) , f(x₂) , ....... , f(x_n)  where the  x_i  are equally spaced     x_i+1 - x_i = h    for i = 1,2,...,n-1

-If n is odd, we use the formula:    §_x1^x3  f(x).dx  ~  h[f(x₁)+4f(x₂)+f(x₃)]/3    over the intervals  [x₁;x₃] ,  [x₃;x₅] , ..... ,  [x_n-2;x_n]
-If n is even,  we use  §_x1^x4  f(x).dx  ~  3h[f(x₁)+3f(x₂)+3f(x₃)+f(x₄)]/8  and the same formula as above over  [x₄;x₆] , ..... ,  [x_n-2;x_n]

-These formulas produce perfect accuracy if f is a polynomial of degree < 4

-We assume  n > 2
 

Data Registers:       •  R00 = h                                                              ( Registers R00 thru Rnn  are to be initialized before executing "SIMP" )
                                  •  R01 = f(x₁)  •  R02 = f(x₂)  .......  •  Rnn = f(x_n)
Flags:  /
Subroutines:  /
 

-Synthetic register M may of course be replaced by any unused data register, for instance R99 ( provided n < 99 )
 

01  LBL "SIMP"
02  STO M
03  4
04  SQRT
05  RCL IND Y
06  CHS
07  LBL 01
08  RCL IND M
09  LASTX
10  X<> T
11  *
12  ST+ Y
13  RDN
14  DSE M
15  GTO 01
16  X<>Y
17  4
18  -
19  X=0?
20  GTO 02
21  CLX                  Lines 21 to 36 are only useful when n is even
22  RCL 02
23  11
24  *
25  RCL 01
26  15
27  *
28  -
29  RCL 03
30  5
31  *
32  -
33  RCL 04
34  +
35  8
36  /
37  LBL 02
38  +
39  RCL 01
40  -
41  RCL 00
42  *
43  3
44  /
45  END

( 64 bytes / SIZE n+1 )
 
 

STACK	INPUT	OUTPUT
X	n	§x1xn  f(x).dx

 where n = the number of points ,  n > 2

Example:   Evaluate  §₀^pi/2  f(x).dx  and  §₀^5pi/12  f(x).dx  from the following values
 

x	0	pi/12	pi/6	pi/4	pi/3	5pi/12	pi/2
f(x)	0	0.2588190	0.5	0.7071068	0.8660254	0.9659258	1

  0  STO 01   0.2588190  STO 02  ................  1  STO 07

  PI  12  /   STO 00

   7   XEQ "SIMP"  >>>>  1.0000263   ( exact result = §₀^pi/2  sin(x).dx = 1 )
   6         R/S          >>>>   0.7412102   ( exact result = §₀^5pi/12  sin(x).dx = 0.7411810 )
 

        a-2) Other Newton-Cotes Formulae
 

-We assume  n = 6k+1  ( where k is an integer )   The  x_i  are still  equally spacedx_i+1 - x_i = h    for i = 1,2,...,n-1

Formula:   §_x1^x7  f(x).dx  ~  h[41f(x₁)+216f(x₂)+27f(x₃)+272f(x₄)+27f(x₅)+216f(x₆)+41f(x₇)]/140

Data Registers:       •  R00 = h                                                                ( Registers R00 thru Rnn  are to be initialized before executing "NC6" )
                                  •  R01 = f(x₁)  •  R02 = f(x₂)  .......  •  Rnn = f(x_n)
Flags:  /
Subroutines:  /
 

01  LBL "NC6"
02  RCL IND X
03  41
04  *
05  DSE Y
06  LBL 01
07  RCL IND Y
08  216
09  *
10  +
11  DSE Y
12  RCL IND Y
13  27
14  *
15  +
16  DSE Y
17  RCL IND Y
18  272
19  *
20  +
21  DSE Y
22  RCL IND Y
23  27
24  *
25  +
26  DSE Y
27  RCL IND Y
28  216
29  *
30  +
31  DSE Y
32  RCL IND Y
33  82
34  *
35  +
36  DSE Y
37  GTO 01
38  RCL 01
39  41
40  *
41  -
42  RCL 00
43  *
44  140
45  /
46  END

( 82 bytes / SIZE nnn+1 )
 
 

STACK	INPUT	OUTPUT
X	n	§x1xn  f(x).dx

 where n = the number of points ;  n = 7 , 13 , 19 , .....

Example:   Evaluate  §₀^pi/2  f(x).dx   from the following values
 

x	0	pi/12	pi/6	pi/4	pi/3	5pi/12	pi/2
f(x)	0	0.2588190	0.5	0.7071068	0.8660254	0.9659258	1

  0  STO 01   0.2588190  STO 02  ................  1  STO 07

  PI  12  /   STO 00
   7   XEQ "NC6"  >>>>  1.0000000   ( the result is correct to 7 decimals! )
 

-Newton-Cotes 10 point formula is also interesting for sets of (9k+1) equally spaced arguments ( if you want to write a "NC9" program ):

     §_x1^x10  f(x).dx  ~  9h[2857(f(x₁)+f(x₁₀))+15741(f(x₂)+f(x₉))+1080(f(x₃)+f(x₈))+19344(f(x₄)+f(x₇))+5778(f(x₅)+f(x₆))]/89600

-Higher degree formulas are theoretically better, but they involve large coefficients of alternate signs which may lead to poor accuracy...
 

        b) Unequally Spaced Abscissas
 

        b-1) Trapezoidal Rule
 

-We assume a function f is only known by n data points
 

x1	x2	...........	xn
f(x1)	f(x2)	...........	f(xn)

-The trapezoidal rule is the easiest formula to evaluate  §_x1^xn  f(x).dx
-We simply add the areas of n trapezoids:  Sum_{i=1,2,...,n-1}   (x_i+1-x_i).[ f(x_i+1)+f(x_i) ]/2
 

Data Registers:       •  R00 = n = number of points ,                              ( Registers R00 thru R2n  are to be initialized before executing "TRAP" )

                                  •  R01 = x₁     •  R03 = x₂      .......  •  R2n-1 = x_n
                                  •  R02 = f(x₁)  •  R04 = f(x₂)  .......  •  R2n   = f(x_n)
Flags:  /
Subroutines:  /
 

01  LBL "TRAP"
02  RCL 00
03  DSE X
04  ST+ X
05  STO M
06  CLX
07  LBL 01
08  RCL M
09  2
10  +
11  X<>Y
12  RCL IND Y
13  RCL IND M
14  +
15  DSE Z
16  DSE M
17  RCL IND Z
18  RCL IND M
19  -
20  *
21  +
22  DSE M
23  GTO 01
24  2
25  /
26  END

( 47 bytes / SIZE 2n+1 )
 
 

STACK	INPUT	OUTPUT
X	/	§x1xn  f(x).dx

 where n = the number of points.

Example:     Calculate  §₁⁸  f(x).dx   from the following values
 

x	1	2.4	4	5.2	7	8
f(x)	1	4	6	5	4	2

Store these 12 numbers into registers   R01  R03  R05  R07  R09  R11   ( x )
                                                          R02  R04  R06  R08  R10  R12   ( f(x) )    respectively

-There are 6 points so,     6   STO 00    XEQ "TRAP"    >>>>    §₁⁸  f(x).dx  ~  29.2
 

        b-2) Connecting Parabolic Segments
 

-More accurate results may be obtained if we use several connected parabolic segments to compute  §_x1^xn  f(x).dx
-However, if n is even,  §_x1^x2  f(x).dx  is calculated by a polynomial of degree 3
-Thus, "ITG" yields the same results as "SIMP"  when the  x_i's are equally spaced.
 

Data Registers:       •  R00 = n = number of points ,        n > 2          ( Registers R00 thru R2n  are to be initialized before executing "ITG" )

                                  •  R01 = x₁     •  R03 = x₂      .......  •  R2n-1 = x_n
                                  •  R02 = f(x₁)  •  R04 = f(x₂)  .......  •  R2n   = f(x_n)
Flags:  /
Subroutines:  /

-Synthetic registers M N O P may be replaced by any unused data registers, for instance R96  R97  R98  R99 if  n < 48
-In this case, replace line 02 ( CLA ) by  CLX  STO 98
 

  01  LBL "ITG"
  02  CLA
  03  RCL 00
  04  2
  05  -
  06  ST+ X
  07  STO M
  08  4
  09  MOD
  10  X#0?                        Lines 12 to 86 are only useful when n is even. A slightly less accurate result may be obtained by replacing lines 12 to 86 by:
  11  GTO 01
  12  RCL 01                     RCL 06
  13  RCL 03                     RCL 05
  14  +                                RCL 01
  15  2                                -
  16  /                                 RCL 03
  17  STO O                      RCL 05
  18  RCL 05                     -
  19  -                                *
  20  STO N                      /
  21  RCL 07                     RCL 03
  22  RCL 01                     RCL 01
  23  -                                -
  24  /                                STO O
  25  RCL 07                     *
  26  RCL 03                     RCL 04
  27  -                                RCL 05
  28  /                                RCL 03
  29  RCL 08                     -
  30  *                                /
  31  RCL O                      +
  32  RCL 07                     RCL 02
  33  -                                RCL 05
  34  ST* N                       RCL 01
  35  RCL 05                     -
  36  RCL 01                     /
  37  -                                -
  38  /                                RCL O
  39  RCL 05                     *
  40  RCL 03                     RCL 02
  41  -                                RCL 04
  42  /                                 +
  43  RCL 06                      3
  44  *                                *
  45  -                                 +
  46  RCL 05                      ST* O
  47  RCL 07
  48  -
  49  /
  50  RCL 03
  51  RCL 01
  52  -
  53  X^2
  54  *
  55  STO O
  56  RCL 02
  57  RCL 01
  58  RCL 05
  59  -
  60  /
  61  RCL 01
  62  RCL 07
  63  -
  64  /
  65  RCL 04
  66  RCL 03
  67  RCL 05
  68  -
  69  /
  70  RCL 03
  71  RCL 07
  72  -
  73  /
  74  +
  75  RCL N
  76  *
  77  ST+ X
  78  RCL 02
  79  +
  80  RCL 04
  81  +
  82  ST+ O
  83  RCL 03
  84  RCL 01
  85  -
  86  ST* O
  87  LBL 01
  88  DSE M
  89  4
  90  RCL M
  91  ST+ Y
  92  2
  93  +
  94  RCL IND M
  95  RCL IND Y
  96  -
  97  STO P                   ( synthetic )
  98  RCL IND Z
  99  LASTX
100  -
101  STO N
102  /
103  PI
104  INT
105  ST+ Z
106  DSE X
107  +
108  RCL IND Y
109  *
110  RCL P
111  RCL N
112  -
113  ST* X
114  LASTX
115  /
116  RCL P
117  /
118  DSE Z
119  DSE Z
120  RCL IND Z
121  *
122  -
123  RCL N
124  RCL P
125  ST- N
126  /
127  2
128  ST- T
129  ST- M
130  +
131  RCL IND Z
132  *
133  +
134  RCL N
135  *
136  ST+ O
137  DSE M
138  GTO 01
139  X<> O
140  6
141  /
142  CLA
143  END

( 192 bytes / SIZE 2n+1 )
 
 

STACK	INPUT	OUTPUT
X	/	§x1xn  f(x).dx

 where n = the number of points.  ( n > 2 )

Example:     Calculate  §₁⁷  f(x).dx  and  §₁⁸  f(x).dx   from the following values
 

x	1	2.4	4	5.2	7	8
f(x)	1	4	6	5	4	2

Store these 12 numbers into registers   R01  R03  R05  R07  R09  R11   ( x )
                                                          R02  R04  R06  R08  R10  R12   ( f(x) )    respectively
 

-Five points:   5   STO 00   XEQ "ITG"  >>>>  §₁⁷  f(x).dx  ~  26.4226
-Six points:     6   STO 00        R/S         >>>>  §₁⁸  f(x).dx  ~  30.5339

-With the modification listed on the right, it gives  §₁⁸  f(x).dx  ~  30.6457

Remark:   We can also use Lagrange interpolation to obtain equally spaced arguments ( cf for instance "Lagrange Interpolation Formula for the HP-41" )

-With this set of 6 data points, it yields:
 

x	1	2	3	4	5	6	7	8
f(x)	1	2.8570	5.3453	6	5.2069	4.3568	4	2

 and "SIMP" now produces   §₁⁸  f(x).dx  ~  29.6997

-In this example, we can even use the Newton-Cotes 8 point formula

    §_x1^x8  f(x).dx  ~  7h[751f(x₁)+3577f(x₂)+1323f(x₃)+2989f(x₄)+2989f(x₅)+1323f(x₆)+3577f(x₇)+751f(x₈)]/17280
    §_x1^x8  f(x).dx  ~  29.6179

 which is probably the best evaluation of this integral without further information about the function.

-For larger sets of data points, one may likewise use Lagrange interpolation to transform unequally spaced abscissas
  into (6k+1) or (9k+1) equally spaced abcsissas, and then apply "NC6" or "NC9" ( cf 2°) a-2) above ).
-Take for instance successive groups of  7 or 10 points with "LAGR"
  but too many points could lead to worse rather than better accuracy!
 

        b-3) Connecting Cubic Segments  ( X-Functions module required )
 

-The following program uses several segments of cubic curves to evaluate  §_x1^xn  f(x).dx
  so that the method produces a perfect accuracy if the function is a polynomial of degree < 4

-Registers R00 to R2n are temporarily modified during the calculations, but their contents are finally restored.
 

Data Registers:       •  R00 = n = number of points ,        n > 3                   ( Registers R00 thru R2n  are to be initialized before executing "ITG3" )

                                  •  R01 = x₁     •  R03 = x₂      .......  •  R2n-1 = x_n
                                  •  R02 = f(x₁)  •  R04 = f(x₂)  .......  •  R2n   = f(x_n)
Flags:  /
Subroutines:  /
 

  01  LBL "ITG3"
  02  CLA
  03  RCL 00
  04  RCL 00
  05   E3
  06  ST/ 00
  07  SIGN
  08  -
  09  3
  10  MOD
  11  -
  12  ST+ 00
  13  LBL 01
  14  ISG 00
  15  FS? 30
  16  GTO 02
  17  RCL 01
  18  RCL 03
  19  +
  20  2
  21  /
  22  STO O
  23  RCL 05
  24  -
  25  STO N
  26  RCL 07
  27  RCL 01
  28  -
  29  /
  30  RCL 07
  31  RCL 03
  32  -
  33  /
  34  RCL 08
  35  *
  36  RCL O
  37  RCL 07
  38  -
  39  ST* N
  40  RCL 05
  41  RCL 01
  42  -
  43  /
  44  RCL 05
  45  RCL 03
  46  -
  47  /
  48  RCL 06
  49  *
  50  -
  51  RCL 05
  52  RCL 07
  53  -
  54  /
  55  RCL 03
  56  RCL 01
  57  -
  58  X^2
  59  *
  60  STO O
  61  RCL 02
  62  RCL 01
  63  RCL 05
  64  -
  65  /
  66  RCL 01
  67  RCL 07
  68  -
  69  /
  70  RCL 04
  71  RCL 03
  72  RCL 05
  73  -
  74  /
  75  RCL 03
  76  RCL 07
  77  -
  78  /
  79  +
  80  RCL N
  81  *
  82  ST+ X
  83  RCL 02
  84  +
  85  RCL 04
  86  +
  87  RCL O
  88  +
  89  RCL 03
  90  RCL 01
  91  -
  92  *
  93  ST+ M
  94  XEQ 00
  95  GTO 01
  96  LBL 00
  97  RCL 00
  98  FRC
  99   E-3
100  ST- Y
101  *
102  ST+ X
103  3
104  LASTX
105  +
106  +
107  REGSWAP
108  RTN
109  LBL 02
110  RCL 00
111  2
112  -
113  3
114  /
115  INT
116  ST- 00
117  LBL 03
118  RCL 04
119  RCL 01
120  RCL 03
121  -
122  STO N
123  LASTX
124  RCL 07
125  -
126  *
127  /
128  RCL 06
129  RCL 01
130  RCL 05
131  -
132  STO O
133  /
134  RCL 05
135  RCL 07
136  -
137  /
138  -
139  RCL 01
140  RCL 07
141  -
142  *
143  RCL 03
144  RCL 05
145  -
146  /
147  RCL 02
148  RCL N
149  /
150  RCL O
151  /
152  -
153  RCL 08
154  RCL 03
155  RCL 07
156  -
157  /
158  RCL 05
159  RCL 07
160  -
161  /
162  -
163  2
164  /
165  RCL O
166  RCL 04
167  *
168  RCL N
169  /
170  RCL 03
171  RCL 05
172  -
173  /
174  RCL 03
175  RCL 07
176  -
177  /
178  -
179  RCL N
180  RCL 06
181  *
182  RCL O
183  /
184  RCL 03
185  RCL 05
186  -
187  /
188  RCL 05
189  RCL 07
190  -
191  /
192  +
193  RCL 01
194  RCL 07
195  -
196  *
197  RCL N
198  1/X
199  RCL O
200  1/X
201  +
202  RCL 02
203  *
204  +
205  RCL 03
206  RCL 07
207  -
208  1/X
209  RCL 05
210  RCL 07
211  -
212  1/X
213  +
214  RCL 08
215  *
216  +
217  RCL 01
218  RCL 07
219  -
220  *
221  RCL 02
222  RCL 08
223  +
224  3
225  *
226  -
227  RCL 01
228  RCL 07
229  -
230  *
231  ST+ M
232  XEQ 00
233  XEQ 00
234  XEQ 00
235  ISG 00
236  GTO 03                 a three-byte GTO
237  REGSWAP
238  RCL 00
239  INT
240  DSE X
241  STO 00
242  X<> M
243  6
244  /
245  CLA
246  END

( 302 bytes / SIZE 2n+1 )
 
 

STACK	INPUT	OUTPUT
X	/	§x1xn  f(x).dx

 where n = the number of points.  ( n > 3 )

Example:     Calculate  §₁⁸  f(x).dx   from the following values
 

x	1	2.4	4	5.2	7	8
f(x)	1	4	6	5	4	2

Store these 12 numbers into registers   R01  R03  R05  R07  R09  R11   ( x )
                                                          R02  R04  R06  R08  R10  R12   ( f(x) )    respectively

-There are 6 points,     6   STO 00    XEQ "ITG3"    >>>>    §₁⁸  f(x).dx  ~  30.2135   ( in 14 seconds )
 

        b-4) Natural Cubic Spline
 

-In the cubic spline integration, we also connect arcs of cubic polynomials,
  but in such a way that the first derivative y' and the second derivative y" are continuous over the whole interval  [ x₁ , x_n ]
-The spline is called "natural" if we set  y"₁ = y"_n = 0

-This leads to a tridiagonal linear system of  (n-2)  equations in  (n-2)  unknowns.
 

  (1/6).( x_k - x_k-1 ).y"_k-1 + (1/3).( x_k+1 - x_k-1 ).y"_k + (1/6).( x_k+1 - x_k ).y"_k+1 =  ( y_k+1 - y_k )/( x_k+1 - x_k ) - ( y_k - y_k-1 )/( x_k - x_k-1 )

-Then    §_xk^xk+1  y.dx = ( x_k+1 - x_k ).( y_k+1 + y_k )/2 - ( x_k+1 - x_k )³.( y"_k+1 + y"_k )/24
 

Data Registers:       •  R00 = n = number of points                                  ( Registers R00 thru R2n  are to be initialized before executing "NCSI" )

                                  •  R01 = x₁     •  R03 = x₂      .......  •  R2n-1 = x_n             R2n+1 thru R6n-7: temp
                                  •  R02 = y₁     •  R04 = y₂      .......   •  R2n   = y_n

  >>>>   If n > 2 , when the program stops:  R5n-8 = y"₁ = 0 ,  R5n-7 = y"₂ , ............... ,  R6n-10 = y"_n-1 , R6n-9 = y"_n = 0

Flags:  /
Subroutine:   none if  n < 4
                        otherwise, "3DLS"  "Tridiagonal Linear Systems"  ( cf "Linear and non-Linear Systems for the HP-41" )
 

  01  LBL "NCSI"
  02  RCL 00
  03  2
  04  X#Y?
  05  GTO 00
  06  RCL 02                         Lines 06 to 14 simply apply the trapezoidal rule if n = 2
  07  RCL 04
  08  +
  09  RCL 03
  10  RCL 01
  11  -
  12  *
  13  2
  14  /
  15  RTN
  16  LBL 00
  17  3
  18  STO M
  19  RCL 00
  20  ST+ X
  21  1
  22  +
  23  RCL 00
  24  5
  25  *
  26  7
  27  -
  28  STO O
  29  DSE X
  30   E3
  31  /
  32  +
  33  STO N
  34  RCL 05
  35  RCL 01
  36  -
  37  ST+ X
  38  STO IND N
  39  RCL 06
  40  RCL 04
  41  -
  42  RCL 05
  43  RCL 03
  44  -
  45  /
  46  RCL 04
  47  RCL 02
  48  -
  49  RCL 03
  50  RCL 01
  51  -
  52  /
  53  -
  54  6
  55  *
  56  STO IND O
  57  LBL 01
  58  ISG N
  59  FS? 30
  60  GTO 02
  61  2
  62  RCL M
  63  ST+ Y
  64  4
  65  +
  66  RCL IND Y
  67  RCL IND M
  68  -
  69  STO P                      ( synthetic )
  70  STO IND N
  71  ISG N
  72  STO IND N
  73  CLX
  74  RCL IND Y
  75  LASTX
  76  -
  77  STO Q
  78  ST+ X
  79  ISG N
  80  STO IND N
  81  ISG M
  82  CLX
  83  CLX
  84  SIGN
  85  -
  86  RCL IND M
  87  RCL IND Y
  88  -
  89  RCL P
  90  ST- Q
  91  /
  92  ISG O
  93  CLX
  94  STO IND O
  95  X<> Q
  96  RCL Y
  97  2
  98  +
  99  RDN
100  RCL IND T
101  RCL IND Z
102  -
103  X<>Y
104  /
105  ST+ IND O
106  6
107  ST* IND O
108  SIGN
109  ST+ M
110  GTO 01
111  LBL 02
112  RCL 00
113  3
114  X#Y?
115  GTO 03
116  CLX                         If n = 3 the "system" has only one equation which is solved directly.
117  STO 09
118  RCL 07
119  ST/ 08
120  7.009
121  GTO 04
122  LBL 03
123  SIGN
124  ST+ O
125  CLX
126  STO IND O
127  RCL M
128  4
129  +
130  STO Y
131  RCL 00
132  LASTX
133  *
134  +
135  11
136  -
137   E3
138  /
139  +
140  XEQ "3DLS"
141  LASTX                    R00 is disturbed by "3DLS"
142  2                              but its content is now restored
143  +                              R00 = n
144  STO 00
145  CLX
146  .999
147  -
148  LBL 04
149  CLA
150  STO N
151  CLX
152  STO IND N
153  RCL 00
154  ST+ X
155   E3
156  /
157  3
158  +
159  STO O
160  LBL 05
161  RCL O
162  2
163  -
164  RCL IND O
165  RCL IND Y
166  -
167  RCL IND N
168  ISG N
169  RCL IND N
170  +
171  RCL Y
172  X^2
173  *
174  ISG Z
175  RCL IND Z
176  ISG O
177  RCL IND O
178  +
179  12
180  *
181  -
182  *
183  ST- M
184  ISG O
185  GTO 05
186  RCL M
187  24
188  /
189  CLA
190  END

( 283 bytes / SIZE 6n-8 )
 
 

STACK	INPUT	OUTPUT
X	/	§x1xn  f(x).dx

 where n = the number of points.

Example:     Calculate  §₁⁸  f(x).dx   from the following values
 

x	1	2.4	4	5.2	7	8
f(x)	1	4	6	5	4	2

Store these 12 numbers into registers   R01  R03  R05  R07  R09  R11   ( x )
                                                          R02  R04  R06  R08  R10  R12   ( f(x) )    respectively

-SIZE 028
-There are 6 points,     6   STO 00    XEQ "NCSI"    >>>>    §₁⁸  f(x).dx  =  29.99938860   ( in 21 seconds )

-And we have  R22 = y"₁ = 0 ,  R23 = y"₂ = -0.237729622 , R24 = y"₃ = -2.456728203 ,
                                                 R25 = y"₄ = 1.365037775 ,  R26 = y"₅ = -1.986381189 ,  R27 = y"₆ = 0
 

        b-5) Integration of the Lagrange Polynomial
 

 n data points are given and stored in registers  R01 thru R2n
 

x1	x2	...........	xn
y1	y2	...........	yn

-The Lagrange polynomial  is the unique polynomial L(x) of degree < n  such that:     L(x_i) = y_i       i = 1 , 2 , ..... , n
-The program below calculates its coefficients c₀ , c₁ , ........... , c_n-1  when L(x) is written:

    L(x) = c₀ + c₁ ( x-x₁ ) + c₂ ( x-x₁ )² + .................. + c_n-1 ( x-x₁ )^n-1

-First, x₁ is subtracted from all other x_i 's
-We have obviously  c₀ = y₁  and  the other y_i 's  are replaced by  ( y_i - y₁ ) / ( x_i - x₁ )
-Then, we find the value of the new Lagrange polynomial at zero, thus producing c₁ which is stored in register R04
-The process is repeated until  c_n-1 is computed and stored in register R2n
-Finally, lines 42 to 57 evaluate  §_x1^xn  L(x).dx
 

Data Registers:       •  R00 = n = number of points                                  ( Registers R00 thru R2n  are to be initialized before executing "LGI" )

                                  •  R01 = x₁     •  R03 = x₂      .......  •  R2n-1 = x_n
                                  •  R02 = y₁     •  R04 = y₂      .......   •  R2n   = y_n

  >>>>    when the program stops,  R02 = c₀ = y₁ ,  R04 = c₁ , ............... ,  R2n = c_n-1

Flags:  /
Subroutine:   "LAGR"  ( cf "Lagrange's Interpolation Formula for the HP-41" )
 

01  LBL "LGI"
02  RCL 00
03  ST+ X
04  STO Y
05  RCL 01
06  LBL 00
07  DSE Z
08  ST- IND Z
09  DSE Z
10  GTO 00
11  CLX
12   E3
13  /
14  ISG X
15  STO 01
16  RCL 02
17  GTO 02
18  LBL 01
19  2
20  ST+ 01
21  RCL 01
22  0
23  XEQ "LAGR"
24  LBL 02
25  ISG Y
26  STO IND Y
27  ISG Y
28  FS? 30
29  GTO 04
30  LBL 03
31  RCL IND Y
32  ISG Z
33  X<>Y
34  ST- IND Z
35  X<>Y
36  ST/ IND Z
37  RDN
38  ISG Y
39  GTO 03
40  GTO 01
41  LBL 04
42  RCL 00
43  ST+ X
44  0
45  LBL 05
46  RCL IND Y
47  RCL 00
48  /
49  +
50  RCL IND 01
51  *
52  2
53  ST- Z
54  RDN
55  DSE 00
56  GTO 05
57  END

( 98 bytes / SIZE 2n+1 )
 
 

STACK	INPUT	OUTPUT
X	/	§x1xn  L(x).dx

Example:     Calculate  §₁⁸  L(x).dx   from the following values
 

x	1	2.4	4	5.2	7	8
y	1	4	6	5	4	2

Store these 12 numbers into registers   R01  R03  R05  R07  R09  R11   ( x )
                                                          R02  R04  R06  R08  R10  R12    ( y )    respectively
 

-There are 6 points:     6   STO 00      XEQ "LGI"     >>>>    §₁⁸  L(x).dx  =   29.61789480  ( in 43 seconds )
 ( this value was already found more laboriously in paragraph 2°) b-2) above )

-We have also the coefficients of the Lagrange polynomial ( as a sum of powers of ( x-x₁ ) )  in registers  R02  R04  .............  R12
-Here, x₁ = 1 and we find:

  L(x) =  1 - 0.362103178 ( x-1 ) + 3.623795356 ( x-1 )² - 1.661873944 ( x-1 )³ + 0.272598127 ( x-1 )⁴ - 0.015381483 ( x-1 )⁵
 

Remarks:

-For large n-values, the coefficients of the Lagrange polynomial are often great in magnitude
  and of alternate signs, thus leading to inaccurate results.
-Moreover, the coefficients themselves are not obtained very accurately if n is too great.

-For instance, with the following data:
 

x	1	3	5	7	9	11	13	15	17	19	21	23	25	27	29	31	33	35	37	39
y	2	4	6	8	10	12	14	16	18	20	22	24	26	28	30	32	34	36	38	40

 "LGI"  yields  797.9971774  wheras the exact result is 798

Execution time t:
 
 

n	t
4	15s
6	43s
10	2mn51s
20	21mn13s

 

        c) Double Integrals  ( Equally Spaced Arguments )
 

-We want to evaluate  §_x1^xn  §_y1^ym  f(x,y).dx.dy    assuming  n & m are odd integers  ( n > 1 , m > 1 )
-The following program uses Simpson's rule in the direction of  x-axis and y-axis

-  x_i  must be  equally spaced:  x_i+1 - x_i = h     for i = 1,2,...,n-1
-  y_j  must be  equally spaced:  y_j+1 - y_j = k    for j = 1,2,...,m-1

Data Registers:       •  R00   = h.k                                              ( Registers R00 thru Rn.m  are to be initialized before executing "INT2" )

                                  •  R01   = f(x₁,y₁)  •  Rn+1 = f(x₁,y₂)  .......  •  Rnm-n+1 = f(x₁,y_m)
                                  •  R02   = f(x₂,y₁)  •  Rn+2 = f(x₂,y₂)  .......  •  Rnm-n+2 = f(x₂,y_m)
                                     .......................................................................................
                                  •  Rnn   = f(x_n,y₁)  •  R2n = f(x_n,y₂)    .......  •  Rnm = f(x_n,y_m)
Flags:  /
Subroutines:  /

-Synthetic registers M N O may be replaced by any unused data registers, for instance  R97  R98  R99 if  n.m < 97
 

01  LBL "INT2"
02  STO N
03  X<>Y
04  STO M
05  *
06  STO O
07  CLX
08  LBL 01
09  RCL O
10  RCL N
11  MOD
12  2
13  X>Y?
14  CLX
15  MOD
16  LASTX
17  X<>Y
18  -
19  ST+ X
20  X<=0?
21  SIGN
22  ABS
23  RCL O
24  1
25  -
26  RCL N
27  /
28  INT
29  1
30  +
31  RCL M
32  MOD
33  2
34  X>Y?
35  CLX
36  MOD
37  LASTX
38  X<>Y
39  -
40  ST+ X
41  X<=0?
42  SIGN
43  ABS
44  *
45  RCL IND O
46  *
47  +
48  DSE O
49  GTO 01
50  RCL 00
51  *
52  9
53  /
54  CLA
55  END

( 77 bytes / SIZE n.m+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	m	/
X	n	§x1xn §y1ym f(x,y).dx.dy

Example:     Calculate   §₂⁶ §₁⁵  f(x,y).dx.dy   from the following f(x,y) values
 

x\y	1	2	3	4	5
2	3	4	7	6	3
4	1	2	4	5	3
6	4	1	3	4	6

                    3  4  7  6  3                                   R01  R04  R07  R10  R13
 We store     1  2  4  5  3        into registers        R02  R05  R08  R11  R14      respectively
                    4  1  3  4  6                                   R03  R06  R09  R12  R15

-Here, we have  h = 2 and  k = 1   whence  1*2 = 2  STO 00

   n = 3  and  m = 5    whence    5  ENTER^
                                                 3  XEQ "INT2"  >>>>   §₂⁶ §₁⁵  f(x,y).dx.dy  ~  56.8889

-Perfect accuracy is achieved if f is a polynomial and  deg(f) < 4
 

        d) Triple Integrals  ( Equally Spaced Arguments )
 

-Now we evaluate  §_x1^xn  §_y1^ym  §_z1^zp  f(x,y,z).dx.dy.dz    assuming  n , m , p are odd integers  ( n > 1 , m > 1 , p > 1 )
-Simpson's rule in the direction of  x- , y- and z-axis is used.

-  x_i  must be  equally spaced:  x_i+1 - x_i = h         for i = 1,2,...,n-1
-  y_j  must be  equally spaced:  y_j+1 - y_j = h'       for j = 1,2,...,m-1
-  z_k  must be  equally spaced:  z_k+1 - y_k = h''    for j = 1,2,...,p-1

Data Registers:       •  R00   = h.h'.h''                                                ( Registers R00 thru Rn.m.p  are to be initialized before executing "INT3" )

                                  •  R01   = f(x₁,y₁,z₁)  •  Rn+1 = f(x₁,y₂,z₁)  .......  •  Rnm-n+1 = f(x₁,y_m,z₁)
                                  •  R02   = f(x₂,y₁,z₁)  •  Rn+2 = f(x₂,y₂,z₁)  .......  •  Rnm-n+2 = f(x₂,y_m,z₁)
                                     ......................................................................................................
                                  •  Rn     = f(x_n,y₁,z₁)  •  R2n = f(x_n,y₂,z₁)    .......  •  Rnm = f(x_n,y_m,z₁)

                                        •  Rnm+1   = f(x₁,y₁,z₂)  •  Rnm+n+1 = f(x₁,y₂,z₂)  .......  •  R2nm-n+1 = f(x₁,y_m,z₂)
                                        •  Rnm+2   = f(x₂,y₁,z₂)  •  Rnm+n+2 = f(x₂,y₂,z₂)  .......  •  R2nm-n+2 = f(x₂,y_m,z₂)
                                           ......................................................................................................
                                        •  Rnm+n   = f(x_n,y₁,z₂)  •  Rnm+2n = f(x_n,y₂,z₂)    .......  •  R2nm = f(x_n,y_m,z₂)

                                               .............................................................................................................
                                                .............................................................................................................
                                                 .............................................................................................................

                                                   •  Rnm(p-1)+1   = f(x₁,y₁,z_p)  •  Rnm(p-1)+n+1 = f(x₁,y₂,z_p)  .......  •  Rnmp-n+1 = f(x₁,y_m,z_p)
                                                   •  Rnm(p-1)+2   = f(x₂,y₁,z_p)  •  Rnm(p-1)+n+2 = f(x₂,y₂,z_p)  .......  •  Rnmp-n+2 = f(x₂,y_m,z_p)
                                                      ......................................................................................................
                                                   •  Rnm(p-1)+n   = f(x_n,y₁,z_p)  •  Rnm(p-1)+2n = f(x_n,y₂,z_p)    .......  •  Rnmp = f(x_n,y_m,z_p)
Flags:  /
Subroutines:  /

-Synthetic registers M N O P Q may be replaced by any unused data registers, for instance  R95 R95 R97  R98  R99 if  n.m.p < 95
 

01  LBL "INT3"
02  STO N
03  X<>Y
04  STO M
05  *
06  STO Q
07  X<>Y
08  STO O
09  *
10  STO P          ( synthetic )
11  CLX
12  LBL 01
13  RCL P
14  RCL N
15  MOD
16  ENTER^
17  SIGN
18  ST- P
19  ST+ X
20  X>Y?
21  CLX
22  MOD
23  LASTX
24  X<>Y
25  -
26  ST+ X
27  X<=0?
28  SIGN
29  ABS
30  RCL P
31  RCL Q
32  /
33  INT
34  ENTER^
35  SIGN
36  +
37  RCL O
38  MOD
39  ENTER^
40  SIGN
41  ST+ X
42  X>Y?
43  CLX
44  MOD
45  LASTX
46  X<>Y
47  -
48  ST+ X
49  X<=0?
50  SIGN
51  ABS
52  *
53  RCL P
54  RCL Q
55  MOD
56  RCL N
57  /
58  INT
59  ENTER^
60  SIGN
61  ST+ P
62  +
63  RCL M
64  MOD
65  ENTER^
66  SIGN
67  ST+ X
68  X>Y?
69  CLX
70  MOD
71  LASTX
72  X<>Y
73  -
74  ST+ X
75  X<=0?
76  SIGN
77  ABS
78  *
79  RCL IND P
80  *
81  +
82  DSE P
83  GTO 01
84  RCL 00
85  *
86  27
87  /
88  CLA
89  END

( 124 bytes / SIZE nmp+1 )
 
 

STACK	INPUTS	OUTPUTS
Z	p	/
Y	m	/
X	n	§x1xn§y1ym§z1zp f(x,y,z)dx.dy.dz

Example:     Evaluate   §₁³ §₁⁵ §₁⁷  f(x,y,z).dx.dy.dz   from the following values

  f(1,1,1) =  4      f(1,3,1) =  6     f(1,5,1) =  8                                                     R01   R04   R07
  f(2,1,1) =  7      f(2,3,1) =  9     f(2,5,1) = 11           to be stored in registers      R02   R05   R08     respectively
  f(3,1,1) = 10     f(3,3,1) = 12    f(3,5,1) = 14                                                    R03   R06   R09

       f(1,1,4) =  64       f(1,3,4) =  96      f(1,5,4) = 128                                                    R10   R13   R16
       f(2,1,4) = 112      f(2,3,4) = 144     f(2,5,4) = 176           to be stored in registers      R11   R14   R17     respectively
       f(3,1,4) = 160      f(3,3,4) = 192     f(3,5,4) = 224                                                    R12   R15   R18

            f(1,1,7) =  196      f(1,3,7) = 294     f(1,5,7) = 392                                                    R19   R22   R25
            f(2,1,7) =  343      f(2,3,7) = 441     f(2,5,7) = 539           to be stored in registers      R20   R23   R26     respectively
            f(3,1,7) =  490      f(3,3,7) = 588     f(3,5,7) = 686                                                    R21   R24   R27

   h.h'.h'' = 1*2*3 = 6     6  STO 00

  n = m = p = 3   whence   3  ENTER^   ENTER^   XEQ "INT3"   >>>>   §₁³ §₁⁵ §₁⁷  f(x,y,z).dx.dy.dz  =  8208

-These values were actually computed from  f(x,y,z) =  (3x+y).z²   and the result is perfectly accurate since f is a polynomial and deg(f) < 4.
 
 

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 C++" - Cambridge University Press - ISBN  0-521-75033-4
 
 

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