Simpson's rule question

[aurimas13]

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:

#include <iostream>
#include <stdio.h>

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;
}

[lastchance]

What does your assignment actually say?

(BTW There aren't actually any "rectangles" anywhere in the problem).

[lastchance]

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.

[aurimas13]

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.

[lastchance]

Have you tried it, @Aurimas?

If you confirm it works then I will explain where it comes from. Yes, the h/12 is intentional.

[aurimas13]

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?:)

[lastchance]

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.

[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?

[lastchance]

Why not bisect the last (uncounted) slice and apply Simpsons rule to those 2 subslices (3 values), and add the result to the rest?

After subtracting off what the code effectively had for this interval already it leads to the alternative additive correction
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.

[48T7M4Gy]

I was thinking this, but it doesn't work :(

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);
    }

[aurimas13]

Thank kemort. Any other ideas how I could get 0.333 for odd n like 7?

[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.

[aurimas13]

Thanks kemort. I'll wait for response.

[48T7M4Gy]

#include <iostream>

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.

[aurimas13]

Thanks again kemort.