The most compact algorithm for Simpson's rule??

[Namir]

Here is a (corrected, based on feedback from site members) pseudo code for what I hope to be the most compact way to implement Simpson's rule:

Code:


Given function f(x), integration interval [A, B] and N divisions, where N is even.
To calculate the integral of f(x) for X=A to x=B using Simpson's Rule

h=(B-A)/N
Sum=f(A)-f(B)
CHS=1
A=A+h
Do
  Sum = Sum + (3+CHS)*f(A)
  CHS=-CHS
  A=A+h
  N=N-1
Loop Until N=0
Area=h/3*Sum

Many implementations use two summations--one for odd terms and one for even terms. In addition these implementations use two loops with slightly different ranges. The above pseudo-coode implements a simple version that uses one loop and a few simple computational tricks.

Namir

[Thomas Klemm]

(12-12-2015 04:24 PM)Namir Wrote:  

Code:

Loop Until B<A

This should be something like \(A\approx B\) or else you may end up with \(A=B+h\).
I'd probably use \(N\) instead to control the loop.

[Namir]

Using N should also work.

[Namir]

(12-12-2015 05:44 PM)Thomas Klemm Wrote:  

(12-12-2015 04:24 PM)Namir Wrote:  

Code:

Loop Until B<A

This should be something like \(A\approx B\) or else you may end up with \(A=B+h\).
I'd probably use \(N\) instead to control the loop.

Yes, A ends up as B+h causing the loop to stop iterating.

[Dieter]

(12-12-2015 04:24 PM)Namir Wrote:  Here is a pseudo code for what I hope to be the most compact way to implement Simpson's rule:

I'm not sure if this counts as "more compact", but the alternating factors k=4, 2, 4, 2,... can be generated by subtraction from the sum of both, i.e. k := 6 – k.

I would also use a simple loop counter instead of continuously adding h (which may lead to roundoff errors).

And finally, the initial value of "sum" of course has to be set to the sum of f(a) and f(b), not the difference.

Edit: I now see this is not an error since you calculate f(b) twice, once before the loop and once within it.
The latter value is added twice, so all this sums up to –f(b)+2f(b) = +f(b), which is correct. Only the exit condition remains wrong.

Also n has to be even, not odd. The number of nodes between a and b is odd, since it equals n–1.

Code:

h = (b-a)/n
sum = f(a) + f(b)
k = 4

for i = 1 to n-1 do
  sum = sum + k * f(a + i*h)
  k = 6 - k
next

result = sum * h / 3

(12-12-2015 04:24 PM)Namir Wrote:  Many implementations use two summations--one for odd terms and one for even terms.

Right – this has a big advantage: If both sums are calculated separately, adding another iteration with 2·n is very easy since half of the function calls have already been evaluated. So you can calculate approximations for n=2, 4, 8, 16, ... until the results agree to the desired accuracy. This way you also get an estimate for the approximation error.

Dieter

[Dieter]

(12-12-2015 08:29 PM)Namir Wrote:  Yes, A ends up as B+h causing the loop to stop iterating.

Since the test if A>B is done *after* the summation this leads to f(B+h) as the last added term – and yields a wrong result. The last term should be f(B–h), resp. f(B) with your method of computing this value twice. As already suggested, all this is better done with a simple loop counter.

Dieter

[Dieter]

(12-12-2015 05:44 PM)Thomas Klemm Wrote:  

(12-12-2015 04:24 PM)Namir Wrote:  

Code:

Loop Until B<A

This should be something like \(A\approx B\) or else you may end up with \(A=B+h\).

Yes, the exact condition would be A=B, but this cannot be used since roundoff errors can and will occur.
A usable exit condition could be

Code:

Loop Until a > b - h / 2

At least as long as h >> ULP(b)...

(12-12-2015 05:44 PM)Thomas Klemm Wrote:  I'd probably use \(N\) instead to control the loop.

Definitely, yes.

Dieter

[Namir]

Thank you all for your feedback. I have corrected my original post accordingly.

Namir

[Namir]

This is a slightly different version based on one of Dieter's remarks.

Code:

h=(B-A)/N
Sum=f(A)-f(B)
k=4
A=A+h
Do
  Sum = Sum + k*f(A)
  k=6-k
  A=A+h
  N=N-1
Loop Until N=0
Area=h/3*Sum

[Namir]

I looked at the Simpson's Rule listings in the Math Pacs for the HP-65 and HP-67. These listings could have used the approach in this thread to shorten the listings for these calculators! Sometimes it takes folks 30 or 40 years to think of shortcuts!

Namir

[Dieter]

