Python, NumPy and MatPlotLib Visualising Numerical Integration

Perquisites

Python

Finding the Area Under a Curve

Let’s visually inspect, how to calculate the area under the curve using a line plot and a bar plot. This forms the basis of integration. Recall that integration can be related to a physical quantity, take the every day example of a car. The function is the speed of a car and the area under the curve is the distance the car has covered.

Constant Term

Lets have a look at the x,y plot of a function that is a constant, this means that the speed is constant and if a car travels at a constant speed, then its distance should just be proportional to the time (x axis) the car has been driving. We will look at x beginning at 0 and ending at 10 (zero order indexing). We will also create a bar plot of x,y.

\displaystyle \text{y}=1

Python

We can see that area under the curve is the length of the vector minus the 2 ends times the constant height of 1. The line only covers half the 0th and end point, so it is actually the length of the vector minus 1 times the constant height. We can compare this with the function trapz() which gives the same result.

Python
lenx=
10
area=
9.0
int1=
9.0

Physically this makes sense as for constant speed, distance equals speed times time.

Linear Term

Lets have a look at the x,y plot of a function where y is equal to x. We will look at x beginning at 0 and ending at 10 (zero order indexing). Once again we will also create a bar plot of x,y. Continuing the analogy of the car, the car speeds up proportionally to the time.

\displaystyle \text{y}=\text{x}

Python

Once again we can calculate the area under the line. As by definition y is equal to x, both lengths are the same. Therefore one can easily see if this plot was square it would be the length of the vector minus 1 squared. However we see that it isn't square, that it is in fact more like a triangle. We can compare this with the function trapz() which gives the same result.

Python
lenx=
10
area2ifsquare=
81
area2iftriangle=
40.5
int2=
40.5

Constant

\displaystyle \text{area}=\frac{1}{\mathbf{1}}{{\left( {\text{lenx}-1} \right)}^{\mathbf{1}}}

Linear Term

\displaystyle \text{area}=\frac{1}{\mathbf{2}}{{\left( {\text{lenx}-1} \right)}^{\mathbf{2}}}

Physically this makes sense also. If the car starts off very slowly, it will cover a small amount of distance and as it speeds up the amount of distance it covers per unit of time will be much larger. When the car's speed is directly proportional to the time, we will see a triangle plot such as the above.

Squared Term

Now lets have a look at the x,y plot of a function where y is equal to x squared. We will look at x beginning at 0 and ending at 10 (zero order indexing). Once again we will also create a bar plot of x,y.

\displaystyle \text{y}={{\text{x}}^{2}}

Python

Once again we can calculate the area under the line. By definition y is equal to x squared and one can easily see if this plot was a rectangle it would be the length of the vector minus 1 cubed. However we see that it isn't square and it is less than a triangle. In fact it covers 1/3 of the area of the rectangle.

Python
lenx=
10
area3ifrectangle=
729
area3iftriangle=
364.5
area3=
243.0
int3=
244.5

Here we can once again see that the trapz function used for integration gives very similar results. There is a slight difference which will be discussed towards the end.

Constant

\displaystyle \text{area}=\frac{1}{\mathbf{1}}{{\left( {\text{lenx}-1} \right)}^{\mathbf{1}}}

Linear Term

\displaystyle \text{area}=\frac{1}{\mathbf{2}}{{\left( {\text{lenx}-1} \right)}^{\mathbf{2}}}

Squared Term

\displaystyle \text{area}=\frac{1}{\mathbf{3}}{{\left( {\text{lenx}-1} \right)}^{\mathbf{3}}}

Cubic Terms

Lets have a look at the x,y plot of a function where y is equal to x cubed. We will look at x beginning at 0 and ending at 10 (zero order indexing). Once again we will also create a bar plot of x,y.

\displaystyle \text{y}={{\text{x}}^{3}}

Python

Once again we can calculate the area under the line. By definition y is equal to x cubed and one can easily see if this plot was a rectangle it would be the length of the vector minus 1 to the power of four. However we see that it isn't square and it is less than a triangle, it is less than 1/3 of the rectangle, in fact it covers 1/4 of the area of the rectangle.

