(PC-12xx~14xx) qromo Romberg Quadrature open intervals - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: Not HP Calculators (/forum-7.html) +--- Forum: Not remotely HP Calculators (/forum-9.html) +--- Thread: (PC-12xx~14xx) qromo Romberg Quadrature open intervals (/thread-16546.html) (PC-12xx~14xx) qromo Romberg Quadrature open intervals - robve - 03-28-2021 01:33 PM 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