Numerical Integration: Gaussian Methods in 1D

Table of Contents

Gauss quadrature Example: Function to integrate Order of integration (number of Gauss points) The final Gauss formula Check the error Generalization to integration over a general interval [a,b] Example: Exercise 1: Exercise 2: Function gaussValues 1D

The Gaussian integration, known also as the method of gaussian quadrature, is a numerical approximation of a definite integral of a function

in a general interval

. We will apply this approach to the computation of the integral of a polynomial function (which, for sufficiently big n, will give an exact result) and of general functions (approximate result).

Gauss quadrature

Gauss quadrature formula is defined for the reference interval

as a weighted sume of the values of

at appropiately chosen points:

We need to choose both points and weights

. The points

are the zeros of Legendre polinomials

. This is a family of polynomial functions represented in the following figure according to the degree (only the first ones). See more information at this Wikipedia link.

For i =1,..,5

The values and weights for these points in the Gauss quadrature formulas are shown in the following table

For example, if we choose n=3 the Gauss approximation formula is

Example:

Now we show an example of the different steps needed to compute the Gauss aprimation value of an integral. As you will see, due to that

is a polynomial function, the result is exact when we choose the proper number of Gauss points.

Function to integrate

As an example, consider the polynomial function

, defined in the reference interval [-1,1]

First we defined the function as an inline function using Matlab

f=@(x) x.^4-3*x.^2+1; %defined as an inline function

Order of integration (number of Gauss points)

When you fix a number of points n, the Gauss integration formula is exact for polynomials of degree 2n-1. In our example the polynomial degree is 4, then if we choose only two points n=2 the integral is still not exact. Therefore, the proper choose is n=3.

n=3; 
% because our function is a polynomial of order 4 we have to choose n=3 
switch (n) %switch executes only one piece of the code according to the value of n
    case 1 %if n=1
        w=2; pG=0;
    case 2 %if n=2
        w=[1,1]; pG=[-1/sqrt(3), 1/sqrt(3)];
    case 3 %if n=3
        w=[5/9, 8/9, 5/9]; pG=[-sqrt(3/5), 0, sqrt(3/5)];
    otherwise
        error('No data are defined for this value');
end
fprintf('n=%d , pG = %f %f %f ',n,pG)
n=3 , pG = -0.774597 0.000000 0.774597 

The final Gauss formula

You can use a sequential formula

sumat = 0;
for i=1:size(pG,2)
    sumat = sumat + w(i)*f(pG(i));
end
intFseq = sumat
intFseq = 0.4000

or equivalently, a more compact form of the same procedure

intFcompact = w*f(pG)'
intFcompact = 0.4000

Check the error

Here, because it is a simple example, we can compute the primitive of the function

and solve the integral using Barrow's rule. Then, we compare the result with the numerical value coming from the Gauss formula.

primitiveF=@(x) x.^5/5-x.^3+x; %just to check the error
barrowRule=primitiveF(1)-primitiveF(-1);
errorInt = abs(barrowRule - intFcompact)
errorInt = 2.2204e-16

Generalization to integration over a general interval [a,b]

When the integration interval is not [-1,1] you have to make a change of variables in the intergral in order to apply the previous formulas to the new interval.

being

and

This means that the Gauss points, initially defined in the interval

, must be transformed into points on the new integration interval

Example:

Compute the integral of the same function, but now in the interval

. So we have to use the change of variables for integrals.

a=2;
b=3;
pGnew = a+(b-a)*(1+pG)/2; %we have already the value of pG
intFcompactNew = (w*f(pGnew)')*(b-a)/2 %the weight values are the same
intFcompactNew = 24.2000

Exercise 1:

According to the table included above, you have to use the function [w,pt] = gaussValues1D(n), included at the end of this script, to see how the approximation of the integral is improved (you must safe the gaussValues1D function in a new file named gaussValues1D.m).

f=@(x) x.^6-x.^4-3*x.^2+1; %defined as an inline function
primitiveF=@(x) x.^7/7-x.^5/5-x.^3+x;  %just to check the error
barrowRuleNew=primitiveF(3)-primitiveF(-2);
a=-2;
b=3;
for n=1:5 %for n=6 must return an error
    [w,pG]=gaussValues1D(n); %function to be implemented
    pGnew = a+(b-a)*(1+pG)/2; %we have already the value of pG
    intFcompactNew = (w*f(pGnew)')*(b-a)/2; %the weight values are the same
    errorInt = abs(barrowRuleNew - intFcompactNew);
    fprintf('for n= %d the error is =%e \n',n,errorInt)
end
for n= 1 the error is =2.446987e+02 
for n= 2 the error is =1.769180e+02 
for n= 3 the error is =2.790179e+01 
for n= 4 the error is =2.842171e-13 
for n= 5 the error is =4.263256e-13 

Exercise 2:

Using the function gaussValues1D, approximate the value of the integral in

of the function

. Use different values of n and show the errors comparing with the true integral value

Sol:

n, error

1 3.1706e-01

2 7.1183e-03

3 6.1578e-05

4 2.8092e-07

5 7.9140e-10

Function gaussValues 1D

function [w,pt]=gaussValues1D(n)    
    switch (n)
        case 1
            w=2; pt=0;
        case 2
            w=[1,1]; pt=[-1/sqrt(3), 1/sqrt(3)];
        case 3
            w=[5/9, 8/9, 5/9]; pt=[-sqrt(3/5), 0, sqrt(3/5)];
        case 4
            w=[(18+sqrt(30))/36, (18+sqrt(30))/36, (18-sqrt(30))/36, (18-sqrt(30))/36]; 
            pt=[-sqrt(3/7-2/7*(sqrt(6/5))), sqrt(3/7-2/7*(sqrt(6/5))),...
                -sqrt(3/7+2/7*(sqrt(6/5))), sqrt(3/7+2/7*(sqrt(6/5)))];
        case 5
            w=[(322+13*sqrt(70))/900, (322+13*sqrt(70))/900, 128/225,...
                (322-13*sqrt(70))/900,(322-13*sqrt(70))/900]; 
            pt=[-1/3*sqrt(5-2*sqrt(10/7)), 1/3*sqrt(5-2*sqrt(10/7)),0,...
                -1/3*sqrt(5+2*sqrt(10/7)), 1/3*sqrt(5+2*sqrt(10/7))];
        otherwise
            error('No data are defined for this value');
    end
end