Numerical Integration Report

cout <<"Enter the Lower limit"<<endl;

cin >> a;

cout <<"Enter the upper limit"<<endl;

cin >> b;

cout <<"Enter the MAX number of strips"<<endl;

cin >> nmax;

The number of strips to start off with was then defined. The trapezium rule, started with n=1, one strip, and Simpson’s Rule started with n=4, 4 number of strips.

A do…while loop, then began which calculated the Trapezium or Simpson’s rule, for each value of n, up to the maximum number of strips specified at the beginning by the user.

Also, dx was defined, right at the beginning of each code, as this is constant each time;

dx = (b-a)/(n);//strip height

Within this loop, a code which calculated the integral was designed. This is where the coding differed for the Trapezium Rule and Simpson’s Rule.

Recalling the formula for the Trapezium rule,

I = dx/2 x (y0 + 2(y1 + y2 + y3….yn-1) + yn);

I decided to calculate the part in the brackets using a do..while loop.

An integer counter was put into place, which ensured that the computer does not go through the loop too often, and stops at the right time. The letter “i”, was used for this.

Below is the part of the code which had the loop.

i=1; //integer counter

        do //loop for trapezium rule

        {

                y = 2*(1/x);

                sum += y;

                x += dx;

                i++;

        }

        while (i<=n);

Each time the y value was calculated, it was stored into sum, thus satisfying the part of the equation;

…..2(y1 + y2 + y3….yn-1)…..

Once the loop has been completed the appropriate number of times, the integral was calculated by simply completing the formula as it is in the Trapezium Rule formula.

        I = (dx/2)* ((1/a) + sum + (1/b));

1/a and 1/b represent the first and the last ordinates, y0 and yn.

The error value was then calculated, by simply subtracting the exact value , ln10, from the value computed in the code.

The number of strips, the integral and the error value were then stored into an external .txt file. This was done using the ofstream function;

trapezium<< n<<"\t"<<I <<"\t"<<error<<endl;

Once the “sum” was reset, one loop was completed.

A value of ten was then added to n, at the end of the loop, in order to increase the value of n by ten each time.

This prompted the start of the loop once more, which only ended once the value of n, was less than the maximum number of strips specified at the beginning.

This part of the code looked as follows;

 

        sum = 0.0; //reset sum

        n+=10; //add ten strips each time

         }

        while (n<nmax);

        return 0;

For Simpson’s Rule, the main body where the integral was calculated presented more of a challenge.

In principle, as you can see from the formula

I = dx/3 x (y0 + yn + 4(y1 + y3 +….. + yn-1) + 2(y2 + y4 +…..+ yn-2)

The Simpson’s Rule operates by calculating the sum of the first and last ordinates, y0 and yn, added to FOUR times the sum of the odd ordinates, plus TWICE the sum of the even ordinates.

Therefore two independent loops were used.

do // loop to calculate sum of all the odd ordinates

        {

                y = 4*(1/x);

                sum1 += y;

                x += 2*dx; // ensures that the next x value used remains odd

                c++;

        }

        while (x<b);

Above is the first loop, which calculated the sum of the odd ordinates, and stores them into sum1.

The integer variable, c , is an integer counter which has no effect on the calculation in the loop. It is purely present as a ‘’check’’ to ensure that the loop is carried out the correct number of times.

The loop is done until x reaches the upper limit of integration, b.

The second loop looks distinctly similar, however the following is added before it starts, to ensure that the first ordinate used is even.

x1 = a+dx+dx; //ensures only even ordinates are used

The loop that follows is

do // loop to calculate sum of all the even ordinates

        {

                z = 2*(1/x1);

                sum2 += z;

                x1 += 2*dx; // ensures that the next x value used remains even

                d++;

        }

        while (x1<b);

As previously, the loop ends once x1 reaches the upper limit of integration, b.

The integer counter in this case is d.

Once the computer has gone through both loops, the final value of the Integral, as well as the error is calculated, as follows.

        I = (dx/3)* ((1/a) + sum1 + sum2 + (1/b)); //Intergral, I

        error = I - log(10.0); //error calculation

Again, 1/a and 1/b represent the first and the last ordinates.

Like with the trapezium rule, the strips, the Integral and the error were all saved into a .txt file.

simpsons <<n<<"\t"<< I<<"\t"<<error<<endl;

Commented out, is a cout, which can be used to check if the loops were carried out the correct number of times.

//cout << c << "\t"<<d<<endl; // cout to see if the right number of ordinates are used

Sum1, and sum2 were then reset, ten was once again added to n, and the calculation for Simpson’s Rule was carried out once more, until n reaches the maximum number of strips specified at the start, nmax.

        sum1 = 0.0;//reset sum1

        sum2 = 0.0;//reset sum2

        n+=10; // add ten strips each time

        

        }

        while (n<=nmax);

The algorithms designed before implementing them into c++ code, can be found in the appendix.

Results, Errors, Discussion;

The relationship between the error and the number of strips used is

Error = 1/nm

Where n = the number of strips and m is equal to a power.

In order to calculate the power m, a log-log graph must be plotted.

By taking the log of both sides, the above equation can be transformed into the form y=mx + c.

ln(Error) = ln (1) – ln(nm)

ln(Error) = o – mln(n)

y            = c + mx

Hence, ln(Error) would be plotted on the y axis, and ln(n) on the x axis, with the negative of the gradient giving the value of m.

