# HP Forums

Full Version: (PC-12xx~14xx) qkron Gauss-Kronrod quadrature
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
GAUSS-KRONROD QUADRATURE
Estimate definite integral over open interval with Gauss-Kronrod (G10,K21) at 21 points
For SHARP PC-12xx to 14xx

See also: https://en.wikipedia.org/wiki/Gauss–Kron...re_formula

Precision: 1E-9 (returns estimated actual error in E)
Points: 21

788 bytes BASIC image (PC-1350)

Example:
300 "F1" Y=4/(1+X*X): RETURN
RUN
f=F1
a=0
b=1
3.141592654
310 "F2" Y=ATN(1/X): RETURN
RADIAN
RUN
f=F2
a=0
b=1
1.131971754

Code:
```' GAUSS-KRONROD QUADRATURE ' Estimate definite integral over open interval with Gauss-Kronrod (G10,K21) at 21 points ' For SHARP PC-12xx to 14xx by Robert van Engelen ' See also: '   https://en.wikipedia.org/wiki/Gauss–Kronrod_quadrature_formula ' Functions to integrate are defined with label "F1", "F2",... should return Y given X ' Algorithm: '   double qkron(double (*f)(double), double a, double b) { '     // abscissas and weights pre-calculated from Legendre Stieltjes polynomials '     static const double abscissas[11] = { '       0.00000000000000000e+00, '       1.48874338981631211e-01, '       2.94392862701460198e-01, '       4.33395394129247191e-01, '       5.62757134668604683e-01, '       6.79409568299024406e-01, '       7.80817726586416897e-01, '       8.65063366688984511e-01, '       9.30157491355708226e-01, '       9.73906528517171720e-01, '       9.95657163025808081e-01, '     }; '     static const double weights[11] = { '       1.49445554002916906e-01, '       1.47739104901338491e-01, '       1.42775938577060081e-01, '       1.34709217311473326e-01, '       1.23491976262065851e-01, '       1.09387158802297642e-01, '       9.31254545836976055e-02, '       7.50396748109199528e-02, '       5.47558965743519960e-02, '       3.25581623079647275e-02, '       1.16946388673718743e-02, '     }; '     static const double gauss_weights[5] = { '       2.95524224714752870e-01, '       2.69266719309996355e-01, '       2.19086362515982044e-01, '       1.49451349150580593e-01, '       6.66713443086881376e-02, '     }; '     double c = (a+b)/2; '     double d = (b-a)/2; '     double r = 0; // kronrod result '     double s = 0; // gauss result '     double fp, fm; '     double e; '     int i; '     r = f(c) * weights[0]; '     for (i = 1; i < 11; i += 2) { '       fp = f(c + d*abscissas[i]); '       fm = f(c - d*abscissas[i]); '       r += (fp + fm) * weights[i]; '       s += (fp + fm) * gauss_weights[i/2]; '     } '     for (i = 2; i < 11; i += 2) '     { '       fp = f(c + d*abscissas[i]); '       fm = f(c - d*abscissas[i]); '       r += (fp + fm) * weights[i]; '     } '     return d*r; // estimated error is fabs(r-s) '   } ' VARIABLES '  A,B       range '  F\$        function label to integrate '  Y         result with error E '  E         estimated error '  C         (a+b)/2 '  D         (b-a)/2 '  H         step size '  R         Kronrod quadrature sum '  S         Gauss quadrature sum '  I,O,U,X   scratch '  A(27..52) scratch abscissas and weights 100 "QKRON" INPUT "f=F";F\$: F\$="F"+F\$ 110 INPUT "a=";A 120 INPUT "b=";B ' init Gauss-Kronrod abscissas (11), the first is unused and not initialized 130 A(27)=.1488743390,A(28)=.2943928627,A(29)=.4333953941 140 A(30)=.5627571347,A(31)=.6794095683,A(32)=.7808177266,A(33)=.8650633667 150 A(34)=.9301574914,A(35)=.9739065285,A(36)=.9956571630 ' init Gauss-Kronrod weights (11) and init quadrature sums 160 A(37)=.1494455540,A(38)=.1477391049,A(39)=.1427759386,A(40)=.1347092173 170 A(41)=.1234919763,A(42)=.1093871588,A(43)=.09312545458,A(44)=.07503967481 180 A(45)=.05475589657,A(46)=.03255816231,A(47)=.01169463887 ' init Gauss weights (5) 190 A(48)=.2955242247,A(49)=.2692667193,A(50)=.2190863625,A(51)=.1494513492 200 A(52)=.06667134431,C=(A+B)/2,D=(B-A)/2,X=C: GOSUB F\$: R=A(37)*Y,S=0 ' loop over Gauss and Kronrod points 210 FOR I=1 TO 5 220 H=D*A(25+I+I),X=C+H: GOSUB F\$: O=Y,X=C-H: GOSUB F\$: R=R+A(36+I+I)*(O+Y),S=S+A(47+I)*(O+Y) 230 H=D*A(26+I+I),X=C+H: GOSUB F\$: O=Y,X=C-H: GOSUB F\$: R=R+A(37+I+I)*(O+Y) 240 NEXT I ' done 250 Y=D*R,E=ABS(R-S): IF E>1E-9 PRINT Y,E: END 260 PRINT Y: END```
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :