(PC-12xx~14xx) qromb Romberg Quadrature

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