The most compact algorithm for Simpson's rule??

12122015, 04:24 PM
(This post was last modified: 12122015 10:58 PM by Namir.)
Post: #1




The most compact algorithm for Simpson's rule??
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:
Many implementations use two summationsone for odd terms and one for even terms. In addition these implementations use two loops with slightly different ranges. The above pseudocoode implements a simple version that uses one loop and a few simple computational tricks. Namir 

12122015, 05:44 PM
Post: #2




RE: The most compact algorithm for Simpson's rule??  
12122015, 08:28 PM
Post: #3




RE: The most compact algorithm for Simpson's rule??
Using N should also work.


12122015, 08:29 PM
(This post was last modified: 12122015 08:29 PM by Namir.)
Post: #4




RE: The most compact algorithm for Simpson's rule??  
12122015, 09:05 PM
(This post was last modified: 12122015 09:25 PM by Dieter.)
Post: #5




RE: The most compact algorithm for Simpson's rule??
(12122015 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 = (ba)/n (12122015 04:24 PM)Namir Wrote: Many implementations use two summationsone 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 

12122015, 09:16 PM
(This post was last modified: 12122015 10:14 PM by Dieter.)
Post: #6




RE: The most compact algorithm for Simpson's rule??
(12122015 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 

12122015, 09:44 PM
(This post was last modified: 12122015 09:44 PM by Dieter.)
Post: #7




RE: The most compact algorithm for Simpson's rule??
(12122015 05:44 PM)Thomas Klemm Wrote:(12122015 04:24 PM)Namir Wrote: 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)... (12122015 05:44 PM)Thomas Klemm Wrote: I'd probably use \(N\) instead to control the loop. Definitely, yes. Dieter 

12122015, 10:54 PM
Post: #8




RE: The most compact algorithm for Simpson's rule??
Thank you all for your feedback. I have corrected my original post accordingly.
Namir 

12122015, 11:02 PM
Post: #9




RE: The most compact algorithm for Simpson's rule??
This is a slightly different version based on one of Dieter's remarks.
Code: h=(BA)/N 

12132015, 01:52 AM
(This post was last modified: 12132015 01:53 AM by Namir.)
Post: #10




RE: The most compact algorithm for Simpson's rule??
I looked at the Simpson's Rule listings in the Math Pacs for the HP65 and HP67. 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 

12132015, 01:27 PM
Post: #11




RE: The most compact algorithm for Simpson's rule??
(12132015 01:52 AM)Namir Wrote: I looked at the Simpson's Rule listings in the Math Pacs for the HP65 and HP67. 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 

12142015, 09:04 PM
Post: #12




RE: The most compact algorithm for Simpson's rule??
(12122015 11:02 PM)Namir Wrote: This is a slightly different version based on one of Dieter's remarks. 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 If you want to count down to n=0 simply add another n=n1 on top of the do loop. Dieter 

12152015, 12:09 AM
(This post was last modified: 12152015 01:14 AM by Namir.)
Post: #13




RE: The most compact algorithm for Simpson's rule??
(12142015 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 

12152015, 07:12 AM
(This post was last modified: 12152015 07:30 AM by Dieter.)
Post: #14




RE: The most compact algorithm for Simpson's rule??
(12152015 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 

12162015, 12:57 AM
(This post was last modified: 12162015 01:53 PM by Namir.)
Post: #15




RE: The most compact algorithm for Simpson's rule??
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. 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 

12162015, 01:25 PM
(This post was last modified: 12162015 01:27 PM by Dieter.)
Post: #16




RE: The most compact algorithm for Simpson's rule??  
12162015, 01:54 PM
Post: #17




RE: The most compact algorithm for Simpson's rule??
(12162015 01:25 PM)Dieter Wrote:(12162015 12:57 AM)Namir Wrote: Here is a version that iterates one less time to calculate the integral: See updated last post. Namir 

« Next Oldest  Next Newest »

User(s) browsing this thread: