(PC-12xx~14xx) qthsh Tanh-Sinh quadrature
TANH-SINH QUADRATURE and COMBINED TANH-SINH, EXP-SINH, SINH-SINH QUADRATURES
Estimate definite integral over open interval with improved and optimized Tanh-Sinh quadrature
Credit: Michalski & Mosig Tanh-Sinh rule, improved in this thread and in the accompanying technical paper
For details, see the technical paper: https://www.genivia.com/files/qthsh.pdf
See also: https://newtonexcelbach.com/2020/10/29/n...ure-v-5-0/
For SHARP PC-12xx to 14xx
373 bytes BASIC image (PC-1350) "qthsh" optimized Tanh-Sinh
936 bytes BASIC image (PC-1350) "quad" optimized combined Tanh-Sinh, Exp-Sinh, Sinh-Sinh
Precision: 1E-9 (adjustable)
Levels: N=6 (adjustable)
Example:
300 "F1" Y=4/(1+X*X): RETURN
RUN
f=F1
a=0
b=1
3.141592653
310 "F2" Y=ATN(1/X): RETURN
RADIAN
RUN
f=F2
a=0
b=1
1.131971754
320 "F3" Y=COS X*LN X: RETURN
RADIAN
RUN
f=F3
a=0
b=1
-9.460830706E-01
330 "F4" Y=X^-.8: RETURN
RUN
f=F4
a=0
b=1
5
Code:
' TANH-SINH QUADRATURE
' Estimate definite integral over open interval with improved and optimized Tanh-Sinh quadrature
' For SHARP PC-12xx to 14xx by Robert van Engelen
' Credit:
' Michalski & Mosig Tanh-Sinh rule
' See also:
' https://www.genivia.com/files/qthsh.pdf
' https://newtonexcelbach.com/2020/10/29/numerical-integration-with-tanh-sinh-quadrature-v-5-0/
' Functions to integrate are defined with label "F1", "F2",... should return Y given X
' Algorithm:
' double qthsh(double (*f)(double), double a, double b, double n, double eps) {
' double c = (a+b)/2; // center (mean)
' double d = (b-a)/2; // half distance
' double s = f(c);
' double e, v, h = 2;
' int k = 0;
' if (n <= 0) // use default levels n=6
' n = 6; // 6 is optimal, 7 just as good taking longer
' if (eps <= 0) // use default eps=1E=9
' eps = 1E-9;
' do {
' double p = 0, q, fp = 0, fm = 0, t, eh;
' h /= 2;
' eh = exp(h);
' t = eh;
' if (k > 0)
' eh *= eh;
' do {
' double u = exp(1/t-t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
' double r = 2*u/(1+u); // = 1 - tanh(sinh(j*h))
' double w = (t+1/t)*r/(1+u); // = cosh(j*h)/cosh(sinh(j*h))^2
' double x = d*r;
' if (a+x > a) { // if too close to a then reuse previous fp
' double y = f(a+x);
' if (isfinite(y))
' fp = y; // if f(x) is finite, add to local sum
' }
' if (b-x < b) { // if too close to b then reuse previous fm
' double y = f(b-x);
' if (isfinite(y))
' fm = y; // if f(x) is finite, add to local sum
' }
' q = w*(fp+fm);
' p += q;
' t *= eh;
' } while (fabs(q) > eps*fabs(p));
' v = s-p;
' s += p;
' ++k;
' } while (fabs(v) > 10*eps*fabs(s) && k <= n);
' e = fabs(v)/(fabs(s)+eps);
' return d*s*h; // result with estimated relative error e
' }
' VARIABLES
' A,B range
' F$ function label to integrate
' Y result with error E
' E estimated relative error of the result
' N levels (up to 6 or 7)
' C (a+b)/2 center
' D (b-a)/2 half distance
' H step size h=2^-k
' K level counter
' P,Q,S quadrature sums
' T exp(j*h)
' L fp=f(a+x)
' M fm=f(b-x)
' J,R,U,X scratch
100 "QTHSH" E=1E-9,N=6: INPUT "f=F";F$: F$="F"+F$
110 INPUT "a=";A
120 INPUT "b=";B
' init
130 C=(A+B)/2,D=(B-A)/2,X=C: GOSUB F$: S=Y,H=1,K=0
' outer loop
140 J=1,P=0,L=0,M=0
' inner loop
150 T=EXP(J*H),U=EXP(1/T-T),R=2*U/(1+U)
160 X=A+D*R: IF X>A GOSUB F$: L=Y
170 X=B-D*R: IF X<B GOSUB F$: M=Y
180 Q=(T+1/T)*R/(1+U)*(L+M),P=P+Q,J=J+1+(K>0)
190 IF ABS Q>E*ABS P GOTO 150
' exit inner loop
200 X=S-P,S=S+P,K=K+1
210 IF ABS X>10*E*ABS S IF K<=N LET H=H/2: GOTO 140
' exit outer loop, output result (and relative error estimate if >E)
220 Y=D*S*H,U=ABS X/(ABS S+E)
230 IF U>E LET E=U: PRINT Y,E: END
240 E=U: PRINT Y: END
' For machines with ON ERROR GOTO such as PC-1475, update the following:
' 100 "QTHSH" ON ERROR GOTO 290: E=1E-9,N=6: INPUT "f=F";F$: F$="F"+F$
' 160 X=A+D*R,Y=L: IF X>A GOSUB F$: L=Y
' 170 X=B-D*R,Y=M: IF X<B GOSUB F$: M=Y
' 290 IF ERN=2 RETURN
Code:
' COMBINED TANH-SINH, EXP-SINH, SINH-SINH QUADRATURE
' Estimate integral over open interval with improved Tanh-Sinh, Exp-Sinh and Sinh-Sinh quadratures
' For SHARP PC-12xx to 14xx by Robert van Engelen
' Credit:
' Michalski & Mosig Tanh-Sinh rule
' See also:
' https://www.genivia.com/files/qthsh.pdf
' https://newtonexcelbach.com/2020/10/29/numerical-integration-with-tanh-sinh-quadrature-v-5-0/
' Functions to integrate are defined with label "F1", "F2",... should return Y given X
' Algorithm:
' double quad(double (*f)(double), double a, double b, double n, double eps) {
' double c = 0, d = 1, s, sign = 1, e, v, h = 2;
' int k = 0, mode = 0; // Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
' if (b < a) { // swap bounds
' v = b;
' b = a;
' a = v;
' sign = -1;
' }
' if (isfinite(a) && isfinite(b)) {
' c = (a+b)/2;
' d = (b-a)/2;
' v = c;
' }
' else if (isfinite(a)) {
' mode = 1; // Exp-Sinh
' c = a;
' v = a+d;
' }
' else if (isfinite(b)) {
' mode = 1; // Exp-Sinh
' c = b;
' v = b-d;
' d = -d;
' sign = -sign;
' }
' else {
' mode = 2; // Sinh-Sinh
' v = 0;
' }
' s = f(v);
' do {
' double p = 0, q, fp = 0, fm = 0, t, eh;
' h /= 2;
' t = eh = exp(h);
' if (k > 0)
' eh *= eh;
' if (mode == 0) { // Tanh-Sinh
' do {
' double u = exp(1/t-t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
' double r = 2*u/(1+u); // = 1 - tanh(sinh(j*h))
' double w = (t+1/t)*r/(1+u); // = cosh(j*h)/cosh(sinh(j*h))^2
' double x = d*r;
' if (a+x > a) { // if too close to a then reuse previous fp
' double y = f(a+x);
' if (isfinite(y))
' fp = y; // if f(x) is finite, add to local sum
' }
' if (b-x < b) { // if too close to b then reuse previous fm
' double y = f(b-x);
' if (isfinite(y))
' fm = y; // if f(x) is finite, add to local sum
' }
' q = w*(fp+fm);
' p += q;
' t *= eh;
' } while (fabs(q) > eps*fabs(p));
' }
' else {
' t /= 2;
' do {
' double r = exp(t-.25/t); // = exp(sinh(j*h))
' double x, y, w = r;
' q = 0;
' if (mode == 1) { // Exp-Sinh
' x = c + d/r;
' if (x == c) // if x hit the finite endpoint then break
' break;
' y = f(x);
' if (isfinite(y)) // if f(x) is finite, add to local sum
' q += y/w;
' }
' else { // Sinh-Sinh
' r = (r-1/r)/2; // = sinh(sinh(j*h))
' w = (w+1/w)/2; // = cosh(sinh(j*h))
' x = c - d*r;
' y = f(x);
' if (isfinite(y)) // if f(x) is finite, add to local sum
' q += y*w;
' }
' x = c + d*r;
' y = f(x);
' if (isfinite(y)) // if f(x) is finite, add to local sum
' q += y*w;
' q *= t+.25/t; // q *= cosh(j*h)
' p += q;
' t *= eh;
' } while (fabs(q) > eps*fabs(p));
' }
' v = s-p;
' s += p;
' ++k;
' } while (fabs(v) > 10*eps*fabs(s) && k <= n);
' e = fabs(v)/(fabs(s)+eps);
' return sign*d*s*h; // result with estimated relative error e
' }
' VARIABLES
' A,B range, -9E99 and 9E99 are -inf and +inf
' F$ function label to integrate
' Y result with error E
' E estimated relative error of the result
' N levels (up to 6 or 7)
' C center
' D half distance (Tanh-Sinh) or D=1 (Exp-Sinh)
' G sign
' H step size h=2^-k
' K level counter
' P,Q,S quadrature sums
' T exp(j*h)
' L f(x) point
' M f(x) point
' J,R,U,X scratch
100 "QUAD" E=1E-9,N=6: INPUT "f=F";F$: F$="F"+F$
110 INPUT "a=";A
120 INPUT "b=";B
' init and swap bounds if b<a
130 D=1,G=1,H=1,K=0: IF B<A LET X=A,A=B,B=X,G=-1
140 IF ABS A<9E99 IF ABS B<9E99 GOTO 180
150 IF ABS A<9E99 LET C=A,X=A+D: GOTO 300
160 IF ABS B<9E99 LET C=B,X=B-D,D=-D,G=-G: GOTO 300
170 GOTO 400
' Tanh-Sinh
180 C=(A+B)/2,D=(B-A)/2,X=C: GOSUB F$: S=Y
' outer loop
190 J=1,P=0,L=0,M=0
' inner loop
200 T=EXP(J*H),U=EXP(1/T-T),R=2*U/(1+U)
210 X=A+D*R: IF X>A GOSUB F$: L=Y
220 X=B-D*R: IF X<B GOSUB F$: M=Y
230 Q=(T+1/T)*R/(1+U)*(L+M),P=P+Q,J=J+1+(K>0)
240 IF ABS Q>E*ABS P GOTO 200
' exit inner loop
250 X=S-P,S=S+P,K=K+1
260 IF ABS X>10*E*ABS S IF K<=N LET H=H/2: GOTO 190
' exit, output result (and relative error estimate if >E)
270 Y=G*D*S*H,U=ABS X/(ABS S+E)
280 IF U>E LET E=U: PRINT Y,E: END
290 E=U: PRINT Y: END
' Exp-Sinh
300 GOSUB F$: S=Y
' outer loop
310 J=1,P=0
' inner loop
320 T=EXP(J*H)/2,U=EXP(T-.25/T)
330 X=C+D/U: IF X=C GOTO 370
340 GOSUB F$: L=Y/U,X=C+D*U: GOSUB F$: M=Y*U
350 Q=(T+.25/T)*(L+M),P=P+Q,J=J+1+(K>0)
360 IF ABS Q>E*ABS P GOTO 320
' exit inner loop
370 X=S-P,S=S+P,K=K+1
380 IF ABS X>10*E*ABS S IF K<=N LET H=H/2: GOTO 310
390 GOTO 270
' Sinh-Sinh
400 X=0: GOSUB F$: S=Y
' outer loop
410 J=1,P=0
' inner loop
420 T=EXP(J*H)/2,U=EXP(T-.25/T)/2,R=U-.25/U
430 X=-R: GOSUB F$: L=Y,X=R: GOSUB F$: M=Y
440 Q=(T+.25/T)*(U+.25/U)*(L+M),P=P+Q,J=J+1+(K>0)
450 IF ABS Q>E*ABS P GOTO 420
' exit inner loop
460 X=S-P,S=S+P,K=K+1
470 IF ABS X>10*E*ABS S IF K<=N LET H=H/2: GOTO 410
480 GOTO 270
' For machines with ON ERROR GOTO such as PC-1475, update the following:
' 100 "QUAD" ON ERROR GOTO 490: E=1E-9,N=6: INPUT "f=F";F$: F$="F"+F$
' 210 X=A+D*R,Y=L: IF X>A GOSUB F$: L=Y
' 220 X=B-D*R,Y=M: IF X<B GOSUB F$: M=Y
' 490 IF ERN=2 RETURN
Edit: major update with the latest version.
"I count on old friends to remain rational"
|