Python
lenx=
10
area4=
1640.25
int4=
1660.5

Here we can once again see that the trapz function used for integration gives very similar results but there is a substantial difference which will once again be discussed towards the end.

Constant

\displaystyle \text{area}=\frac{1}{\mathbf{1}}{{\left( {\text{lenx}-1} \right)}^{\mathbf{1}}}

Linear Term

\displaystyle \text{area}=\frac{1}{\mathbf{2}}{{\left( {\text{lenx}-1} \right)}^{\mathbf{2}}}

Squared Term

\displaystyle \text{area}=\frac{1}{\mathbf{3}}{{\left( {\text{lenx}-1} \right)}^{\mathbf{3}}}

Cubed Term

\displaystyle \text{area}=\frac{1}{\mathbf{4}}{{\left( {\text{lenx}-1} \right)}^{\mathbf{4}}}

Multiplication by a Scalar

Let's now look at the effect of multiplication by a scalar, in this case 2.

\displaystyle \text{y}=2\text{x}

Python
Python
lenx=
10
area2iftriangle=
81.0
int5=
81.0

As we can see multiplication of the function by a scalar simply just multiplies the area by the same scalar:

\displaystyle \text{y}=\mathbf{2}\text{x}

\displaystyle \text{area}=\mathbf{2}\times \frac{1}{2}{{\text{x}}^{2}}

Addition of Terms

Let's now have a look at the addition of a squared term and a cubic term.

Python

Here we can see that the area under the curve of a function consisting of an x squared and x cubed term is the sum of the individual areas of these terms.

Python
int6a=
244.5
int6b=
1660.5
int6c=
1905.0
int6a+int6b=
1905.0

Uncertainty in Results

As we seen looking at the squared term there is a small error when measuring the area under the curve, this is higher for the cubic term. The reason for this is the bars are not completely below or above the curve. Due to the thickness of each bar, some of the bars are above the red line and there is some empty white space below the red line. This is in essence the source of the error. We can modify the code for the squared term by adding a bar width (line 1, 2, 20 and 23)

Python
lenx=
10
area3=
243.0
int3=
244.5

Now we can make the barwidth smaller going from a size of 1.0 to 0.5 (line 1).

lenx=
20
area3=
285.79166666666663
int3=
286.1875

Note the values are closer to each other now that the bars are half height. The reason both values are higher is due to 0 order indexing. We first went from 0 to 10 which gave use the values [0,1,2,3,4,5,6,7,8,9], with 20 steps we are instead going from from 0 to 20 which gives [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19], this additional 19th point at the end of the curve makes the value slightly higher (in essence we closer approach the upper limit of 10). As can be seen, the space occupied by bars above the curve and the space not occupied by bars below the curve is smaller. To rectify this inconstency, the code will be modified to stop at 9+barwidth (line 9).

Python
lenx=
19
area3=
243.0
int3=
243.375

Lets now try to make the barwidth 0.1 wide (line 1).

lenx=
91
area3=
243.0
int3=
243.01500000000001

These two values become much more similar and the space occupied by bars above the curve and the space not occupied by bars below the curve is much smaller. Let's now try a barwidth of 0.01, now each individual bar cannot be distinguished.

lenx=
901
area3=
243.0
int3=
243.00015

The area3 and int3 now match up to 3 decimal places. Integration in theory uses an infinitely small width of bar but the smaller the width of bar, the longer it will take a computer to make the calculation. Finally computers have a precision limit and can only make bars so small however the precision on a computer is generally sufficient for numerical integration of most problems.

We can repeat this for the cubed term with a barwidth of 1, 0.5, 0.1 and 0.01:

Python
lenx=
10
area4=
1640.25
int4=
1660.5
lenx=
19
area4=
1640.25
int4=
1645.3125
lenx=
91
area4=
1640.25
int4=
1640.4525000000003
lenx=
901
area4=
1640.25
int4=
1640.2520250000002