Gauss-Kronrod

# 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

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

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: