## 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 2

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 44

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.