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
		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++)            
		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;
	dx = (b-a)/n;       
	for ( int i=1; i<=n; i++ )     //loop sum = yo+y1+y2+y3+...+yn              
	{	x+=dx;
	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
	for ( int i=1; i<n/2; i++ )     //loop q = y2+y4+y6+...+y|n-2
	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 by peter_budo: Keep it Organized - For easy readability, always wrap programming code within posts in [code] (code blocks)

Attachments .png 31.63 KB
7 Years
Discussion Span
Last Post by DavidB

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;
{   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 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.
Take the time to help us to help you. Please be thoughtful and detailed and be sure to adhere to our posting rules.