Adaptive Simpson and Romberg integration methods

[toml_12953]

Very nice routine! Here's how to do it in a modern BASIC

Code:

DECLARE EXTERNAL FUNCTION qasi

PRINT qasi(0, 1, 1e-9)

END 

!
! Define the function to be integrated here
!
EXTERNAL FUNCTION f(x)

LET f = x^3-2*x^2+x

END FUNCTION


EXTERNAL FUNCTION qasi(a, b, ep)

DECLARE EXTERNAL FUNCTION qas_, f

LET fa = f(a)
LET fm = f((a+b)/2)
LET fb = f(b)
LET v  = (fa + 4*fm + fb)*(b-a)/6
LET qasi = qas_(a, b, fa, fm, fb, v, ep, 20)

END FUNCTION



EXTERNAL FUNCTION qas_(a, b, fa, fm, fb, v, ep, n)

DECLARE EXTERNAL FUNCTION f

LET h  = (b-a)/2
LET f1 = f(a + h/2)
LET f2 = f(b - h/2)
LET sl = h*(fa + 4*f1 + fm)/6
LET sr = h*(fm + 4*f2 + fb)/6
LET s  = sl+sr
LET d  = (s-v)/15
LET m  = a+h
IF (n <= 0 OR ABS(d) < ep) THEN
   LET qas_ = t + s + d ! convergence or max depth exceeded
   EXIT FUNCTION
END IF 
LET t = qas_(a, m, fa, f1, fm, sl, ep/2, n-1)
LET qas_ = qas_(m, b, fm, f2, fb, sr, ep/2, n-1)

END FUNCTION