# hp41programs

NumInt Numerical Integration for the HP-41

Overview

1°) Continuous Data

a)  Gauss-Legendre 2-point & 3-point Formulae
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
b-2) Connecting Parabolic Segments ( 2 programs )
b-3) Connecting Cubic Segments
b-4) Natural Cubic Spline
b-5) Integration of the Lagrange Polynomial

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

Latest Update:   Gauss-Legendre 2-point formula

"Double Integrals for the HP-41"
"Triple Integrals for the HP-41"
"Multiple Integration for the HP-41"

-Integrals are denoted §      §ab f(x).dx

1°) Continuous Data

a) Gauss-Legendre 2-Point & 3-Point  Formulae

-Gauss-Legendre formulas are defined by  §-11 f(x).dx  ~  w1.x1 + w2.x2 + ....... + wn.xn
where wi and xi  are choosen so that the formula produces perfect accuracy when  f(x) is a polynomial of degree < 2n
-The  xi  are the zeros of the Legendre polynomials and they can't be easily computed for large n-values.

-The 2-point Gauss-Legendre formula is   §-11 f(x).dx  ~  f(-1/sqrt(3)) + f(+1/sqrt(3))

-The 3-point Gauss-Legendre formula is   §-11 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 "GL2" )

•  R01 = a
•  R02 = b
•  R03 = n = number of subintervals over which the formula is applied.          R04 =  §ab 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 "GL2"  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 3  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 RCL 08            27 +  28 XEQ IND 00  29 ST+ 04  30 RCL 06  31 ST+ 07  32 DSE 05 33 GTO 01  34 RCL 04            35 *  36 2  37 /  38 STO 04  39 END

( 55 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS X / §ab f(x).dx

Example:
Evaluate  §13  e-x^2 dx

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

n =  4     4   STO 03  XEQ "GL2"  >>>>  0.139404286
n = 16   16  STO 03       R/S          >>>>  0.139383300

Notes:

-The exact value is  0.139383215  (rounded to 9 decimals )

-"GL2" is one of the shortest program to compute numerical integrations
-But 3-point formula will be more accurate:

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 =  §ab 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 )

 STACK INPUTS OUTPUTS X / §ab f(x).dx

Example:
Evaluate  §13  e-x^2 dx

-With the same subroutine,

"T"  ASTO 00
1    STO 01
3    STO 02

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 wi and xi  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 )

 STACK INPUTS OUTPUTS X / §ab f(x).dx

Example:     Evaluate  §03  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

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"
replace line 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:    §ab  ((x-a)(x-b))-1/2  f(x).dx  ~  (pi/n) [ f(x1) + f(x2) + ..... + f(xn) ]
where   xi = (a+b)/2 + ((b-a)/2).cos((2i-1)pi/(2n))
-Unlike "GL3" and "GL16", "GC" determines the coefficients wi and xi 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       14  ASIN 15  RCL 05 16  ST+ X 17  DSE X    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 )

 STACK INPUTS OUTPUTS X / §ab f(x).dx

