Simpson's Rule Implementation trick?

[Namir]

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

[Dieter]

(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

[Namir]

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

[Vtile]

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.

[Namir]

(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

[Vtile]

I see, back to the books then.

[ttw]

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.

[Dieter]

(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

[Paul Dale]

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

[Namir]

(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

[Namir]

(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

[Dieter]

(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