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

[Namir]

Very clever algorithm! Hats off! The reusing of summations from previous calculations is very clever.

Namir

PS: And thank you for the BASIC code.

[Dieter]

(09-28-2016 08:04 AM)Pekis Wrote:  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.

As it is the case with most clever ideas, you indeed reinvented the wheel. ;-)
Which of course doesn't make the method less clever.

(09-28-2016 08:04 AM)Pekis Wrote:  What do you think of it ?

It's such an obvious improvement that some time ago I posted a HP41 version of this idea. This program even calculates an improved approximation from the two last Simpson results, using the error estimate associated with this method. For instance, if you have the results for 4 and 8 intervals, you can get a new result from these two, having an accuracy comparable to the result with 16 intervals. What about implemeting this idea in your BASIC program?

A long time ago a much better mathematician had a similar idea based on the trapezoid method combined with an extraordinarily clever extrapolation scheme, leading to what is now known as the Romberg method.

Dieter

[Namir]

(09-28-2016 12:25 PM)Dieter Wrote:  

(09-28-2016 08:04 AM)Pekis Wrote:  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.

As it is the case with most clever ideas, you indeed reinvented the wheel. ;-)
Which of course doesn't make the method less clever.

(09-28-2016 08:04 AM)Pekis Wrote:  What do you think of it ?

It's such an obvious improvement that some time ago I posted a HP41 version of this idea. This program even calculates an improved approximation from the two last Simpson results, using the error estimate associated with this method. For instance, if you have the results for 4 and 8 intervals, you can get a new result from these two, having an accuracy comparable to the result with 16 intervals. What about implemeting this idea in your BASIC program?

A long time ago a much better mathematician had a similar idea based on the trapezoid method combined with an extraordinarily clever extrapolation scheme, leading to what is now known as the Romberg method.

Dieter

A few years ago I replaced the Trapezoidal rule use in Romberg with Simpson's method (and also tried using other quadrature algorithms). You can check it here and click on the link A New Face of Romberg Integration . There is a link next to it that allows you to download a related Excel demo file by Graeme Dennes.

Namir

