### Simpson's rule question

Hey guys, I know that Simpson's rule has to be applied for even n rectangles, however in my assignment I am asked to calculate in 7 rectangles between a = 0 and b = 1. I have read forum links on this problem, but couldn't find the answer. Therefore, I am asking what I should change in my code to get a correct answer of integral = 0.333 for 7 rectangles (Currently I get 0.293). When I calculate for 12, 10, 8, etc with a = 0 and b = 1 everything is fine - I get integral = 0.333.

Here is my code:

 ``1234567891011121314151617181920212223242526272829303132333435363738394041424344454647`` ``````#include #include double simpson(double a, double b, int n); double fk(double x); int main() { // Database double a, b, f; int n; printf("a = "); scanf("%lf", &a); printf("b = "); scanf("%lf", &b); printf("n = "); scanf("%d", &n); f = simpson(a,b,n); // Integral // Result printf("Integral = %3.3f \n", f); return 0; } double simpson(double a, double b, int n) { double sum_even = 0, sum_odd = 0; double h = (b - a) / n; // finding h base length // summation of all sums for (int i = 1; i < n; i+=2) { sum_odd += 4*fk(i*h); } for (int i = 2; i < n; i+=2) { sum_even += 2*fk(i*h); } // counting of integral double result = h / 3 * (fk(a) + sum_odd + sum_even + fk(b)); return result; } double fk(double x) { return x * x; }``````
Last edited on
What does your assignment actually say?

(BTW There aren't actually any "rectangles" anywhere in the problem).
If you desperately want your code to work with odd values of n then add the line

` if ( n % 2 != 0 ) result += ( h / 12 ) * ( fk(b) + 4 * fk(b-h) - fk(b-2*h) );`

between lines 40 and 41.
Last edited on
lastchance, thank you. Why do you get h / 12? The simpson's rule formula says it's h/3. Also why you are using fk(b-h) and then minus fk(b-2*h)?

My asssignment says to create a program that calculates integral by simpson's rule which I have done. And then they ask to calculate integral for n = 7, a = 0 and b = 0.
Have you tried it, @Aurimas?

If you confirm it works then I will explain where it comes from. Yes, the h/12 is intentional.
Yeah. I tried it. It works well. I get 0.333. Can you now give me the answers to my questions in the last post?:)
One way of deriving Simpson's rule is to fit a quadratic through 3 points and then integrate over a double interval of length 2h. If n is odd then you have one interval of length h left over. In this case, if you fit a quadratic through the last 3 points (b, b-h and b-2h) but only integrate it over the last h interval you get (after a lot of algebra) a contribution
(h/12)(5f(b)+8f(b-h)-f(b-2h))
Then I have to subtract off what your code already includes for this interval. The way that you have done the sums gives
(h/3)(f(b-h)+f(b))
The first f comes from half of the sum_even contribution and the second from the end of the interval.

If you subtract these (carefully) you will get the additive correction that I gave.

WARNING: this depends very much on how your original code treats the last interval, and different implementations will require different multipliers.

I have no idea whether this is the recommended technique - I simply hashed it out on the train on the way to work. However, it is still second-order and will give the exact result for a quadratic test function ... as you have here.

The reason I asked about your assignment is that I think it slightly more like that you might be required to do a check that n was even.
closed account (48T7M4Gy)
Why not bisect the last (uncounted) slice and apply Simpsons rule to those 2 subslices (3 values), and add the result to the rest?
 Why not bisect the last (uncounted) slice and apply Simpsons rule to those 2 subslices (3 values), and add the result to the rest?

` if ( n % 2 != 0 ) result += ( h / 6 ) * ( -fk(b-h) + 4 * fk(b-0.5*h) - fk(b) );`
and it produces the correct answer.

Depends whether this now constitutes 8 (unequal) intervals instead of 7, I guess.
Last edited on
closed account (48T7M4Gy)
I was thinking this, but it doesn't work :(

 ``123456789`` ``````if(n % 2 == 0) { f = simpson(a,b,n); // Integral } else { double dx = (b - a)/n; f = simpson(a, b-dx, n-1) + simpson(b-dx, b, 2); }``````
Thank kemort. Any other ideas how I could get 0.333 for odd n like 7?
closed account (48T7M4Gy)
Yeah, I'll post it in a second. It is similar to lastchance but a slightly different approach in that I take one slice off and put it back again as a new slice the same width as all the others, but it is separately split in two.

Thanks kemort. I'll wait for response.
closed account (48T7M4Gy)
 ``12345678910111213141516171819202122232425262728293031323334353637383940`` ``````#include double simpson(double, double, int); double y(double); int main() { double a = 0.0, b = 1.0, f = 0; int n = 3; double dx = (b-a)/n; if(n%2 == 0) f = simpson(a, b, n); else f = simpson(a, b-dx, n - 1) + dx/6*( y(b-dx) + 4*y((2*b-dx)/2) + y(b) ); // Result std::cout << "Integral = " << f << '\n'; return 0; } double simpson(double x0, double x2, int n) { double sum_even = 0, sum_odd = 0; double h = (x2 - x0) / n; // finding h base length // summation of all sums for (int i = 1; i < n; i += 2) { sum_odd += 4*y(i*h); } for (int i = 2; i < n; i += 2) { sum_even += 2*y(i*h); } // counting of integral return h / 3 * (y(x0) + sum_odd + sum_even + y(x2)); } double y(double x) { return x * x; }``````

(Note: n >= 2)

4*y((2*b-dx)/2) or more tidily should be (obviously) 4*y(b-dx/2) because the slice width is just dx/2.

Last edited on
Thanks again kemort.
Topic archived. No new replies allowed.