(12-13-2015 01:52 AM)Namir Wrote:  I looked at the Simpson's Rule listings in the Math Pacs for the HP-65 and HP-67. These listings could have used the approach in this thread to shorten the listings for these calculators! Sometimes it takes folks 30 or 40 years to think of shortcuts!

It doesn't take 30 or 40 years for this. Just a bit of thinking and careful programming.

My impression of HP's software pacs, be it on cards for the 67/97, in ROM modules for the 41 or in their solution books, is that there are rarely examples that match contemporary programming standards. Often there is much, much room for improvement. One more reason for taking a look at what the users' libraries had to offer. ;-)

Dieter

[Dieter]

(12-12-2015 11:02 PM)Namir Wrote:  This is a slightly different version based on one of Dieter's remarks.

Code:

h=(B-A)/N
Sum=f(A)-f(B)
k=4
A=A+h
Do
  Sum = Sum + k*f(A)
  k=6-k
  A=A+h
  N=N-1
Loop Until N=0
Area=h/3*Sum

I still don't get why you want to compute f(b) twice and generate an additional function call that is not required. Why don't you do it the classic way?

Code:

h = (b - a) / n
sum = f(a) + f(b)
k = 4
Do
  a = a + h
  sum = sum + k * f(a)
  k = 6 - k
  n = n - 1
Loop Until n = 1
result = sum * h / 3

If you want to count down to n=0 simply add another n=n-1 on top of the do loop.

Dieter

[Namir]

(12-14-2015 09:04 PM)Dieter Wrote:  I still don't get why you want to compute f(b) twice and generate an additional function call that is not required. Why don't you do it the classic way?

To get good results the value of N has to be high. In this case one more call to f(x) does not make much of a difference. By contrast, saving one call to f(x) when N is low is not to help reducing the error.

Namir

[Dieter]

(12-15-2015 12:09 AM)Namir Wrote:  To get good results the value of N has to be high. In this case one more call to f(x) does not make much of a difference. By contrast, saving one call to f(x) when N is low is not to help reducing the error.

Sorry, if I'm a bit stubborn here, but the question is not whether one more call hurts or not, but where is the benefit of modifiying the original method so that another call is required, so that a value is added to the sum and subtracted again to compensate this error. One could just as well work with factors 8 and 4 instead of 4 and 2 and apply a correction by a factor of 0,5 afterwards. Doesn't make much sense either. ;-)

As far as accuracy is concerned: The error of Simpson's rule (or: Kepler's, to be more precise) can be estimated quite well. Here a method similar to the Romberg scheme can be applied: Calculate the integral for n and 2n intervals, then determine an extrapolated result via (16·I_2n – I_n)/15. This way an accuracy is achieved that is comparable to a Simpson approximation with a much higher number of iterations.

Example 1:
f(x) = 1/x, a=1, b=2
Exact integral = ln 2 = 0,6931471806

Result for n=24: I = 0,6931472743...
Result for n=48: I = 0,6931471864...
Extrapolation: I = 0,6931471806

A comparable accuracy would require a simple Simpson approximation with n=150...200.

Example 2:
f(x) = sin(x), a=0, b=Pi
Exact integral = 2

Result for n=40: I = 2,000000423...
Result for n=80: I = 2,000000026...
Extrapolation: I = 2,0000000000

A comparable accuracy would require a simple Simpson approximation with n=500.

Dieter

[Namir]

Here is a version that iterates one less time to calculate the integral:

Code:

Given function f(x), integration interval [A, B] and N divisions, where N is even.
To calculate the integral of f(x) for X=A to x=B using Simpson's Rule

h=(B-A)/N
Sum=f(A)+f(B)
C=4
Do
  A=A+h
  Sum = Sum + C*f(A)
  C=6-C
  N=N-1
Loop Until N=1
Area=h/3*Sum

The last iteration multiplies f(A) by 4 and leaves C with the value of 2. The initial value of Sum is initialized as f(A)+f(B). By controlling the loop with a counter, we can have a compact version that does not not do extra calculations.

Namir

[Dieter]

(12-16-2015 12:57 AM)Namir Wrote:  Here is a version that iterates one less time to calculate the integral:

Now remove the initial a = a+h and move the one within the loop right behind the Do, and you get what I suggested in my post #12, saving one line of code. ;-)

Dieter

[Namir]

(12-16-2015 01:25 PM)Dieter Wrote:  

(12-16-2015 12:57 AM)Namir Wrote:  Here is a version that iterates one less time to calculate the integral:

Now remove the initial a = a+h and move the one within the loop right behind the Do, and you get what I suggested in my post #12, saving one line of code. ;-)

Dieter

See updated last post.

Namir