 The most compact algorithm for Simpson's rule?? - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: The most compact algorithm for Simpson's rule?? (/thread-5318.html) The most compact algorithm for Simpson's rule?? - Namir - 12-12-2015 04:24 PM 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 RE: The most compact algorithm for Simpson's rule?? - Thomas Klemm - 12-12-2015 05:44 PM (12-12-2015 04:24 PM)Namir Wrote:   Code: `Loop Until BB 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 RE: The most compact algorithm for Simpson's rule?? - Dieter - 12-12-2015 09:44 PM (12-12-2015 05:44 PM)Thomas Klemm Wrote:   (12-12-2015 04:24 PM)Namir Wrote:   Code: `Loop Until B 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 RE: The most compact algorithm for Simpson's rule?? - Namir - 12-12-2015 10:54 PM Thank you all for your feedback. I have corrected my original post accordingly. Namir RE: The most compact algorithm for Simpson's rule?? - Namir - 12-12-2015 11:02 PM 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``` RE: The most compact algorithm for Simpson's rule?? - Namir - 12-13-2015 01:52 AM 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 RE: The most compact algorithm for Simpson's rule?? - Dieter - 12-13-2015 01:27 PM (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 RE: The most compact algorithm for Simpson's rule?? - Dieter - 12-14-2015 09:04 PM (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 RE: The most compact algorithm for Simpson's rule?? - Namir - 12-15-2015 12:09 AM (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 RE: The most compact algorithm for Simpson's rule?? - Dieter - 12-15-2015 07:12 AM (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·I2n – In)/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 RE: The most compact algorithm for Simpson's rule?? - Namir - 12-16-2015 12:57 AM 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 RE: The most compact algorithm for Simpson's rule?? - Dieter - 12-16-2015 01:25 PM (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 RE: The most compact algorithm for Simpson's rule?? - Namir - 12-16-2015 01:54 PM (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