Example:     Compute  §13  ex(-x2+4.x-3)-1/2 dx  =  §13  ex((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

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:    §ab  f(x).cos(k.x).dx   and   §ab  f(x).sin(k.x).dx

§ab  f(x).cos(k.x).dx  ~  h.[ A(k.h) (  f(x2n).sin(k.x2n)  -  f(x0).sin(k.x0) ) + B(k.h) C2n + C(k.h) C2n-1 ]                 a = x0 ;  b = x2n ;  h = (b-a)/(2n)
§ab  f(x).sin(k.x).dx   ~  h.[ A(k.h) ( -f(x2n).cos(k.x2n) + f(x0).cos(k.x0) ) + B(k.h) S2n + C(k.h) S2n-1 ]

where   C2n = f(x0).cos(k.x0) + f(x2).cos(k.x2) + ......... + f(x2n).cos(k.x2n)  -  (1/2) ( f(x2n).cos(k.x2n) + f(x0).cos(k.x0) )
S2n = f(x0).sin(k.x0) + f(x2).sin(k.x2) + ......... + f(x2n).sin(k.x2n)  -  (1/2) ( f(x2n).sin(k.x2n) + f(x0).sin(k.x0) )

C2n-1 = f(x1).cos(k.x1) + f(x3).cos(k.x3) + ......... + f(x2n-1).cos(k.x2n-1)
S2n-1 = f(x1).sin(k.x1) + f(x3).sin(k.x3) + ......... + f(x2n-1).sin(k.x2n-1)

and     A(µ) = 1/µ + (sin(2µ))/(2µ2) - 2(sin2µ)/µ3
B(µ) =  2(1+cos2µ)/µ2 - 2(sin(2µ))/µ3
C(µ) =  4(sinµ)/µ3 - 4(cosµ)/µ2

-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

>>>  At the end,    R04 =   §ab  f(x).cos(k.x).dx  = X-register  ,  R05 =   §ab  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 = §16  Ln(x).cos(10x).dx    and   J = §16  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:  x2 + y2 = R2     (  ds = ( dx2+dy2)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 )

 STACK INPUTS OUTPUTS X / §(C)  f(x,y).ds

Example:   Evaluate    §(C)  Ln(3+x.y).ds    where  (C):  x2 + y2 = 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  §x1xn  f(x).dx  but we only know   f(x1) , f(x2) , ....... , f(xn)  where the  xi  are equally spaced     xi+1 - xi = h    for i = 1,2,...,n-1

-If n is odd, we use the formula:    §x1x3  f(x).dx  ~  h[f(x1)+4f(x2)+f(x3)]/3    over the intervals  [x1;x3] ,  [x3;x5] , ..... ,  [xn-2;xn]
-If n is even,  we use  §x1x4  f(x).dx  ~  3h[f(x1)+3f(x2)+3f(x3)+f(x4)]/8  and the same formula as above over  [x4;x6] , ..... ,  [xn-2;xn]

-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(x1)  •  R02 = f(x2)  .......  •  Rnn = f(xn)
Flags:  /
Subroutines:  /

-Synthetic register M may of course be replaced by any unused data register, for instance R99 ( provided n < 99 )
-Lines 21 to 36 are only useful when n is even

 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      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  §0pi/2  f(x).dx  and  §05pi/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 = §0pi/2  sin(x).dx = 1 )
6         R/S          >>>>   0.7412102   ( exact result = §05pi/12  sin(x).dx = 0.7411810 )

a-2) Other Newton-Cotes Formulae

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

Formula:   §x1x7  f(x).dx  ~  h[41f(x1)+216f(x2)+27f(x3)+272f(x4)+27f(x5)+216f(x6)+41f(x7)]/140

Data Registers:       •  R00 = h                                                               ( Registers R00 thru Rnn  are to be initialized before executing "NC6" )

•  R01 = f(x1)  •  R02 = f(x2)  .......  •  Rnn = f(xn)
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  §0pi/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 ):

§x1x10  f(x).dx  ~  9h[2857(f(x1)+f(x10))+15741(f(x2)+f(x9))+1080(f(x3)+f(x8))+19344(f(x4)+f(x7))+5778(f(x5)+f(x6))]/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  §x1xn  f(x).dx
-We simply add the areas of n trapezoids:  Sumi=1,2,...,n-1   (xi+1-xi).[ f(xi+1)+f(xi) ]/2

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

•  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
•  R02 = f(x1)  •  R04 = f(x2)  .......  •  R2n   = f(xn)
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  §18  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"    >>>>    §18  f(x).dx  ~  29.2

b-2) Connecting Parabolic Segments  ( 2 programs )

-More accurate results may be obtained if we use several connected parabolic segments to compute  §x1xn  f(x).dx
-However, if n is even,  §x1x2  f(x).dx  is calculated by a polynomial of degree 3
-Thus, "ITG" yields the same results as "SIMP"  when the  xi'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 = x1     •  R03 = x2      .......  •  R2n-1 = xn
•  R02 = f(x1)  •  R04 = f(x2)  .......  •  R2n   = f(xn)
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
-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:

RCL 06
RCL 05
RCL 01
-
RCL 03
RCL 05
-
*
/
RCL 03
RCL 01
-
STO O
*
RCL 04
RCL 05
RCL 03
-
/
+
RCL 02
RCL 05
RCL 01
-
/
-
RCL O
*
RCL 02
RCL 04
+
3
*
+
ST* O

 01  LBL "ITG"     02  CLA   03  RCL 00   04  2   05  -   06  ST+ X   07  STO M   08  4   09  MOD   10  X#0?     11  GTO 01   12  RCL 01    13  RCL 03     14  +     15  2     16  /       17  STO O    18  RCL 05   19  -         20  STO N     21  RCL 07    22  RCL 01   23  -        24  /       25  RCL 07   26  RCL 03    27  -           28  /       29  RCL 08 30  *         31  RCL O    32  RCL 07   33  -        34  ST* N     35  RCL 05          36  RCL 01   37  -       38  /       39  RCL 05    40  RCL 03    41  -         42  /         43  RCL 06   44  *       45  -       46  RCL 05   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     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  §17  f(x).dx  and  §18  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"  >>>>  §17  f(x).dx  ~  26.4226
