I have done approximations of I=∫dx/x using Simpson's rule and Trapezium rule, and I output the results to excel. And I plotted a graph. My problem is, why are there fluctation at a point of n instead of showing linearity all the way(Picture below)? I have been thinking this all the time and still having no idea. Please help Thanks! My program is below.
//Numerical Integration of Simpson’s rule and Trapezium rule
#include <fstream>
#include <cmath>
using namespace std;
double Trap(double, double, int);
double Simp(double, double, int);
double f(double);
int main()
{
ofstream file1("Error of trap.txt");
double a=1.0, b=10.0, ET, ES;
double M = log(b) - log(a); //acutal value by integration
double nterm=1.0; //no. of strp
file1 << "log(n)" << "\t" << "log(error)" << "\t" << endl;
for (int i=0; i<10; i++) //use loop to increase no. of strips
{
nterm*=10;
ET = 100*(Trap(a,b,nterm)-M)/M;
file1 << log10(nterm) << "\t" << log10(abs(ET)) << "\t" << endl;
}
ofstream file2("Error of sim.txt");
double mterm=1.0;
file2 << "log(n)" << "\t" << "log(error)" << "\t" << endl;
for (int i=0; i<10; i++)
{
mterm*=10;
ES = 100*(Simp(a,b,mterm)-M)/M;
file2 << log10(mterm) << "\t" << log10(abs(ES)) << "\t" << endl;
}
return 0;
}
//functions of Trapezium rule and Simpson's rule
double Trap(double a, double b, int n) //function of Trapezium rule
{
double sum, y, x, dx;
x=a;
sum=f(x);
dx = (b-a)/n;
for ( int i=1; i<=n; i++ ) //loop sum = yo+y1+y2+y3+...+yn
{ x+=dx;
sum+=f(x);
}
y = dx*(sum - f(a)/2 - f(b)/2); //Trapezium rule
return y;
}
double Simp(double a, double b, int n) //functions of Simpsons's rule
{
if((n%2!=0) || (n<4.0)) {exit(99);} //limitation
double p=0.0, q=0.0; //loops for Simpson's rule
double dx = (b-a)/n;
for ( int i=0; i<n/2; i++ ) //loop p = y1+y3+y5+...+y|n-1
{p+=f(a+(2*i+1)*dx);}
for ( int i=1; i<n/2; i++ ) //loop q = y2+y4+y6+...+y|n-2
{q+=f(a+2*i*dx);}
double y = dx*( f(a) + 4*p + 2*q + f(b) ) /3; //Simpson's rule
return y;
}
double f(double x) //The integrand
{ return 1/x; }