Simpson revival

[Namir]

Here is the Basic code from Pekis with the improved calculations using Diter's approach

Here is version 1 with the improvement applied only to the final results:

Code:

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 GOTO 160
153 PRINT "Integral value: ";Z
155 PRINT "Improved Integral: "; (16*Z-W)/15
157 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

Here is version 2 where the improvment is applied to each iteration:

Code:

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 GOTO 160
153 PRINT "Integral value: ";Z
155 PRINT "Improved Integral: "; (16*Z-W)/15
157 STOP 'OK with expected tolerance
160 IF W=1E99 THEN W=Z ELSE W=(16*Z-W)/15
165 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

Version 1 gave better results. These results seem counter intuitive, since I expected version 2 to give a better final result.

Any or all improvements and corrections are welcome.

Namir