Table of Contents

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 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

.

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.

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

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)

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

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

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

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

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

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 [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

(c) Numerical Factory 2023