[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

[Dieter]

(09-28-2016 03:50 PM)Namir Wrote:  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.

I am not sure how your code works in detail, but here is an example to compare with.

Edit: I think I now know where your error is. The improved result is calculated from the previous (W) and current (Z) Simpson approximation. But in your second program W is replaced with the last improved approximation, cf. line 160. So the next improved value is calculated from the the previous improved value and the current Simpson approximation. Here seems to be the problem. You must not overwrite W. Leave it as it is and only print/display the improved value.
Or, even better, calculate the previous and current improved approximations WI and ZI and stop the iteration as soon as these agree within the specified tolerance.

Finally, here is my example:
Let f(x) = 1/x and integrate it from 1 to 2.
The exact result, rounded to 12 digits, is ln(2)=0,693147180560.

Code:

   n      Standard Simpson    Error        Improved            Error
 ----------------------------------------------------------------------
   2      0,694444444444      1,3 E-03   
   4      0,693253968254      1,1 E-04     0,693174603175      2,7 E-05
   8      0,693154530655      7,4 E-06     0,693147901481      7,2 E-07
  16      0,693147652819      4,7 E-07     0,693147194297      1,4 E-08
  32      0,693147210290      3,0 E-08     0,693147180788      2,3 E-10
  64      0,693147182421      1,9 E-09     0,693147180564      3,6 E-12
 128      0,693147180676      1,2 E-10     0,693147180560      5,6 E-14

The "improved" column holds the extrapolated results, i.e. [16*Simpson(a, b, n) – Simpson(a, b, n/2)] / 15.
The "Error" columns show the absolute error of the calculated integrals.
You can see that the "improved" values are roughly as accurate as the standard Simpson results with twice the number of intervals, and sometimes even better.
From step to step the error of the standard method improves by a factor approaching 16, while it's finally factor 64 for the improved method.

Here is the VBA code I used to calculate Simpson(1, 2, n):

Code:

Function Simpson(a, b, n)
   h = (b - a) / n
   k = 4
   s = 0
   For i = 1 To n - 1
      s = s + k * f(a + i * h)
      k = 6 - k
   Next
   Simpson = (f(a) + s + f(b)) * h / 3
End Function

Function f(x)
   f = 1 / x
End Function

What results do you get?

Dieter

[Namir]

So my version 1 of my Basic code is the way to go since it does not update the older area value in W.

To be honest, I am a bit disappointed by your algorithm, since the general expectation for improving a calculated value in an iteration should be used in the next iteration and cause enhancing the result and/or reduce the number of iteration. An example is Ostrowski's root solving algorithm where he first uses Newton's method to refine the guess for the root and then (in the same iteration) refines that guess again. The result is that Ostrokswki's method matches the efficiency of Halley's root finding method.

In version 1, the (error of the refined area/area of calculated result) is about 100! It seems that applying the refinement on two good estimates for the integral works well to enhance the final result only.

Namir

[Dieter]

(09-29-2016 04:44 PM)Namir Wrote:  So my version 1 of my Basic code is the way to go since it does not update the older area value in W.

IMHO the way to go is the approach in my HP41 program:

You calculate the Simpson approximation for n=2, 4, 8, ... intervals. After each step the improved value is calculated from the last two approximations. This value is returned to the user. If last two improved (!) approximations agree within the stated tolerance, the iteration exits.

In the given example two "improved" steps yield a similar accuracy as three "standard" Simpson steps, so the final result is obtained faster. If only one single iteration is saved, the total execution time is reduced by 50%.

(09-29-2016 04:44 PM)Namir Wrote:  To be honest, I am a bit disappointed by your algorithm, since the general expectation for improving a calculated value in an iteration should be used in the next iteration and cause enhancing the result and/or reduce the number of iteration.

Just try it the way described above. ;-)

EDIT: Here is a short VBA program that implements the suggested method.

Code:

Function Simpson(a, b, Optional errlimit = 0.000001, Optional nmax = 2000)
   
   fab = f(a) + f(b)
   s2 = 0
   s4 = f((a + b) / 2)
   simp_new = (fab + 4 * s4) * (b - a) / 6
   improved_new = simp_new
   n = 4
Do
   h = (b - a) / n
   s2 = s2 + s4
   s4 = 0
   
   For i = 1 To n - 1 Step 2
      s4 = s4 + f(a + i * h)
   Next
   
   simp_old = simp_new
   simp_new = (fab + 2 * s2 + 4 * s4) * h / 3
   improved_old = improved_new
   improved_new = (16 * simp_new - simp_old) / 15
   n = 2 * n
Loop Until (Abs(improved_new - improved_old) <= errlimit) Or (n > nmax)

Simpson = improved_new

End Function

Function f(x)
   f = 1 / x
End Function

Dieter

[Namir]

A modified version of the BASIC program works if we have two tests:

If Abs(Z - W) <= P Then
....

And

If Abs(Z - (16 * Z - W) / 15) <= P Then
...

I tested f(x)=1/X for integrals in [1,2], [1,10], and [1,100]. In each case the absolute value of the ratio of the error of the calculated integral divided by the error of the improved integral is in the order of 1000!

I suggest omitting the first test altogether!

Here is an updated Basic version that reflects Dieter's suggestion and my above comments:

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
145 Y=(16*Z-W)/15 ' Calculate improved integral
150 IF ABS(Z-Y)<=P THEN PRINT "Integral value: ";Y: 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

Namir

[Dieter]

(09-29-2016 10:27 PM)Namir Wrote:  A modified version of the BASIC program works if we have two tests:

If Abs(Z - W) <= P Then

This compares the old and new Simpson result. As you already said, this test is obsolete.

