+- 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) qgaus Gauss-Legendre quadrature (/thread-16548.html)
(PC-12xx~14xx) qgaus Gauss-Legendre quadrature - robve - 03-28-2021 01:47 PM
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