Hello guys, I am working on a Simpson's rule problem. I am said to integrate from a=0 to b=1 at n=7 steps and should get an integral equal to 0.333. Also need to integrate from a=1 to b=0 with n=7 steps and get -0.333 for integral (Calculated manually). However, at the moment when I Integrate from 0 to 1 at 7 steps I get -0.486 or integrating from 1 to 0 at 7 steps give me integral -1.728. Please review my program if you can. I think the problem is somewhere in the function double simpson or when I call the function double simpson in main().
I looked over all the internet on how to solve the simpson's rule problem for uneven steps, but couldn't find anything programmed in C++. I know that it converts parabola in small squares, which means it should be even steps, but would love and need to know how to do it in uneven steps (in this case 7) for C++.
#include <iostream>
#include <cmath>
double simpson(double a, double b, int n);
double fk(double x);
int main()
{
// Data
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
printf("Integralas = %3.3f \n", f); // Rezultatas
return 0;
}
double simpson(double a, double b, int n)
{
double sum = fk(0);
double h = (b - a) / n; // base length
double h2 = h / 2; // middle base length between b and (b-a)/n
for (int i = 0; i < n; i++)
{
sum += (fk(a) + 4 * fk(a + h2) + fk(a + h)); // integral summation
}
// integral calculation
return h / 3 * sum - (fk(b) / 2);
}
double fk(double x)
{
return (x * x);
}
P.S. The difficulty arises when I have to compute uneven steps (in this case 7). I have no problem in solving the integral of f(x)=x^2 for even steps (2,4,6,8,etc.), but am struggling with uneven steps :/
Changed the simpson function, but the interval still is negative even tho it gets closer to 0.333 or -0.333 (-0.321 for 0 to 1 at 7 steps & -1.281 for 1 to 0 at 7 steps). Here how it looks:
1 2 3 4 5 6 7 8 9 10 11 12 13
double simpson(double a, double b, int n)
{
double sum = fk(0);
double h = (b - a) / n; // base length
double h2 = h / 2; // middle base length between b and (b-a)/n
for (int i = 0; i < n; i++)
{
if (i >= 1 && i <= 7) {sum += (fk(a) + 4 * fk(a + i*h2) + fk(a + i*h));} // integral summation
else {sum += (fk(a) + 4 * fk(a + h2) + fk(a + h));} //integral summation
}
// integral calculation
return h / 3 * sum - (fk(b) / 2);for (int i = 0; i < n; i++);
}
double simpson( double a, double b, int n )
{
double sum = 0.0;
double h = ( b - a ) / n; // Base length - assumes n "intervals". Weighting: 1:4:1 for start, middle, end of interval.
for ( int i = 0; i < n; i++ )
{
double left = a + i * h;
sum += fk( left ) + 4 * fk( left + 0.5 * h ) + fk( left + h ); // weighted sum for that interval
}
// Integral calculation
return sum * h / 6; // 6 = 1 + 4 + 1
}