(PC-12xx~14xx) qasi Adaptive Simpson quadrature

[robve]

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

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

[Albert Chan]

Hi, robve

I translated your BASIC code for HP71B
Example is for ∫(surd(x,-3), x=-2..3) ≈ 0.739024156626

10 DEF FNF(X)=SGN(X)*ABS(X)^(-1/3)
20 A=-2 @ B=3 @ E=.0001 @ DIM Z(210)
30 T2=TIME @ I=FNF(A) @ J=FNF((A+B)/2) @ K=FNF(B)
40 U=(B-A)*(I+4*J+K)/6 @ M=A @ N=B @ T=0 @ O=1

100 H=(N-M)/2 @ P=FNF(M+H/2) @ Q=FNF(N-H/2)
110 L=H*(I+4*P+J)/6 @ R=H*(J+4*Q+K)/6 @ S=L+R @ D=(S-U)/15
120 IF E<1.E-15 OR ABS(D)<E THEN T=T+S+D @ GOTO 160
130 Z(O+1)=N @ N=M+H @ Z(O)=N @ Z(O+2)=J @ Z(O+3)=Q @ Z(O+4)=K @ Z(O+5)=R
140 E=E/2 @ Z(O+6)=E @ K=J @ J=P @ U=L @ O=O+7
150 GOTO 100
160 O=O-7 @ IF O<1 THEN DISP T,TIME-T2 @ END
170 M=Z(O) @ N=Z(O+1) @ I=Z(O+2) @ J=Z(O+3) @ K=Z(O+4) @ U=Z(O+5) @ E=Z(O+6)
180 GOTO 100

>RUN
 .739024205098      10.99

HP71B can do recursive calls, without us handle recursive calling stacks.
I think this version is simpler, which I posted a few days ago.

[robve]

(04-15-2022 01:45 AM)Albert Chan Wrote:  Hi, robve

I translated your BASIC code for HP71B
Example is for ∫(surd(x,-3), x=-2..3) ≈ 0.739024156626

>RUN
 .739024205098      10.99

HP71B can do recursive calls, without us handle recursive calling stacks.
I think this version is simpler, which I posted a few days ago.

Very nice. A bit surprised to see that the recursive version on the HP 71B takes much longer, about 1.7 times the execution time of the stack version:

>RUN
.73902427614 18.76

The quadrature should be comparable, but not necessarily identical as the summation order differs (in-order traversal sum versus a bottom-up sum in the recursive version that is typically more precise).

Applying tail-recursion optimization to the recursive version should speed this up some.

- Rob

[Albert Chan]

(04-16-2022 01:45 PM)robve Wrote:  A bit surprised to see that the recursive version on the HP 71B takes much longer,
about 1.7 times the execution time of the stack version:

If we replace stack version FNF(x) as subroutine, 11s goes up to 14s.
About half of the slowdown is because of the more efficient FNF(x).

SUB SIMPSON() does not see FNF(x), thus required slower SUB F(x,y)

The other half is due to different recurse conditions. (see below)

Quote:The quadrature should be comparable, but not necessarily identical as the summation order differs
(in-order traversal sum versus a bottom-up sum in the recursive version that is typically more precise).

Also, recurse condition of two versions are close, but not identical.

Stack version recurse if E = [1E-15, ABS((S-U)/15)]
SUB SIMPSON recurse if E = [1E-15/16, ABS((S-U)/16)]

For our test example, with undefined value at x=0, SIMPSON version recurse deeper.
With SIMPSON version matching stack version E recurse condition, it produce similar result.
> RUN
 .73902420507      15.44

Compare apples to apples, both using SUB F(x,y), timing difference is now minor.