(PC-12xx~14xx) qgaus Gauss-Legendre quadrature

[robve]

GAUSS-LEGENDRE 10 POINT QUADRATURE
Estimate definite integral over open interval with Gauss Legendre
For SHARP PC-12xx to 14xx

349 bytes BASIC image (PC-1350)

Credit: Numerical Recipes qgaus
See also: https://en.wikipedia.org/wiki/Gaussian_quadrature

Precision: 1E-3 to 1E-10
Points: 10 (fixed)

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
RUN
f=F2
a=0
b=1
1.131971754

Code:

' GAUSS-LEGENDRE 10 POINT QUADRATURE
' Estimate definite integral over open interval with Gauss Legendre
' For SHARP PC-12xx to 14xx by Robert van Engelen
' Credit:
'   Numerical Recipes qgaus
' See also:
'   https://en.wikipedia.org/wiki/Gaussian_quadrature

' No error tolerance and no error estimation, error is as low as 1E-9 for smooth functions

' Functions to integrate are defined with label "F1", "F2",... should return Y given X

' Algorithm:
'   double qgaus(double (*f)(double), double a, double b) {
'     static const double w[5] = { .2955242247, .2692667193, .2190863625, .1494513491, .0666713443 };
'     static const double x[5] = { .1488743389, .4333953941, .6794095682, .8650633666, .9739065285 };
'     double s = 0;
'     double xm = (b+a)/2;
'     double xr = (b-a)/2;
'     int i;
'     for (i = 0; i < 5; i++) {
'       double h = xr*x[i];
'       s += w[i]*(f(xm+h)+f(xm-h));
'     }
'     return xr*s;
'   }

' VARIABLES
'  A,B       range
'  F$        function label to integrate
'  Y         result with unknown error
'  C         (a+b)/2
'  D         (b-a)/2
'  H         step size
'  I,U,X     scratch
'  A(27..36) scratch weights and abscissas

100 "QGAUS" INPUT "f=F";F$: F$="F"+F$
110 INPUT "a=";A
120 INPUT "b=";B
' init sum=0, xm=(a+b)/2, xr=(b-a)/2, weights and abscissas
130 U=0,C=(A+B)/2,D=(B-A)/2,A(27)=.2955242247,A(28)=.1488743389
140 A(29)=.2692667193,A(30)=.4333953941,A(31)=.2190863625,A(32)=.6794095682
150 A(33)=.1494513491,A(34)=.8650633666,A(35)=.0666713443,A(36)=.9739065285
' integrate over 10 points
160 FOR I=27 TO 35 STEP 2: H=D*A(I+1),X=C+H: GOSUB F$: U=U+A(I)*Y,X=C-H: GOSUB F$: U=U+A(I)*Y: NEXT I
170 Y=D*U: PRINT Y: END