ROMBERG QUADRATURE
Estimate definite integral over closed interval with Romberg trapezoidal rule
For SHARP PC-12xx to 14xx
355 bytes BASIC image (PC-1350)
Credit: Numerical Recipes qromb + polint
See also:
https://en.wikipedia.org/wiki/Romberg%27s_method
Precision: 1E-9 (adjustable)
Subdivision steps: up to N=10 (adjustable)
Example:
300 "F1" Y=4/(1+X*X): RETURN
RUN
f=F1
a=0
b=1
3.141592654
Code:
' ROMBERG QUADRATURE
' Estimate definite integral over closed interval with Romberg trapezoidal rule
' For SHARP PC-12xx to 14xx by Robert van Engelen
' Credit:
' Numerical Recipes qromb + 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 qromb(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;
' Ro[0] = (f(a)+f(b))*h/2;
' for (i = 1; i < n; ++i) {
' unsigned long long k = 1UL << i; // k = 2^(i-1)
' unsigned long long s = 1;
' double sum = 0;
' double *Rt;
' h /= 2;
' for (j = 1; j < k; j += 2)
' sum += f(a+j*h);
' Ru[0] = h*sum + Ro[0]/2;
' for (j = 1; j <= i; ++j) {
' s *= 4;
' Ru[j] = (s*Ru[j-1] - Ro[j-1])/(s-1);
' }
' if (i > 2 && 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 (=10)
' I iteration step
' U current row
' O previous row
' J,S,X scratch
' A(27..26+2*N) scratch auto-array
100 "QROMB" E=1E-9,N=10: INPUT "f=F";F$: F$="F"+F$
110 INPUT "a=";A
120 INPUT "b=";B
' init and first trapezoidal step
130 H=B-A,X=A: GOSUB F$: S=Y,X=B: GOSUB F$: O=27,U=O+N,A(O)=H*(S+Y)/2,I=1
' next trapezoidal step
140 H=H/2,S=0
150 FOR J=1 TO 2^I STEP 2: X=A+J*H: GOSUB F$: S=S+Y: NEXT J
' integrate and combine with previous results
160 A(U)=H*S+A(O)/2,S=1
170 FOR J=1 TO I: S=4*S,A(U+J)=(S*A(U+J-1)-A(O+J-1))/(S-1): NEXT J
' loop until convergence
180 IF I>2 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