 Simpson's Rule Implementation trick? - 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: Simpson's Rule Implementation trick? (/thread-5886.html) Simpson's Rule Implementation trick? - Namir - 03-18-2016 12:11 AM I found C code in an article on the Internet that implements Simpson's rule using a single for loop. Here is the pseudo-code for the C code: Code: ``` Simpson's Rule for Numerical Integration Ver 1 ====================================== Given f(x) and limit [a, b} with n points if n modulo 2 = 0 then n = n + 1 h = (b - a) / (n + 1) sum = f(a) + f(b) + 4*f(a+h) for i=2 to n step 2   sum = sum + 2*f(a+i*h) + 4*f(a + (i+1)*h) next area = h/3*sum``` Here is a version that I modified to avoid using a+i*h and a+(i+1)*h: Code: ``` Simpson's Rule for Numerical Integration Ver 2 ====================================== Given f(x) and limit [a, b} with n points if n modulo 2 = 0 then n = n + 1 h = (b - a) / (n + 1) x = a+h sum = f(a) + f(b) + 4*f(x) for i=2 to n step 2   x = x + h   sum = sum + 2*f(x)   x = x + h   sum = sum + 4*f(x) next area = h/3*sum``` Consolidating the for loop to have a single update to sum: Code: ```Simpson's Rule for Numerical Integration Ver 3 ============================================== Given f(x) and limit [a, b} with n points if n modulo 2 = 0 then n = n + 1 h = (b - a) / (n + 1) x = a+h sum = f(a) + f(b) + 4*f(x) c = 2 for i=2 to n   x = x + h   sum = sum + c*f(x)   c = 6 - c next area = h/3*sum``` Enjoy! Namir RE: Simpson Rule Implementation trick? - Dieter - 03-18-2016 01:14 PM (03-18-2016 12:11 AM)Namir Wrote:  I found C code in an article on the Internet that implements Simpson's rule using a single for loop. Fine. But what's the improvement over another solution with a simple for loop as it has been discussed in December in this thread? Dieter RE: Simpron Rule Implementation trick? - Namir - 03-18-2016 08:47 PM By initializing the sum as: Code: `sum = f(a) + f(b) + 4*f(a+h)` Instead of the usual expression: Code: ``` sum = f(a) + f(b)``` We can then iterate over the terms multiplied by 2 and then by 4. We can do this: Code: ```for i=2 to n step 2   sum = sum + 2*f(a+i*h) + 4*f(a + (i+1)*h) next``` Or, borrowing from the old threads: Code: ```c = 2 for i=2 to n   x = x + h   sum = sum + c*f(x)   c = 6 - c next``` RE: Simpron Rule Implementation trick? - Vtile - 03-18-2016 11:15 PM Why Simpson rule, why not ie. runge-kutta as isn't that the way when calculating numerical solution for integrals with high accuracy? I know must be rather idiotic question (since obviously it wouldn't discussed here if it would not have use), but need to ask to learn. RE: Simpron Rule Implementation trick? - Namir - 03-18-2016 11:53 PM (03-18-2016 11:15 PM)Vtile Wrote:  Why Simpson rule, why not ie. runge-kutta as isn't that the way when calculating numerical solution for integrals with high accuracy? I know must be rather idiotic question (since obviously it wouldn't discussed here if it would not have use), but need to ask to learn. The Gauss-Krnorod quadrature is among the best numerical integration methods. This thread just comments on a little trick that I noticed in an article that discusses basic numerical integration methods. The code that initializes the sum of terms caught my attention .... that's all. BTW, The Runge-Kutta methods are for solving ordinary differential equations. Namir RE: Simpron Rule Implementation trick? - Vtile - 03-19-2016 01:14 AM I see, back to the books then. RE: Simpron Rule Implementation trick? - ttw - 03-19-2016 04:59 AM Romberg's method is pretty good. It can be of arbitrarily high order if the integrand is smooth enough and can be adapted to giving error estimates. RE: Simpron Rule Implementation trick? - Dieter - 03-19-2016 06:45 AM (03-19-2016 04:59 AM)ttw Wrote:  Romberg's method is pretty good. It can be of arbitrarily high order if the integrand is smooth enough and can be adapted to giving error estimates. This is also possible with Simpson's method. Two approximations with n and 2n intervals give a new and improved estimate. Take a look a this thread in the HP41 software library. Dieter RE: Simpron Rule Implementation trick? - Paul Dale - 03-19-2016 06:50 AM I'm a fan of the Gauss-Kronrod quadratures. The 34S started with one of these for its integrate function. It wasn't adaptive and failed too often and people complained so I went back to a Romberg method. For the 43s I'd like to switch back to an adaptive GK method. It should have the memory to handle it. - Pauli RE: Simpron Rule Implementation trick? - Namir - 03-19-2016 02:21 PM (03-19-2016 06:45 AM)Dieter Wrote:   (03-19-2016 04:59 AM)ttw Wrote:  Romberg's method is pretty good. It can be of arbitrarily high order if the integrand is smooth enough and can be adapted to giving error estimates. This is also possible with Simpson's method. Two approximations with n and 2n intervals give a new and improved estimate. Take a look a this thread in the HP41 software library. Dieter I have implemented a variant of the Romberg method that uses Simpson's rule to start with instead of teh Trapezoidal rule. You can find this variant method on my web site. Namir RE: Simpron Rule Implementation trick? - Namir - 03-19-2016 02:22 PM (03-19-2016 06:50 AM)Paul Dale Wrote:  I'm a fan of the Gauss-Kronrod quadratures. The 34S started with one of these for its integrate function. It wasn't adaptive and failed too often and people complained so I went back to a Romberg method. For the 43s I'd like to switch back to an adaptive GK method. It should have the memory to handle it. - Pauli I think the Integrate function of the HP-34C and HP-15C (and later models?) use a method based on Romberg's method. Namir RE: Simpson Rule Implementation trick? - Dieter - 03-19-2016 05:05 PM (03-19-2016 02:21 PM)Namir Wrote:  I have implemented a variant of the Romberg method that uses Simpson's rule to start with instead of teh Trapezoidal rule. You can find this variant method on my web site. Yes, I have seen this some time ago and I tried some of the VBA code (after some modifications for my older Excel version). Interesting results, the Romberg-Simpson method looks quite promising. However, on a calculator with limited ressources I prefer the Simpson method with additional extrapolation. Now the next step would be an improvement for periodic functions, i.e. something that uses nodes that are not evenly distributed. Such a method is also used in the original 34C Integrate method. BTW... Namir, could you please correct the typo in the subject of your initial post? The way it is written now ("Simpron") noone will be able to find this thread while searching the forum. Dieter