Simpson's Method not producing correct results

[gdsturgill]

Hi, everyone!

I am trying to create a function that calculates an integral using Simpson's composite method but I am hitting a problem. When I run the program, I get an answer that is very different that what I think I should be getting.

I am using the equation f(x) = x^5 - 2x^2 + 3 on the interval [0,50], with n = 100 to 10,000 intervals. This integral should equal 1479316.6666667, but my results are as follows:

Actual Value = 1479316.6666667
For n = 100:
Simpsons = 539458
For n = 500:
Simpsons = 525717
For n = 1000:
Simpsons = 523967
For n = 5000:
Simpsons = 522562
For n = 10000:
Simpsons = 522386

Here is the code I am using:

float simpsons(float a, float b, int n)
{
	float sum1 = 0, sum2 = 0;
	//Find length of interval
	float h = (b - a) / n;

	//Find summations for even terms
	for (int i = 1; i <= (n / 2) - 1; i++)
	{
		if ((i % 2) == 0)
		{
			sum1 = sum1 + f(a + (i * h));
		}
	}
	//Find summation for odd terms
	for (int i = 1; i <= (n / 2); i++)
	{
		if ((i % 2) == 1)
		{
			sum2 = sum2 + f(a + 2 * (i * h));
		}
	}
	//Compute Simpson's Rule
	float s = (h / 3) * (f(a) + f(b) + (2 * sum1) + (4 * sum2));
	return s;
}

Any insight you can give is very much appreciated.

Thank you!

[TheIdeasMan]

Hi,

Prefer double rather than float.

Try using the debugger, see where things go wrong.

Please show the code for the f function.

Also it's a good idea to provide code we can compile :+)

[TheIdeasMan]

A float can only do 6 or 7 significant figures, so it's precision is easily exceeded.

On an interval of 50.0 and n = 100.0 the step size is 0.5. the first term is x⁵ , so when x > 1.0 it is already approaching the limit of the precision. So this has the effect of making the area of an individual interval a rectangle, rather than a trapezium reflecting the increase in the y value.

And it's get worse as n increases.

[gdsturgill]

Good call on checking the function. I didn't realize that it increased sharply starting,at about x=5. My function must not be hitting the higher values it needs to create an accurate integral. Need to change equation. Thanks for putting me on the right track!

[48T7M4Gy]

Another way using a single loop;

#include <iostream>


// function is x^5 -2x^2 + 3, x: 0->50
// integral is (1/6)x^6 - 2/3(x^3) + 3x
// (50^6)/6 - (2*(50^3)/3) + (3*50) = 2604083483.33

double f(const double);
double simpsons(double, double, int);

int main()
{
    double integral = simpsons( 0, 50, 100);
    std::cout << integral << '\n';
    
    return 0;
}

double f(const double x){
    double y = x*x * x*x * x - 2 * x*x  + 3;
    return y;
}

double simpsons(double a, double b, int n)
{
    double sum_even = 0, sum_odd = 0;
    
    //Find length of interval
    double h = (b - a) / n;

    for (int i = 1; i < n - 1; i += 2)
    {
        sum_odd += f(i * h);
        sum_even += f( (i+1) * h);
    }
    
    //Compute Simpson's Rule
    double s = h / 3 * ( f(a) + (4 * sum_odd) + (2 * sum_even) + f(b) );
    return s;
}

Definitely need double's!

[TheIdeasMan]

Good call on checking the function. I didn't realize that it increased sharply starting,at about x=5.

It's not the size of the y value that is the problem per say, it's the precision of the values both for low and higher values.

At x = 1.75, the integral y is 13.610066732 (9dp) , but a float can only do 13.61006 (7sf).

From kemorts code 2604083483.33, a float can only do 2604083000 (7sf).

When n = 10'000 , a lot of the individual strip areas will be the same.

So one can see how this affects the individual areas, and the final result.

[48T7M4Gy]

Another small point about Simpsons rule is it doesn't do much of a job, partly because it uses a quadratic approximation to the integral vs a 5th power in the function.

But there again, it was developed at a time when 100 slices would be a monumental hand-calculation and 10 slices would be tops giving an answer of "near-enough is good enough". We've moved on ...

_{If you ask a mathematician what 1+2 is you'll get a 12 page proof that the answer is 3.
If you ask a physicist the answer is 3.00 +- .01
If you ask an engineer the answer is "about 6".}