(PC-12xx~14xx) qromo Romberg Quadrature open intervals

[robve]

ROMBERG QUADRATURE FOR OPEN INTERVALS
Estimate definite integral over open interval with Romberg midpoint rule
For SHARP PC-12xx to 14xx

369 bytes BASIC image (PC-1350)

Credit: Numerical Recipes qromo + midpnt + polint
See also: https://en.wikipedia.org/wiki/Romberg%27s_method

Precision: 1E-9 (adjustable)
Subdivision steps: up to N=7 (adjustable)

Example:
300 "F1" Y=4/(1+X*X): RETURN
RUN
f=F1
a=0
b=1
3.141592656
310 "F2" Y=ATN(1/X): RETURN
RUN
f=F2
a=0
b=1
1.131971753

Code:

' ROMBERG QUADRATURE FOR OPEN INTERVALS
' Estimate definite integral over open interval with Romberg midpoint rule
' For SHARP PC-12xx to 14xx by Robert van Engelen
' Credit:
'   Numerical Recipes qromo + midpnt + polint
' See also:
'   https://en.wikipedia.org/wiki/Romberg%27s_method

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

' Algorithm:
'   double qromo(double (*f)(double), double a, double b, int n, double eps) {
'     double R1[n], R2[n];
'     double *Ro = &R1[0], *Ru = &R2[0];
'     double h = b-a;
'     int i, j;
'     unsigned long long k = 1;
'     Ro[0] = f((a+b)/2)*h;
'     for (i = 1; i < n; ++i) {
'       unsigned long long s = 1;
'       double sum = 0;
'       double *Rt;
'       k *= 3;
'       h /= 3;
'       for (j = 1; j < k; j += 3)
'         sum += f(a+(j-1)*h+h/2) + f(a+(j+1)*h+h/2);
'       Ru[0] = h*sum + Ro[0]/3;
'       for (j = 1; j <= i; ++j) {
'         s *= 9;
'         Ru[j] = (s*Ru[j-1] - Ro[j-1])/(s-1);
'       }
'       if (i > 1 && fabs(Ro[i-1]-Ru[i]) <= eps*fabs(Ru[i])+eps)
'         return Ru[i];
'       Rt = Ro;
'       Ro = Ru;
'       Ru = Rt;
'     }
'     return Ro[n-1]; // no convergence, return best result, error is fabs((Ru[n-2]-Ro[n-1])/Ro[n-1])
'   }

' VARIABLES
'  A,B           range
'  F$            function label to integrate
'  Y             result
'  E             relative error: integral = Y with precision E (attempts E = 1E-10)
'  H             step size
'  N             max number of Romberg steps (=7)
'  I             iteration step
'  U             current row
'  O             previous row
'  J,S,X         scratch
'  A(27..26+2*N) scratch auto-array

100 "QROMO" E=1E-9,N=7: INPUT "f=F";F$: F$="F"+F$
110 INPUT "a=";A
120 INPUT "b=";B
' init and first midpoint step
130 H=B-A,X=A+H/2: GOSUB F$: O=27,U=O+N,A(O)=H*Y,I=1
' next midpoint step
140 H=H/3,S=0
150 FOR J=1 TO 3^I STEP 3: X=A+(J-.5)*H: GOSUB F$: S=S+Y,X=A+(J+1.5)*H: GOSUB F$: S=S+Y: NEXT J
' integrate and combine with previous results
160 A(U)=H*S+A(O)/3,S=1
170 FOR J=1 TO I: S=9*S,A(U+J)=(S*A(U+J-1)-A(O+J-1))/(S-1): NEXT J
' loop until convergence
180 IF I>1 LET Y=A(U+I): IF ABS(Y-A(O+I-1))<=E*ABS Y+E PRINT Y: END
190 J=O,O=U,U=J,I=I+1: IF I<N GOTO 140
' no convergence, output result with error estimate
200 E=ABS(Y-A(U+N-2))/(ABS Y+E): PRINT Y,E: END