-Six points:     6   STO 00        R/S         >>>>  §18  f(x).dx  ~  30.5339

-With the modification listed above the listing, it gives  §18  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.857 5.3453 6 5.2069 4.3568 4 2

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

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

§x1x8  f(x).dx  ~  7h[751f(x1)+3577f(x2)+1323f(x3)+2989f(x4)+2989f(x5)+1323f(x6)+3577f(x7)+751f(x8)]/17280
§x1x8  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!

Note:

-We can simplify the program if n is odd ( at least 3 points )

Data Registers:          R00 to R03: temp        n odd and > 2               ( Registers Rbb thru Ree  are to be initialized before executing "ITG" )

•  Rbb = x1        •  Rb+2 = x2      .......  •  Ree-1 = xn               bbb > 003
•  Rb+1 = f(x1)  •  Rb+3 = f(x2)   .......  •  Ree   = f(xn)
Flags:  /
Subroutines:  /

 01 LBL "ITG"  02 2  03 +  04 STO 01  05 CLX  06 STO 00            07 LBL 01  08 2  09 CHS  10 2  11 RCL 01  12 INT 13 ST+ Y  14 ST+ Z  15 RCL IND Z  16 RCL IND Y  17 -  18 STO 03  19 RCL IND Z  20 LASTX  21 -  22 STO 02            23 /  24 3 25 ST+ 01  26 ST+ Z  27 DSE X  28 +  29 RCL IND Y  30 *  31 RCL 03            32 RCL 02  33 -  34 X^2  35 RCL 02  36 / 37 RCL 03  38 /  39 DSE Z  40 DSE Z  41 RCL IND Z  42 *  43 -  44 RCL 02            45 RCL 03  46 ST- 02  47 /  48 2  49 ST- T 50 +  51 RCL IND Z  52 *  53 +  54 RCL 02  55 *  56 ST+ 00  57 ISG 01  58 GTO 01  59 RCL 00            60 6  61 /  62 END

( 88 bytes / SIZE ee+1 )

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

where n = the number of points.  ( n odd > 2 )     and    bbb >  003

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

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

-If you store for instance these 10 numbers into registers  R11 thru R20

11.020   XEQ "ITG"  >>>>    §17  f(x).dx  =  26.42261905

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

-The following program uses several segments of cubic curves to evaluate  §x1xn  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 = x1     •  R03 = x2      .......  •  R2n-1 = xn
•  R02 = f(x1)  •  R04 = f(x2)  .......  •  R2n   = f(xn)
Flags:  /
Subroutines:  /

-Line 236 is a three-byte  GTO 03

 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  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  §18  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"    >>>>    §18  f(x).dx  ~  30.2135   ( in 14 seconds )

b-4) Natural Cubic Spline

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

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

(1/6).( xk - xk-1 ).y"k-1 + (1/3).( xk+1 - xk-1 ).y"k + (1/6).( xk+1 - xk ).y"k+1 =  ( yk+1 - yk )/( xk+1 - xk ) - ( yk - yk-1 )/( xk - xk-1 )

