Gauss-Kronrod Integration for the HP-41
Overview
-When we use a Gauss-Legendre m-point formula to evaluate I = §ab
f(x) dx , there is no error estimate.
-To get an error estimate, we can subdivide the interval into 2 parts
and take the difference between the 2 results as an error-estimate.
-Gauss-Kronrod integration adds new abscissas to those already employed
by a Gauss-Legendre m-point formula to build a more accurate formula
-Thus, less function-evaluations are required to get an error estimate.
-The following programs use a 7-point Gauss formula embedded in a 15-point
Kronrod formula.
-So, 15 function evaluations are computed instead of 21.
-Higher order formulae may be found on the web.
Program Listing1
-The interval [a;b] may be divided into n sub-intervals
[ a+k.(b-a)/n ; a+(k+1).(b-a)/n ] k = 0 , 1 , ......
, n-1
- n must be stored in R03
Data Registers: • R00 = function name ( Registers R00 thru R03 & R12 thru R30 are to be initialized before executing "GK7" )
• R01 = a
R04 = Gauss7-Integral
R06 to R11 & R31: temp
• R02 = b
R05 = Kronrod15-Integral
• R03 = n
• R12 = 0.4179591837
• R17 = 0.1903505781
• R22 = 0.1406532597
• R27 = 0.06309209263
• R13 = 0.2094821411
• R18 = 0.4058451514
• R23 = 0.7415311856
• R28 = 0.9491079123
• R14 = 0.2044329401
• R19 = 0.1690047266
• R24 = 0.1047900103
• R29 = 0.02293532201
• R15 = 0.2077849550
• R20 = 0.5860872355
• R25 = 0.8648644234
• R30 = 0.9914553711
• R16 = 0.3818300505
• R21 = 0.2797053915
• R26 = 0.1294849662
Flags: /
Subroutine: A program that returns f(x)
in X-register, assuming x is in X-register upon entry.
| 01 LBL "GK7"
02 CLX 03 STO 04 04 STO 05 05 RCL 02 06 RCL 01 07 STO 08 08 - 09 RCL 03 10 STO 31 11 ST+ X 12 / 13 STO 07 14 LBL 01 15 ST+ 08 |
16 30.013
17 STO 09 18 SIGN 19 STO 06 20 LBL 02 21 RCL 08 22 RCL IND 09 23 RCL 07 24 * 25 STO 10 26 - 27 XEQ IND 00 28 STO 11 29 RCL 08 30 RCL 10 |
31 +
32 XEQ IND 00 33 RCL 11 34 + 35 STO 10 36 DSE 09 37 RCL IND 09 38 * 39 ST+ 05 40 RCL 06 41 CHS 42 STO 06 43 X<0? 44 GTO 02 45 RCL 10 |
46 DSE 09
47 RCL IND 09 48 * 49 ST+ 04 50 LBL 02 51 DSE 09 52 GTO 02 53 RCL 08 54 XEQ IND 00 55 STO 10 56 RCL IND 09 57 * 58 ST+ 05 59 DSE 09 60 CLX |
61 RCL 10
62 RCL IND 09 63 * 64 ST+ 04 65 RCL 07 66 ST+ 08 67 DSE 31 68 GTO 01 69 ST* 04 70 ST* 05 71 RCL 05 72 RCL 04 73 END |
( 112 bytes / SIZE 032 )
| STACK | INPUTS | OUTPUTS |
| Y | / | I15 |
| X | / | I7 |
with I7 = Gauss7-Integral
= R04
and I15 = Kronrod15-Integral
= R05
Example: Evaluate §04 Sin ( x2 ) dx
LBL "T"
X^2
SIN
RTN
END
"T" ASTO 00
XEQ "RAD"
0 STO 01
4 STO 02
-With n = 2 2 STO 03
XEQ "GK7" >>>> I7 =
0.747111478
---Execution time = 33s---
RDN I15 = 0.747133845
-With n = 4 4 STO 03
R/S
>>>> I7 = 0.747133834
---Execution time = 66s---
RDN I15 = 0.747133846
Variant:
-Add the following instructions after line 01
| ASTO 00
STO 02 X<>Y STO 01 |
1
STO 03 LBL 03 2 |
ST* 03
XEQ 00 VIEW 04 RND |
X<>Y
RND X#Y? GTO 03 |
LASTX
CLD RTN LBL 00 |
-Thus, neither R00 nor R01 ..... nor R03 needs to be initialized:
-Place the function name in the alpha "register", a ENTER^ b
XEQ "GK7"
-The loop is executed until 2 rounded approximations are equal
-The accuracy is then controlled by the display settings
-With the same example §04 Sin ( x2 ) dx
alpha T alpha
FIX 9
0 ENTER^
4 XEQ "GK7" >>>> I = 0.747133846
( n = R03 = 16 )
Program Listing2
-Instead of storing the abscissas and weights into R12 to R30, they
may be inserted in the program itself.
-The version below employs this method: less data registers are used,
but the routine occupies 344 bytes instead of 112 !
Data Registers: • R00 = function name ( Registers R00 thru R03 are to be initialized before executing "GK7" )
• R01 = a
R04 = Gauss7-Integral
R06 to R25: temp
• R02 = b
R05 = Kronrod15-Integral
• R03 = n
Flags: /
Subroutine: A program that returns f(x)
in X-register, assuming x is in X-register upon entry.
| 01 LBL "GK7"
02 .9914553711 03 STO 25 04 .9491079123 05 STO 24 06 .8648644234 07 STO 23 08 .7415311856 09 STO 22 10 .5860872355 11 STO 21 12 .4058451514 13 STO 20 14 .207784955 15 STO 19 16 11.018 17 CLRGX 18 RCL 02 19 RCL 01 20 STO 07 21 - 22 RCL 03 |
23 STO 04
24 ST+ X 25 / 26 STO 06 27 LBL 01 28 ST+ 07 29 25.018 30 STO 08 31 18 32 STO 05 33 LBL 02 34 RCL 07 35 RCL IND 08 36 RCL 06 37 * 38 STO 09 39 - 40 XEQ IND 00 41 STO 10 42 RCL 07 43 RCL 09 44 + |
45 XEQ IND 00
46 RCL 10 47 + 48 ST+ IND 05 49 DSE 05 50 DSE 08 51 GTO 02 52 RCL 07 53 XEQ IND 00 54 ST+ 11 55 RCL 06 56 ST+ 07 57 DSE 04 58 GTO 01 59 RCL 11 60 .2094821411 61 * 62 RCL 12 63 .2044329401 64 * 65 + 66 RCL 13 |
67 .1903505781
68 * 69 + 70 RCL 14 71 .1690047266 72 * 73 + 74 RCL 15 75 .1406532597 76 * 77 + 78 RCL 16 79 .1047900103 80 * 81 + 82 RCL 17 83 15.84984676 84 / 85 + 86 RCL 18 87 43.60087029 88 / |
89 +
90 * 91 STO 05 92 RCL 11 93 .4179591837 94 * 95 RCL 13 96 .3818300505 97 * 98 + 99 RCL 15 100 .2797053915 101 * 102 + 103 RCL 17 104 .1294849662 105 * 106 + 107 RCL 06 108 * 109 STO 04 110 END |
( 344 bytes / SIZE 026 )
| STACK | INPUTS | OUTPUTS |
| Y | / | I15 |
| X | / | I7 |
with I7 = Gauss7-Integral
= R04
and I15 = Kronrod15-Integral
= R05
Example: The same one: compute §04 Sin ( x2 ) dx
LBL "T"
X^2
SIN
RTN
END
"T" ASTO 00
XEQ "RAD"
0 STO 01
4 STO 02
-With n = 2 2 STO 03
XEQ "GK7" >>>> I7 =
0.747111478
---Execution time = 32s---
RDN I15 = 0.747133845
-With n = 4 4 STO 03
R/S
>>>> I7 = 0.747133835
---Execution time = 65s---
RDN I15 = 0.747133845
Variant:
-The same modification may be used:
-Add the following instructions after line 01
| ASTO 00
STO 02 X<>Y STO 01 |
1
STO 03 LBL 03 2 |
ST* 03
XEQ 00 VIEW 04 RND |
X<>Y
RND X#Y? GTO 03 |
LASTX
CLD RTN LBL 00 |
-Place the function name in the alpha "register", a ENTER^ b
XEQ "GK7"
-The loop is executed until 2 rounded approximations are equal
-The accuracy is controlled by the display settings
-With the same example §04 Sin ( x2 ) dx
alpha T alpha
FIX 9
0 ENTER^
4 XEQ "GK7" >>>> I = 0.747133845
( n = R03 = 8 , roundoff-errors are slightly different... )
Reference:
[1] http://en.wikipedia.org/wiki/Gauss–Kronrod_quadrature_formula