03-28-2021, 01:40 PM
ADAPTIVE SIMPSON QUADRATURE
Estimate definite integral over closed interval with adaptive Simpson's rule
For SHARP PC-12xx to 14xx
465 bytes BASIC image (PC-1350)
See also: https://en.wikipedia.org/wiki/Adaptive_S...27s_method
Precision: 1E-9 (adjustable) recommended 1E-5
Recursive subdivision steps: up to 21 (adjustable at line 160)
Example:
300 "F1" Y=4/(1+X*X): RETURN
RUN
f=F1
a=0
b=1
3.141592654
Estimate definite integral over closed interval with adaptive Simpson's rule
For SHARP PC-12xx to 14xx
465 bytes BASIC image (PC-1350)
See also: https://en.wikipedia.org/wiki/Adaptive_S...27s_method
Precision: 1E-9 (adjustable) recommended 1E-5
Recursive subdivision steps: up to 21 (adjustable at line 160)
Example:
300 "F1" Y=4/(1+X*X): RETURN
RUN
f=F1
a=0
b=1
3.141592654
Code:
' ADAPTIVE SIMPSON QUADRATURE
' Estimate definite integral over closed interval with adaptive Simpson's rule
' For SHARP PC-12xx to 14xx by Robert van Engelen
' See also:
' https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method
' Functions to integrate are defined with label "F1", "F2",... should return Y given X
' Algorithm:
' double qasi(double (*f)(double), double a, double b, double eps) {
' double fa = f(a);
' double fm = f((a+b)/2);
' double fb = f(b);
' double v = (fa + 4*fm + fb)*(b-a)/6;
' return qas_(f, a, b, fa, fm, fb, v, eps, 20, 0); // last 0 for tail-recursive version
' }
' double qas_(double (*f)(double), double a, double b, double fa, double fm, double fb, double v, double eps, int n) {
' double h = (b-a)/2;
' double f1 = f(a + h/2);
' double f2 = f(b - h/2);
' double sl = h*(fa + 4*f1 + fm)/6;
' double sr = h*(fm + 4*f2 + fb)/6;
' double s = sl+sr;
' double d = (s-v)/15;
' double m = a+h;
' if (n <= 0 || fabs(d) < eps)
' return s + d; // convergence or max depth exceeded
' return qas_(f, a, m, fa, f1, fm, sl, eps/2, n-1) + qas_(f, m, b, fm, f2, fb, sr, eps/2, n-1);
' }
' Optimized version with tail recursion:
' double qas_(double (*f)(double), double a, double b, double fa, double fm, double fb, double v, double eps, int n, double t) {
' double h = (b-a)/2;
' double f1 = f(a + h/2);
' double f2 = f(b - h/2);
' double sl = h*(fa + 4*f1 + fm)/6;
' double sr = h*(fm + 4*f2 + fb)/6;
' double s = sl+sr;
' double d = (s-v)/15;
' double m = a+h;
' if (n <= 0 || fabs(d) < eps)
' return t + s + d; // convergence or max depth exceeded
' t = qas_(f, a, m, fa, f1, fm, sl, eps/2, n-1, t);
' return qas_(f, m, b, fm, f2, fb, sr, eps/2, n-1, t);
' }
' Non-recursive version with a stack stk[] and stack pointer sp:
' double qas_(double (*f)(double), double a, double b, double fa, double fm, double fb, double v, double eps, int n, double t, int sp, double stk[]) {
' entry:
' double h = (b-a)/2;
' double f1 = f(a + h/2);
' double f2 = f(b - h/2);
' double sl = h*(fa + 4*f1 + fm)/6;
' double sr = h*(fm + 4*f2 + fb)/6;
' double s = sl+sr;
' double d = (s-v)/15;
' double m = a+h;
' // convergence or max depth exceeded?
' if (eps < 1E-15 || fabs(d) < eps) {
' // exit if stack is empty
' if (sp <= 0)
' return t + s + d;
' // pop parameters for the second tail-recursive call
' eps = stk[--sp]; s = stk[--sp]; fb = stk[--sp]; fm = stk[--sp]; fa = stk[--sp]; b = stk[--sp]; a = stk[--sp];
' goto entry;
' }
' eps /= 2;
' // save the parameters for the second "recursive" call on the stack to pop when "returning"
' stk[sp++] = m; stk[sp++] = b; stk[sp++] = fm; stk[sp++] = f2; stk[sp++] = fb; stk[sp++] = sr; stk[sp++] = eps;
' // the first "recursive" call
' b = m; fb = fm; fm = f1; s = sl;
' goto entry;
' }
' A,B range
' D (s-v)/15
' E error 1E-9, changed afterwards
' F$ function label to integrate
' H interval
' I fa=f(a)
' J fm=f((a+b)/2)
' K fb=f(b)
' L sl=h*(fa+4*f1+fm)/6
' M a
' N b
' O stack pointer
' P f1=f(a+h/2)
' Q f2=f(b-h/2)
' R sr=h*(fm+4*f2+fb)/6
' S s=sl+sr
' T total sum
' U partial sum v
' X scratch
' A(27..166) auto-array stack for up to 21 recursive calls with 7 parameters per call (E=2^-20*1E-9<1E-15 at line 160)
' for 11 recursive calls: E=2^-10*1E-9<1E-12 at line 160
100 "QASI" E=1E-9: INPUT "f=F";F$: F$="F"+F$
110 INPUT "a=";A
120 INPUT "b=";B
' init fa=f(a),fb=f(b),fm=f((a+b)/2)
130 X=A: GOSUB F$: I=Y,X=(A+B)/2: GOSUB F$: J=Y,X=B: GOSUB F$: K=Y
140 U=(B-A)*(I+4*J+K)/6,M=A,N=B,T=0,O=27
' recursive call entry
150 H=(N-M)/2,X=M+H/2: GOSUB F$: P=Y,X=N-H/2: GOSUB F$: Q=Y
160 L=H*(I+4*P+J)/6,R=H*(J+4*Q+K)/6,S=L+R,D=(S-U)/15: IF E<1E-15 OR ABS D<E LET T=T+S+D: GOTO 190
' save the parameters for the second "recursive" call on the stack to pop when "returning"
170 E=E/2,A(O)=M+H,A(O+1)=N,A(O+2)=J,A(O+3)=Q,A(O+4)=K,A(O+5)=R,A(O+6)=E,O=O+7
' the first "recursive" call
180 N=M+H,K=J,J=P,U=L: GOTO 150
' exit if stack is empty
190 IF O<=27 LET Y=T: PRINT Y: END
' if stack not empty then pop parameters for the second tail-recursive call
200 O=O-7,M=A(O),N=A(O+1),I=A(O+2),J=A(O+3),K=A(O+4),U=A(O+5),E=A(O+6): GOTO 150