+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: Not HP Calculators (/forum-7.html)
+--- Forum: Not remotely HP Calculators (/forum-9.html)
+--- Thread: (PC-12xx~14xx) qkron Gauss-Kronrod quadrature (/thread-16550.html)
(PC-12xx~14xx) qkron Gauss-Kronrod quadrature - robve - 03-28-2021 07:06 PM
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–Kronrod_quadrature_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