(09-29-2016 10:27 PM)Namir Wrote:  If Abs(Z - (16 * Z - W) / 15) <= P Then
...

This compares the last Simpson result and the last improved approximation. Does this make sense?

(09-29-2016 10:27 PM)Namir Wrote:  I suggest omitting the first test altogether!

The test should compare the last and previous improved approximations, cf. the posted VBA code. In your program that's the current and previous value of Y.

Dieter

[Dieter]

(09-30-2016 06:00 AM)I Wrote:  The test should compare the last and previous improved approximations, cf. the posted VBA code. In your program that's the current and previous value of Y.

It looks like even this remaining error can be estimated and an even more accurate result can be calculated from the last two improved values. Here is a more structured experimental version.

Code:

Function Simpson(a, b, Optional errlimit = 0.000001, Optional nmax = 3000)
   
   fab = f(a) + f(b)
   s2 = 0
   s4 = 0
   simp_new = 0
   improved_new = 0
   n = 2
Do
   h = (b - a) / n
   s2 = s2 + s4
   s4 = 0
   
   For i = 1 To n - 1 Step 2
      s4 = s4 + f(a + i * h)
   Next
   
   simp_old = simp_new
   improved_old = improved_new
   simp_new = (fab + 2 * s2 + 4 * s4) * h / 3
   
   If n = 2 Then
      improved_new = simp_new
      converged = False
   Else
      improved_new = (16 * simp_new - simp_old) / 15
      converged = Abs(improved_new - improved_old) < 63 * errlimit
   End If
   
   n = 2 * n

Loop Until converged Or n > nmax

Simpson = (64 * improved_new - improved_old) / 63

End Function

Function f(x)
   f = 1 / x
End Function

What do you think?

The following diagram may give an impression for the integral of f(x) = 1/x from a=1 to b=2. Here X is the number of iterations (i.e. n=2^x intervals) and Y is the accuracy, i.e. the number of valid digits. The red line is the classic Simpson method, the blue one represents the "improved" method, and finally the green graph shows the further extrapolated result as shown in the last program.

   

Dieter

[Namir]

(10-02-2016 03:29 PM)Dieter Wrote:  

(09-30-2016 06:00 AM)I Wrote:  The test should compare the last and previous improved approximations, cf. the posted VBA code. In your program that's the current and previous value of Y.

It looks like even this remaining error can be estimated and an even more accurate result can be calculated from the last two improved values. Here is a more structured experimental version.

Code:

Function Simpson(a, b, Optional errlimit = 0.000001, Optional nmax = 3000)
   
   fab = f(a) + f(b)
   s2 = 0
   s4 = 0
   simp_new = 0
   improved_new = 0
   n = 2
Do
   h = (b - a) / n
   s2 = s2 + s4
   s4 = 0
   
   For i = 1 To n - 1 Step 2
      s4 = s4 + f(a + i * h)
   Next
   
   simp_old = simp_new
   improved_old = improved_new
   simp_new = (fab + 2 * s2 + 4 * s4) * h / 3
   
   If n = 2 Then
      improved_new = simp_new
      converged = False
   Else
      improved_new = (16 * simp_new - simp_old) / 15
      converged = Abs(improved_new - improved_old) < 63 * errlimit
   End If
   
   n = 2 * n

Loop Until converged Or n > nmax

Simpson = (64 * improved_new - improved_old) / 63

End Function

Function f(x)
   f = 1 / x
End Function

What do you think?

The following diagram may give an impression for the integral of f(x) = 1/x from a=1 to b=2. Here X is the number of iterations (i.e. n=2^x intervals) and Y is the accuracy, i.e. the number of valid digits. The red line is the classic Simpson method, the blue one represents the "improved" method, and finally the green graph shows the further extrapolated result as shown in the last program.



Dieter

You latest code works in general better than Pekis' version. The ratio of errors generated by your code to that of Pekis range in the orders [1, 10]. Good work!

