3

This code is using Simpson's rule to calculate the integral of x*sin(x) with the bounds of (1,2). The problem I am having is, while its getting very close to the actual value. Even with 999 iterations, it still doesn’t hit the point. While I have a separate program that uses trapezoidal rule for the same thing, and it DOES hit the point exactly after 1000 iterations. The point it should hit is "1.440422"

Is that what should be happening for Simpson's rule? Or is there something wrong with my code?

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
double f(double x);
int main()
{
  double x,result,y,z,h,s1,s2;
  s1 = 0;
  s2 = 0;
  int i,n;
  printf("\nHow many points to you want it evaluated at (odd number)? And what are the bounds? lower bound,upper bound >\n");
  scanf("%d %lf,%lf",&n,&y,&z);
  h = (z-y)/n;
  result = 0;
  if(n%2!=0)
    {
      for(i=0;i<n;i++)
    {
      if(i%2==0)
        {
          s1 = s1+f(y+i*h);
        }
      else
        {
          s2 = s2+f(y+i*h);
        }
    }
      result = (h/3)*(f(y)+f(z)+4*s2+2*s1);
      printf("\nThe value is %lf with %d interations\n",result,i);
    }
  else 
    {
      printf("\n The number of points has to be odd, try again\n");
    }
}


double f(double x)
{
  return(x*sin(x));
}
user3908631
  • 63
  • 1
  • 5
  • Welcome to the world of approximations and floating point arithmetic. If you use an approximate formula, you get an approximate result. It is probably a fluke that the trapezium rule ends up with 'exactly' the right answer. – Jonathan Leffler Oct 14 '15 at 15:32
  • Yeah thats the part that was confusing me was that I figured Simpson's approximations would have less error than the trapezoidal. But I guess it could just be a fluke with the function I am using. Thank you – user3908631 Oct 14 '15 at 15:36
  • 1
    With 99999 iterations I get your required value of `1.440422`. Some functions converge faster than others, for example there are several ways to calculate `pi` but not all are useful. – Weather Vane Oct 14 '15 at 16:43

1 Answers1

3

The problem you are seeing is probably because of the format string used to read the numbers.

scanf("%d %lf,%lf",&n,&y,&z);
         // ^^^ Is the , there on purpose?

Try removing that , from the format string and see if the problem goes away.

It cannot be emphasized enough - Always check the return value of scanf.

if ( scanf("%d %lf %lf", &n, &y, &z) != 3 )
{
   // Deal with error.
}

To make sure that the numbers read are accurate, add a line that echos the input back to stdout.

printf("n: %d, y: %lf, z: %lf\n", n, y, z);

I noticed couple of errors in your code:

  1. The interval h is not right. Since you are using n points, there are n-1 intervals. Hence, h has to be:

    h = (z-y)/(n-1);

  2. Since you are adding f(y) and f(z) in the very last statement, the loop has to be:

    // Not good to use for(i=0;i<n;i++)
for(i=1;i<n-1;i++)
{

With those fixes, I get the output of 1.440422 using n = 1001.

R Sahu
  • 204,454
  • 14
  • 159
  • 270