-Then    §xkxk+1  y.dx = ( xk+1 - xk ).( yk+1 + yk )/2 - ( xk+1 - xk )3.( 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 = x1     •  R03 = x2      .......  •  R2n-1 = xn             R2n+1 thru R6n-7: temp
•  R02 = y1     •  R04 = y2      .......   •  R2n   = yn

>>>>   If n > 2 , when the program stops:  R5n-8 = y"1 = 0 ,  R5n-7 = y"2 , ............... ,  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" )

-Lines 06 to 14 simply apply the trapezoidal rule if n = 2
-If n = 3 the "system" has only one equation which is solved directly.

 01  LBL "NCSI"   02  RCL 00   03  2   04  X#Y?   05  GTO 00   06  RCL 02    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     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   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     142  2           143  +     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  §18  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"    >>>>    §18  f(x).dx  =  29.99938860   ( in 21 seconds )

-And we have  R22 = y"1 = 0 ,  R23 = y"2 = -0.237729622 , R24 = y"3 = -2.456728203 ,
R25 = y"4 = 1.365037775 ,  R26 = y"5 = -1.986381189 ,  R27 = y"6 = 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(xi) = yi       i = 1 , 2 , ..... , n
-The program below calculates its coefficients c0 , c1 , ........... , cn-1  when L(x) is written:

L(x) = c0 + c1 ( x-x1 ) + c2 ( x-x1 )2 + .................. + cn-1 ( x-x1 )n-1

-First, x1 is subtracted from all other xi 's
-We have obviously  c0 = y1  and  the other yi 's  are replaced by  ( yi - y1 ) / ( xi - x1 )
-Then, we find the value of the new Lagrange polynomial at zero, thus producing c1 which is stored in register R04
-The process is repeated until  cn-1 is computed and stored in register R2n
-Finally, lines 42 to 57 evaluate  §x1xn  L(x).dx

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

•  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
•  R02 = y1     •  R04 = y2      .......   •  R2n   = yn

>>>>    when the program stops,  R02 = c0 = y1 ,  R04 = c1 , ............... ,  R2n = cn-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  §18  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"     >>>>    §18  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-x1 ) )  in registers  R02  R04  .............  R12
-Here, x1 = 1 and we find:

L(x) =  1 - 0.362103178 ( x-1 ) + 3.623795356 ( x-1 )2 - 1.661873944 ( x-1 )3 + 0.272598127 ( x-1 )4 - 0.015381483 ( x-1 )5

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  §x1xn  §y1ym  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

-  xi  must be  equally spaced:  xi+1 - xi = h     for i = 1,2,...,n-1
-  yj  must be  equally spaced:  yj+1 - yj = 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(x1,y1)  •  Rn+1 = f(x1,y2)  .......  •  Rnm-n+1 = f(x1,ym)
•  R02   = f(x2,y1)  •  Rn+2 = f(x2,y2)  .......  •  Rnm-n+2 = f(x2,ym)
.......................................................................................
•  Rnn   = f(xn,y1)  •  R2n = f(xn,y2)    .......  •  Rnm = f(xn,ym)
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   §26 §15  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
Store     1  2  4  5  3        in        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"  >>>>   §26 §15  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  §x1xn  §y1ym  §z1zp  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.

-  xi  must be  equally spaced:  xi+1 - xi = h         for i = 1,2,...,n-1
-  yj  must be  equally spaced:  yj+1 - yj = h'       for j = 1,2,...,m-1
-  zk  must be  equally spaced:  zk+1 - yk = 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(x1,y1,z1)  •  Rn+1 = f(x1,y2,z1)  .......  •  Rnm-n+1 = f(x1,ym,z1)
•  R02   = f(x2,y1,z1)  •  Rn+2 = f(x2,y2,z1)  .......  •  Rnm-n+2 = f(x2,ym,z1)
......................................................................................................
•  Rn     = f(xn,y1,z1)  •  R2n = f(xn,y2,z1)    .......  •  Rnm = f(xn,ym,z1)

•  Rnm+1   = f(x1,y1,z2)  •  Rnm+n+1 = f(x1,y2,z2)  .......  •  R2nm-n+1 = f(x1,ym,z2)
•  Rnm+2   = f(x2,y1,z2)  •  Rnm+n+2 = f(x2,y2,z2)  .......  •  R2nm-n+2 = f(x2,ym,z2)
......................................................................................................
•  Rnm+n   = f(xn,y1,z2)  •  Rnm+2n = f(xn,y2,z2)    .......  •  R2nm = f(xn,ym,z2)

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

•  Rnm(p-1)+1   = f(x1,y1,zp)  •  Rnm(p-1)+n+1 = f(x1,y2,zp)  .......  •  Rnmp-n+1 = f(x1,ym,zp)
•  Rnm(p-1)+2   = f(x2,y1,zp)  •  Rnm(p-1)+n+2 = f(x2,y2,zp)  .......  •  Rnmp-n+2 = f(x2,ym,zp)
......................................................................................................
•  Rnm(p-1)+n   = f(xn,y1,zp)  •  Rnm(p-1)+2n = f(xn,y2,zp)    .......  •  Rnmp = f(xn,ym,zp)
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  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   §13 §15 §17  f(x,y,z).dx.dy.dz   from the following values

-Store:

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         in        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           in        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             in         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"   >>>>   §13 §15 §17  f(x,y,z).dx.dy.dz  =  8208

-These values were actually computed from  f(x,y,z) =  (3x+y).z2   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