Post Reply 
(PC-12xx~14xx) qasi Adaptive Simpson quadrature
03-28-2021, 01:40 PM
Post: #1
(PC-12xx~14xx) qasi Adaptive Simpson quadrature
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

HP 71B,Prime G2;Ti VOY200,Nspire CXII CAS;Casio fx-CG50,fx-115ES+2;Sharp PC-G850VS,E500S,1475,1450,1360,1350,2500,1262,1500A
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 




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