The code was compiled and built, with the maximum number of strips set to 1 million, in order to give a large spread of data.

The results were stored into the .txt files and then plotted onto a graph.

The gradient of the trapezium rule graph came out to near enough 1.

This means that the relationship between the error and the number of strips used is inverse.

From the graph it is visible that up to roughly 2980 strips, there is a straight line graph. However after that it can be seen that the points begin to appear far more random, as the number of strips goes way past 3000, the more the number of strips used, the larger the uncertainty in the error becomes. There is near enough no correlation after about 3000 strips.

For the first part of the graph, the gradient calculated comes to roughly 3.75.

From this, it can be deduced that Simpson’s rule is a far more accurate method of numerical integration.

This can be checked by as follows.

For 2000 strips, the error value for the trapezium rule will be a factor of;

Error = 1/20001

= 0.0005

For Simpson’s rule, this will be;

Error = 1/20003.75

= 4x10-13

There is an enormous difference between these two accuracies.

However as mentioned above, once the number of strips exceeds roughly 3000, there begins to be no correlation between the number of strips and the error.

To examine this more closely, a graph showing the number of strips between 1,000,000 and 10,000,000 will be plotted.

There is no correlation from this graph. However it is evident that the accuracy, once the number of strips exceeds roughly 3000 becomes very random, and as the number of strips increases, the margin in which the error is found gradually increases also.

Conclusion:

The Trapezium Rule and Simpson’s Rule are two methods of numerical integration which are fairly good estimates. The more the number of strips used, the more accurate the result is. A thorough examination of both methods was carried out for large values of n, number of strips, and graphs were plotted in order to compare the degree of accuracy of the two methods.

It was found that Simpson’s Rule is a far more accurate method of numerical integration. However, for extremely large values of n, exceeding 3000, the error can not be accurately computed, as it becomes more and more random for larger number of strips. Nevertheless, even though this is the case, the result is still significantly more accurate than for The Trapezium Rule.

Further study could be carried out and even larger values of “n” could be used to further test this theory.

Bibliography:

 

 

Appendix:

TRAPEZIUM RULE

//Trapezium rule

#include <iostream>

#include <cmath>

#include <fstream>

 using namespace std;

int main ()

{

        

        double sum = 0.0;

        double a, b,dx, x, I, n,y,nmax,error;

        int i;

        cout <<"Enter the Lower limit"<<endl;

        cin >> a;

        cout <<"Enter the upper limit"<<endl;

        cin >> b;

        cout <<"Enter the MAX number of strips"<<endl;

        cin >> nmax;

        ofstream trapezium("trapezium.txt");

        trapezium.precision(20);//precision to 20dp

        n=1; //start with n=1 strips

        do

        {

        dx = (b-a)/(n);//strip height

        x=a+dx;

        i=1; //integer counter

        do //loop for trapezium rule

        {

                y = 2*(1/x);

                sum += y;

                x += dx;

                i++;

        }

        while (i<=n);

        

        I = (dx/2)* ((1/a) + sum + (1/b));

        

        error = I - log(10.0);//error calculation

        

        trapezium<< n<<"\t"<<I <<"\t"<<error<<endl;

        sum = 0.0; //reset sum

        n+=10; //add ten strips each time

         }

        while (n<nmax);

        return 0;

        }

SIMPSONS  RULE

//Simpsons Rule

#include <iostream>

#include <cmath>

#include <fstream>

 using namespace std;

int main ()

{

        ofstream simpsons("simpsons rule.txt");

        simpsons.precision(20);//make precision to 20digits

        double sum1 = 0.0;

        double sum2 = 0.0;

        double a, b,dx, x, I,nmax,y,x1, z,error;

        int n,c,d;

        cout <<"Enter the Lower limit"<<endl;

        cin >> a;

        cout <<"Enter the upper limit"<<endl;

        cin >> b;

        cout <<"Enter the max number of strips"<<endl;

        cin >> nmax;

        n= 4; // start with n=4 strips

        do //loop for continuing number of strips up to N, max number of strips

        {

        dx = (b-a)/(n); //strip height

        x=a+dx;

        c=0; //integer counter check. to see if it is taking the right number of ODD ordinates

        do // loop to calculate sum of all the odd ordinates

        {

                y = 4*(1/x);

                sum1 += y;

                x += 2*dx; // ensures that the next x value used remains odd

                c++;

        }

        while (x<b);

        x1 = a+dx+dx; //ensures only even ordinates are used

        d = 0; ////integer counter check. to see if it is taking the right number of EVEN ordinates

        do // loop to calculate sum of all the even ordinates

        {

                z = 2*(1/x1);

                sum2 += z;

                x1 += 2*dx; // ensures that the next x value used remains even

                d++;

        }

        while (x1<b);

        I = (dx/3)* ((1/a) + sum1 + sum2 + (1/b)); //Intergral, I

        error = I - log(10.0); //error calculation

        //cout << c << "\t"<<d<<endl; // cout to see if the right number of ordinates are used

        simpsons <<n<<"\t"<< I<<"\t"<<error<<endl;

        sum1 = 0.0;//reset sum1

        sum2 = 0.0;//reset sum2

        n+=10; // add ten strips each time

        

        }

        while (n<=nmax);

        return 0;

}

}

Trapezium Rule Flowchart

Simpson’s Rule

Flowchart