1D Gaussian Quadratures

Gaussian integration (quadrature) methods are numerical approximation of definite integral in a general interval [a,b].

We will apply this to the computation of the integral of a polynomial funcion (exact) and general functions (approximate).

$$\int_{a}^{b}{f(x) dx}= \sum \omega_i f(p_i)$

Gauss quadrature
Function to integrate
Order of integration (number of Gauss points)
The final formula
Exercise 1:
Exercise 2:

Gauss quadrature

Gauss quadrature is defined for the reference interval [-1,1] and the choosed points are the zeros of Legendre polinomials $$P_i(x)$ represented in the following figure

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

Function to integrate

Let's consider the function

$$ f(x)=x^4-3x^2+1$

Defined in the reference interval [-1,1]

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

or it can be also defined as a file function

function y=f(x)

y=x.^4-3*x.^2+1;

Order of integration (number of Gauss points)

As you increment n the formula is exact for polynomials of degree 2n-1.

n=3;
% because our function is a ploynomial of order 4 we have to choose n=3
switch (n)
    case 1
        w=2; pG=0;
    case 2
        w=[1,1]; pG=[-1/sqrt(3), 1/sqrt(3)];
    case 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

The final 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 =

   4.0000e-01

or equivalently, a compact form of the same value

intFcompact = sum(w.*f(pG))

intFcompact =

   4.0000e-01

Check the error against the actual value

primitiveF=@(x) x.^5/5-x.^3+x;
barrowRule=primitiveF(1)-primitiveF(-1);
errorInt = abs(barrowRule - intFcompact)

errorInt =

   2.2204e-16

Exercise 1:

make the assignment of the Gauss points and weight values a Matlab function for n=1,..5 (check n <=5, bigger values are not allow, use the following sentence).

error('No data are defined for this value');

According to the table included above you have to implement and get the results:

function [w,pt] = gaussValues1D(n)

for n=1:6 %for n=6 must return an error
    n
    [w,pt]=gaussValues1D(n) %function to be implemented
end

n =

     1


w =

     2


pt =

     0


n =

     2


w =

     1     1


pt =

  -5.7735e-01   5.7735e-01


n =

     3


w =

   5.5556e-01   8.8889e-01   5.5556e-01


pt =

  -7.7460e-01            0   7.7460e-01


n =

     4


w =

   6.5215e-01   6.5215e-01   3.4785e-01   3.4785e-01


pt =

  -3.3998e-01   3.3998e-01  -8.6114e-01   8.6114e-01


n =

     5


w =

   4.7863e-01   4.7863e-01   5.6889e-01   2.3693e-01   2.3693e-01


pt =

  -5.3847e-01   5.3847e-01            0  -9.0618e-01   9.0618e-01


n =

     6

Error using gaussValues1D (line 19)
No data are defined for this value

Error in Gauss1D (line 83)
    [w,pt]=gaussValues1D(n) %function to be implemented

Exercise 2:

Using the function gaussValues1D, approximate the value of the integral in [-1,1] of the function f(x)=cos(x). Use different values of n and show the errors comparing with the true value 2*sin(1)

Sol:

n, error

1   3.1706e-01

2   7.1183e-03

3   6.1578e-05

4   2.8092e-07

5   7.9140e-10

(c) Numerical Factory 2016