Problem with Numerical Integration

[faat99]

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

[Zjarek]

Double has precision close to 10^-16, your weird results are caused only by this. You also lose precision on every addition, specially in Trap function (that's why for n=10^10 you get e=10^-5).

[faat99]

thanks zjarek. Thats what I think as well. And now I attempt to design the program to use the max n that error is smallest. So I insert the Error of each n term into an array. Then Use a while loop, wanting to stop at when Next error starts to rebound. So..

double Sarr[20]; double mterm=1.0; i=1; Sarr[0]=10000;
do
{   mterm*=10;
    ES = 100*(Simp(a,b,mterm)-M)/M;
    Sarr[i-1]=ES; i++;
}while (Sarr[i]-Sarr[i-1]<0);

But this doesn't work at all. I guess it is because the compiler can't regonise the "i" in the while condition. I want to sychronise the condition i and the loop i. Anyone please tell me how? Thank you!

[DavidB]

Right now you define "i" on the fly, within the for-loop.

Why not just define it outside the loop, say, at the beginning of the subroutine where you define the rest of the variables.