Simpson revival
09-28-2016, 03:50 PM (This post was last modified: 09-28-2016 05:17 PM by Namir.)
Post: #5 Namir Senior Member Posts: 822 Joined: Dec 2013
RE: Simpson revival
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
 « Next Oldest | Next Newest »

 Messages In This Thread Simpson revival - Pekis - 09-28-2016, 08:04 AM RE: Simpson revival - Namir - 09-28-2016, 10:40 AM RE: Simpson revival - Dieter - 09-28-2016, 12:25 PM RE: Simpson revival - Namir - 09-28-2016, 02:00 PM RE: Simpson revival - Namir - 09-28-2016 03:50 PM RE: Simpson revival - Dieter - 09-28-2016, 10:00 PM RE: Simpson revival - Albert Chan - 07-31-2018, 02:57 PM RE: Simpson revival - Namir - 09-29-2016, 04:44 PM RE: Simpson revival - Dieter - 09-29-2016, 06:23 PM RE: Simpson revival - Namir - 09-29-2016, 10:27 PM RE: Simpson revival - Dieter - 09-30-2016, 06:00 AM RE: Simpson revival - Dieter - 10-02-2016, 03:29 PM RE: Simpson revival - Namir - 10-02-2016, 04:48 PM RE: Simpson revival - Namir - 10-09-2016, 03:14 AM RE: Simpson revival - Albert Chan - 08-04-2018, 05:29 PM

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