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
f) 2 Other Formulae
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: §1°) f)
-For continuous data, see also
"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
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
-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
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
f) 2 Other Formulae
-Here is another formula to compute §ab f(x) dx without computing f(a) & f(b)
§-11 f(x).dx ~ { 275 [ f(-4/5) + f(4/5) ] + 100 [ f(-2/5) + f(2/5) ] + 402 f(0) } / 576 exact formula if f is a polynomial with deg(f) < 6
Data Registers: • R00 = Function name ( Registers R00 thru R03 are to be initialized before executing "IG5" )
• R01 = a
• R02 = b
• R03 = n = number of subintervals over which the formula is applied.
R04 = §ab f(x).dx
R05 thru R09: temp
Flags: /
Subroutine: A program which computes
f(x) assuming x is in X-register upon entry.
| 01 LBL "IG5" 02 RCL 02 03 RCL 01 04 - 05 RCL 03 06 STO 05 07 / 08 STO 06 09 5 10 / 11 STO 08 12 CLX |
13 STO 04 14 LBL 01 15 RCL 05 16 .5 17 - 18 RCL 06 19 * 20 RCL 01 21 + 22 STO 07 23 RCL 08 24 ST+ X |
25 STO 09 26 - 27 XEQ IND 00 28 X<> 09 29 RCL 07 30 + 31 XEQ IND 00 32 RCL 09 33 + 34 2.75 35 * 36 ST+ 04 |
37 RCL 07 38 RCL 08 39 - 40 XEQ IND 00 41 ST+ 04 42 RCL 07 43 RCL 08 44 + 45 XEQ IND 00 46 ST+ 04 47 RCL 07 48 XEQ IND 00 |
49 4.02 50 * 51 ST+ 04 52 DSE 05 53 GTO 01 54 RCL 04 55 RCL 06 56 * 57 11.52 58 / 59 STO 04 60 END |
( 92 bytes / SIZE 010 )
| STACK | INPUT | OUTPUT |
| X | / | §ab f(x).dx |
Example: Evaluate §13
e-x^2 dx
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 "IG5" >>>> 0.139383245
n = 8 8 STO 03 R/S >>>> 0.139383216
Notes:
-The exact value is 0.139383215 (rounded to 9 decimals )
-Here is yet another formula to compute §ab f(x) dx without computing f(a) & f(b)
§-11 f(x).dx ~ { 24745 [ f(-6/7) + f(6/7) ] + 882 [ f(-4/7) + f(4/7) ] + 56007 [ f(-2/7) + f(2/7) ] - 25028 f(0) } / 69120 exact formula if f is a polynomial with deg(f) < 8
Data Registers: • R00 = Function name ( Registers R00 thru R03 are to be initialized before executing "IG7" )
• R01 = a
• R02 = b
• R03 = n = number of subintervals over which the formula is applied.
R04 = §ab f(x).dx
R05 thru R09: temp
Flags: /
Subroutine: A program which computes
f(x) assuming x is in X-register upon entry.
| 01 LBL "IG7" 02 RCL 02 03 RCL 01 04 - 05 RCL 03 06 STO 05 07 / 08 STO 06 09 7 10 / 11 STO 08 12 CLX 13 STO 04 14 LBL 01 15 RCL 05 16 .5 |
17 - 18 RCL 06 19 * 20 RCL 01 21 + 22 STO 07 23 RCL 08 24 3 25 * 26 STO 09 27 - 28 XEQ IND 00 29 X<> 09 30 RCL 07 31 + 32 XEQ IND 00 |
33 RCL 09 34 + 35 24745 36 * 37 ST+ 04 38 RCL 07 39 RCL 08 40 ST+ X 41 STO 09 42 - 43 XEQ IND 00 44 X<> 09 45 RCL 07 46 + 47 XEQ IND 00 48 RCL 09 |
49 + 50 882 51 * 52 ST+ 04 53 RCL 07 54 RCL 08 55 - 56 XEQ IND 00 57 STO 09 58 RCL 07 59 RCL 08 60 + 61 XEQ IND 00 62 RCL 09 63 + 64 56007 |
65 * 66 ST+ 04 67 RCL 07 68 XEQ IND 00 69 25028 70 * 71 ST- 04 72 DSE 05 73 GTO 01 74 RCL 04 75 RCL 06 76 * 77 138240 78 / 79 STO 04 80 END |
( 124 bytes / SIZE 010 )
| STACK | INPUT | OUTPUT |
| X | / | §ab f(x).dx |
Example: Evaluate §13
e-x^2 dx
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 "IG7" >>>> 0.139383215
n = 8 8 STO 03 R/S >>>> 0.139383215
-The exact value is 0.139383215 (rounded to 9 decimals )
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
)
-If n is always odd, the program becomes smaller:
| 01 LBL "SIMP" 02 RCL IND X 03 DSE Y 04 LBL 01 05 RCL IND Y 06 4 07 * 08 + 09 DSE Y 10 RCL IND Y 11 ST+ X 12 + 13 DSE Y 14 GTO 01 15 RCL 01 16 - 17 RCL 00 18 * 19 3 20 / 21 END |
( 38 bytes / SIZE n+1 )
| STACK | INPUT | OUTPUT |
| X | n | §x1xn f(x).dx |
where n = the number of points , n > 2
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.8570 | 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!
-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
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