[Namir]

Thank you Pekis for your code. I had fun translating it into programs for the HP-71B, HP-85B, HP-41C and HP-67. I like your approach because it gives accurate results. Of course, the calculations take more time since the program repeats the integral calculations at ever decreasing increment values.

Namir

[Albert Chan]

(09-28-2016 03:50 PM)Namir Wrote:  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.

Hi, Namir

From the book SICP, I think I know why your version 2 does not give expected better estimate.
The trick is not just apply correction based on 2 iterations, but from all previous iterations.

(09-28-2016 10:00 PM)Dieter Wrote:  Let f(x) = 1/x and integrate it from 1 to 2.
The exact result, rounded to 12 digits, is ln(2)=0,693147180560.

Code:

   n      Standard Simpson    Error        Improved            Error
 ----------------------------------------------------------------------
   2      0,694444444444      1,3 E-03   
   4      0,693253968254      1,1 E-04     0,693174603175      2,7 E-05
   8      0,693154530655      7,4 E-06     0,693147901481      7,2 E-07
  16      0,693147652819      4,7 E-07     0,693147194297      1,4 E-08
  32      0,693147210290      3,0 E-08     0,693147180788      2,3 E-10
  64      0,693147182421      1,9 E-09     0,693147180564      3,6 E-12
 128      0,693147180676      1,2 E-10     0,693147180560      5,6 E-14

The "improved" column holds the extrapolated results, i.e. [16*Simpson(a, b, n) – Simpson(a, b, n/2)] / 15.

The first estimate is correct: T1 + (T1-T0)/(16-1) = 0.693174603175

However, second estimate should apply correction again, from previous estimate.

T1 + (T1-T0)/(64-1) = 0.693147901481 + (-2.67017e-05)/63 = 0.693147477645

Third estimate need 3 corrections (/15, /63, /255), ... This is the revised table.

Code:

   n      Standard Simpson    Error        Improved            Error
 ----------------------------------------------------------------------
   2      0,694444444444      1,3 E-03   
   4      0,693253968254      1,1 E-04     0,693174603175      2,7 E-05
   8      0,693154530655      7,4 E-06     0,693147477645      3,0 E-07
  16      0,693147652819      4,7 E-07     0,693147181917      1,4 E-09
  32      0,693147210290      3,0 E-08     0,693147180562      2,4 E-12
  64      0,693147182421      1,9 E-09     0,693147180560      1,6 E-15

Edit:
After reading the whole thread, it seems the code is not trying to use Romberg's method.
All it wanted was the first approximation of Romberg's extrapolation.
Sorry for the noise.

Anyway, this is what Romberg's method would look like ...

[Albert Chan]

In Python, Romberg's method corrections can be simplified with generator.

Code:

from __future__ import division

def simpson(f, a, b, n=2):          # Simpson generator
    h = (b-a) / 2
    t0 = (f(a)+f(b)) * h
    while 1:
        s = sum(f(a + i*h) for i in xrange(1, n, 2))
        t1 = t0/2 + s*h             # simpson's partial sum
        yield t1 + (t1-t0)/3        # simpson's formula
        t0, n, h = t1, 2*n, h/2

Code:

def _romberg(gen, y, k):            # Romberg generator
    for z in gen:
        yield z + (z-y)/k
        y = z

Code:

def romberg(f, a, b, verbal=False, k=15):
    "Romberg's Method integration"
    gen = simpson(f, a, b, n=2)
    x = gen.next()
    while 1:
        y = gen.next()
        x = y + (y-x) / k
        if verbal: print x
        if x == y: return x         # stop if converged
        gen = _romberg(gen, y, k)   # Romberg corrections
        k = 4 * k + 3

Redo previous post example:

>>> romberg(lambda x: 1/x, 1, 2, verbal=True)
0.693174603175
0.693147477645
0.693147181917
0.693147180562
0.69314718056
0.69314718056
0.69314718055994506