(PC-12xx~14xx) qkron Gauss-Kronrod quadrature

[robve]

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