(PC-12xx~14xx) qromb Romberg Quadrature
03-28-2021, 01:24 PM
Post: #1
 robve Senior Member Posts: 398 Joined: Sep 2020
(PC-12xx~14xx) qromb Romberg Quadrature
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

"I count on old friends" -- HP 71B,Prime|Ti VOY200,Nspire CXII CAS|Casio fx-CG50...|Sharp PC-G850,E500,2500,1500,14xx,13xx,12xx...
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)