To find a numerical value for the integral of a real valued function of a real variable over a specific range over the real line. This means to evaluate
Geometrically, this integral represents the area under from to
We can always approximate the area by dividing it in equal width strips and then sum the areas of all the strips.
In general, there will always be an error in the estimate of the area using this method. The error will become smaller the more strips we use (which implies a smaller strip width). Hence we can write
Where is the error between the actual area and the approximated area using the above method of numerical integration. above is the number of strips or can also be refereed to as the number of integration points.
Instead of keep referring to the ’width of the strip’ all the time, we will call this quantity the weight that we will multiply the value of the function with to obtain the area. Hence the above becomes
Using implied summation on indices the above becomes
In the above we divided the range of the integration (the distance between the upper and lower limits of integration) into equal intervals. We can decide to evaluate at the middle of the strip or at the start of the strip or at the end of the strip. In the diagram above we evaluated the at the left end of the strip.
Our goal is to evaluate this integral such as the error is minimum and using the smallest number of integration points. In a sense this can be considered an optimization problem with constraints: minimize the error of integration using the smallest possible number of points.
To be able to do this minimization, we need to consider what are the variables involved. We see that there are two degrees of freedom to this problem. One is the width of the strip or . We do not have to use a fixed value of width, we can use different width for different strips if the resulting integral gives better approximation.
The second degree of freedom is the point at which to evaluate associated with each strip. In the example above we choose to evaluate at left end point of the strip. We can choose to select a different point if this will result in a better approximation.
This is the main idea of Gauss Quadrature numerical integration. It is to be able to choose specific values for these two degrees of freedom, the and the .
It turns out that if the function is a polynomial, then there is an optimal solution. There is an optimal for each polynomial of order .
We can determine these degrees of freedom such that the error is zero, and with the least possible number of integration points. We are able to tabulate these two degrees of freedom for each polynomial of specific order. In other words, if the function is a polynomial of order then we know before the computation starts what these 2 degrees of freedom should be. We know the locations of and we know weight that we need to multiply with to obtain the area with minumum error.
You might ask how can this method of integration know the locations of the integration points beforehand without being given the integration range of the function to integrate? It turns out that we will map into a new known and specific range of integration (from to ) for the method that we will now discuss.
From now on we will assume the function to be integrated is a polynomial in of some order .
Gauss quadrature is optimal when the function is a polynomial
The main starting point is to represent the function as a combination of linearly independent basis.
Instead of using strips of equal width, we assume the width can vary from one strip to the next. Let us call the width of the strip . Instead of taking the height of the strip to be the value of the function at the left edge of the strip, let us also be flexible on the location of the associated with strip and call the height of the strip as where is to be determined. Hence the above integration becomes
So our goal is to determine and such as the error is minimized in the above equation. We would really like to find and such that the error is zero with the smallest value for .
It seems as if this is a very hard problem. We have unknown quantities to determine. different widths, and associated points to evaluate the height of each specific strip at. And we only have as an input and the limits of integration, and we need to determine these quantities such that the error in integration is zero.
In other words, the problem is to find such that
One way to make some progress is to expand as a series. We can approximate as convergent power series for example. If happens to be a polynomial instead, we can represent it exactly using a finite sequence of Legendre polynomials. It is in this second case where this method makes the most sense to use due to the advantages we make from the second representation.
We show both methods below.
Expanding as convergent power series over the range gives
Substituting (2) into (1) gives
Substituting the above into (3) results in
But from (2) we see that
Substituting the above into (4) gives
Equating coefficients on both sides results in
Since we have unknown quantities to solve for ( weights and points ) we need equations. In other words, we need to have . The above set of simultaneous equations can now be solved for the unknown and this will give us the numerical integration we wanted.
The above assumed that the function can be represented exactly by the power series expansion with terms.
We now consider the representation of as a series of Legendre polynomials. We do this since when itself is a polynomial. We can represent exactly by a finite number of Legendre polynomials. But since Legendre polynomials are defined over we start by transforming to this new range and then we can expand the mapped (which we will call ) in terms of the Legendre polynomials.
An easy way to find this mapping is to align the ranges over each others and take the ratio between as the scale factor. This diagram shows this for a general case where we map defined over to a new range defined over
We see from the diagram that
is the same ratio as
The above is called the Jacobin of the transformation. Now, From the diagram we see that
Hence (6) becomes
Hence we obtain that
For the specific case when and the above expressions become
Before we proceed further, It will be interesting to see the effect of this transformation on the shape of some functions. Below I plotted some functions under this transformation. The left plots are the original functions plotted over some range, in this case and the right side plots show the new shape (the function ) over the new range
We must remember that in the following analysis, we are integrating now the function over and not over . Hence to obtain the required integral we need to transform back the final result. We will show how to do this below.
We can approximate any function as a series expansion in terms of weighted sums of a set of basis functions. We do this for example when we use Fourier series expansion.
Hence we write
We can give an intuitive justification to the above formulation as follows. If we think in terms of vectors. A vector in an n-dimensional space is written in terms of its components as follows
Where is the basis vectors in that space.
If we now consider a general infinite dimensional vector space where each point in that space is a function, then we see that we can also represent that function as a weighted series of a basis functions just as we did for a normal vector.
There are many sets of basis functions we can choose to represent the function with. The main requirements for the basis functions is that they are complete (This means they span the whole space) and there is defined an inner product on them.
For our purposes here, we are interested in the class of function that are polynomials in . The basis we will select to use are the Legendre basis as explained above. To do this, we transform to as shown above and then now our integral becomes
This is because we found that from above.
If we call as the integral becomes
Since is the Jacobian of the transformation, we write the integral as
Since the Jacobian is constant in this case, we will only worry about and we just need to remember to scale the result at the end by . This is the integral we will now numerically integrate.
Equation (7) is now written as
Where is the Legendre polynomial of order and is a polynomial in or some order .
Now we carry the same analysis we did when we expressed as a power series. The difference now is that the limit of integration is symmetrical and the basis are now the Legendre polynomials instead of the polynomials as the case was in the power series expansion. So now equation (1) becomes
And equation (2) becomes
Substituting (9) into (8) we get the equivalent of equation (3)
The above occurs since the integral of any Legendre polynomial of order greater than zero will be zero over
Substituting the above into (10) we obtain
But from (9) we see that
Substituting the above in (11) gives
Rearranging results in
Equating coefficients gives
If we select the points to be the roots of we can write the above as