Simpson revival

[Pekis]

Hello,

The problem with Simpson's rule for calculating a definite integral is that you have to specify the number of intervals. I surely reinvented the wheel but I found a way to deal with a given tolerance (e.g. 1E-5), WITHOUT wasting any calculation in the process. It's very useful when you have to call integration in the middle of a root solving problem.

Beginning with 2 intervals, it doubles the number of intervals at each step and recycles previous results.
That lead to this compact program (here in BASIC (for Namir ):

What do you think of it ?

Variables
A,B: End points
f(X): Function to integrate
P: Expected Tolerance
Z: Integral value
W: Z from previous step
N: Number of intervals
L: N from previous step
H: Size of one interval ( (B-A)/N )
S: Sum of 4*f(x) for new (odd) points for this step
R: S from previous step
U: Sum of 2*f(x) for even points, known from previous steps
T: U from previous step, initialized with f(A)+f(B)

Program
10 INPUT A
20 INPUT B
30 P=1E-5 'Expected tolerance
40 R=0: L=1: N=1: W=1E99 'Initializing
50 X=A: GOSUB 1000: T=Y 'Call Y=f(X) subroutine with X=A
60 X=B: GOSUB 1000: T=T+Y 'T=f(A)+f(B) now
70 N=N*2: H=(B-A)/N 'New step: Double the number of intervals
80 X=A+H: S=0 'Start with 1st new (odd) point
90 FOR I=1 TO L 'Calculate Sum of f(x) for new (odd) points
100 GOSUB 1000: S=S+Y
110 X=X+2*H 'Next new (odd) point
120 NEXT I
130 S=S*4: U=T+R/2 'Odd points have a coefficient of 4; Even points come from old odd points but with a coefficient of 2 instead of 4
140 Z=H*(S+U)/3 'New Integral value at this step
150 IF ABS(Z-W)<=P THEN PRINT "Integral value: ";Z: STOP 'OK with expected tolerance
160 W=Z: T=U: R=S: L=N 'Prepare new step
170 GOTO 70 'New step, please

1000 Y=f(X) 'Calculate f(X) subroutine
1010 RETURN