title: Errors in Composite Simpson's Rule in C++
body: <p>I am writing a small program to approximate integrals in C++ using the composite Simpson's rule. It seems to work, but when I calculate the integral for different numbers of intervals (hence the <code>for</code> loop), some errors appear, as it seems completely at random. Here is my code.</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;cmath&gt;
  
  double f(double x)
  {
      double y = x;
      return y;
  }
  
  int main()
  {
      for(double n=1; n&lt;20; n+=1)
      {
          double a = 1;
          double b = 4;
          double d = (b-a)/n;
          double i = a;
          double area = 0;
          while(i&lt;b)
          {
              double j = i + d;
              area += (d/6) * ( f(i) + 4*f((i+j)/2) + f(j) );
              i = i + d;
          }
          std::cout &lt;&lt; "Order " &lt;&lt; n &lt;&lt; " : " &lt;&lt; area &lt;&lt; std::endl;
      }
      return 0;
  }
  </code></pre>
  
  <p>The output is as follows.</p>
  
  <pre><code>Order 1 : 7.5
  Order 2 : 7.5
  Order 3 : 7.5
  Order 4 : 7.5
  Order 5 : 7.5
  Order 6 : 7.5
  Order 7 : 9.30612
  Order 8 : 7.5
  Order 9 : 7.5
  Order 10 : 8.745
  Order 11 : 8.6281
  Order 12 : 7.5
  Order 13 : 7.5
  Order 14 : 7.5
  Order 15 : 7.5
  Order 16 : 7.5
  Order 17 : 8.22145
  Order 18 : 8.18056
  Order 19 : 7.5
  </code></pre>
  
  <p>I think this might have something to with variable types or storing problems, although I can't figure out what the problem is. Any help is appreciated.</p>