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

Edited 6 Years Ago by peter_budo: Keep it Organized - For easy readability, always wrap programming code within posts in [code] (code blocks)

Attachments .png 31.63 KB

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).

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!

Edited 3 Years Ago by Nick Evan: Fixed formatting

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.

This article has been dead for over six months. Start a new discussion instead.