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 by peter_budo: Keep it Organized - For easy readability, always wrap programming code within posts in [code] (code blocks)

Attachments
3
Contributors
3
Replies
4
Views
8 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;
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 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 topic has been dead for over six months. Start a new discussion instead.
Have something to contribute to this discussion? Please be thoughtful, detailed and courteous, and be sure to adhere to our posting rules.