Subject : The following demonstrates the convergence of the Gaussian method of estimating integrals of continuous functions.
Authors : Nathan Collier, Autar Kaw, Loubna Guennoun
Date : 27 October 2005

Introduction

auss Quadrature Rule is another method of estimating an integral. The theory behind the two point Gauss Quadrature Rule is to approximate the integral by taking the area  under a straight line connecting any two points on the curve that are not predetermined as a and b, but as unknowns x1 and x2. For n-points rules, the general form to approximate the integral is

where ci and xi are the weighting factors and function arguments used in Gauss Quadrature formulas, respectively. However, these factors and arguments are already defined to approximate any integral from -1 to 1. To be able to use them, the limits of the integral of the function f(x) need to be changed to [-1,1].

NOTE: Weighting factors c and function arguments x used in Gauss Quadrature Rule have already been defined  in the textbook for  up to six points.
The following procedure will illustrate the Gauss Quadrature Rule of integration. The user may enter any function f(x), the lower and upper limit for the function, and the number of points n in the data section (up to six points). By entering this data, the program will calculate the exact value of the integral, followed by the results using the Gauss Quadrature Rule with n points. The program will also display the true error, the absolute relative true percentage error, the approximate error, the absolute relative approximate percentage error, and the number of significant digits that are at least correct.

Inputs

Integrand f[x] = 0

f[x_] := 300.0*x/(1.0 + Exp[x]);

Lower and upper limit of the integral, a and b respectively.

a = 0.0;
b = 10.0;

Maximum number of points. This number must be between 1 and 6.

(Subscript[n, maximum]) = 6;

Weighting Factors and Functions Arguments for up to 6 Points

Weighting Factors for 1-6 Points

(Subscript[C, 1, 0]) = 2.0;
(Subscript[C, 1, 1]) = 1.0;
(Subscript[C, 2, 1]) = 1.0;
(Subscript[C, 1, 2]) = 0.5555556;
(Subscript[C, 2, 2]) = 0.8888889;
(Subscript[C, 3, 2]) = 0.5555556;
(Subscript[C, 1, 3]) = 0.347854845;
(Subscript[C, 2, 3]) = 0.652145155;
(Subscript[C, 3, 3]) = 0.652145155;
(Subscript[C, 4, 3]) = 0.347854845;
(Subscript[C, 1, 4]) = 0.236926885;
(Subscript[C, 2, 4]) = 0.478628670;
(Subscript[C, 3, 4]) = 0.568888889;
(Subscript[C, 4, 4]) = 0.478628670;
(Subscript[C, 5, 4]) = 0.236926885;
(Subscript[C, 1, 5]) = 0.171324492;
(Subscript[C, 2, 5]) = 0.360761573;
(Subscript[C, 3, 5]) = 0.467913935;
(Subscript[C, 4, 5]) = 0.467913935;
(Subscript[C, 5, 5]) = 0.360761573;
(Subscript[C, 6, 5]) = 0.171324492;

Function Arguements for 1-6 Points

(Subscript[X, 1, 0]) = 0.0;
(Subscript[X, 1, 1]) = -0.5773503;
(Subscript[X, 2, 1]) = 0.5773503;
(Subscript[X, 1, 2]) = -0.774596669;
(Subscript[X, 2, 2]) = 0.0;
(Subscript[X, 3, 2]) = 0.774596669;
(Subscript[X, 1, 3]) = -0.861136312;
(Subscript[X, 2, 3]) = -0.339981044;
(Subscript[X, 3, 3]) = 0.339981044;
(Subscript[X, 4, 3]) = 0.861136312;
(Subscript[X, 1, 4]) = -0.906179846;
(Subscript[X, 2, 4]) = -0.538469310;
(Subscript[X, 3, 4]) = 0.0;
(Subscript[X, 4, 4]) = 0.538469310;
(Subscript[X, 5, 4]) = 0.906179846;
(Subscript[X, 1, 5]) = -0.932469514;
(Subscript[X, 2, 5]) = -0.661209386;
(Subscript[X, 3, 5]) = -0.238619186;
(Subscript[X, 4, 5]) = 0.238619186;
(Subscript[X, 5, 5]) = 0.661209386;
(Subscript[X, 6, 5]) = 0.932469514;

curve = 
        Plot[f[x], {x, a, b}, 
          PlotLabel -> "Entered function on given interval", 
          TextStyle -> {FontSize -> 
                11}];
General::spelll  Off;

True Solution

This is the solution found by Mathematica

Actual = Integrate[f[x], {x, a, b}]

Value of integral as a function of iterations

The integral given above has the limits of[a,b].It needs to be converted into an integral with limits[-1,1] fnew(x) is the new function that will be used for evaluating the integral using the Gauss Quadrature rule

fnew[x_] := 
        f[(b - a)*x/2 + (b + a)/2]*(b - a)/2;

Gaussian

Array[AV, (Subscript[n, maximum])];

For[i = 1,  
      i ≤ (Subscript[n, maximum]), i++, 
AV[i] = 
        Sum[(Subscript[C, j, i - 1])*fnew[(Subscript[X, j, i - 1])], {j, 1, i}]]

Absolute true error

Array[Et, (Subscript[n, maximum])];

For[i = 1, i <= (Subscript[n, maximum]), i++, 
      Et[i] = Abs[Actual - AV[i]]]

Absolute relative true error

Array[ϵt, (Subscript[n, maximum])];

For[i = 1, 
      i <= (Subscript[n, maximum]), i++, ϵt[i] = 
        Abs[Et[i]/Actual*100]]

Absolute approximate error

Array[Ea, (Subscript[n, maximum])];

For[i = 1, i <= (Subscript[n, maximum]), i++, 
      If[i <= 1,  Ea[i] =  0, Ea[i] = Abs[AV[i] - AV[i - 1]]]]

Absolute relative approximate error

Array[ϵa, (Subscript[n, maximum])];

For[i = 1, i <= (Subscript[n, maximum]), i++, 
      If[i <= 1, ϵa[i] = 0, ϵa[i] = 
          Abs[Ea[i]/AV[i]*100]]]

Significant digits at least correct

Array[sigdig, (Subscript[n, maximum])];

For[i = 1, i <= (Subscript[n, maximum]), i++, 
      If[(ϵa[i] ≥ 5)  ||  (i <= 1), 
        sigdig[i] = 0, 
        sigdig[i] = 
          Floor[(2 - Log[10, Abs[ϵa[i]/0.5]])]]]

Graphs

xrplot = Table[AV[i], {i, 1, (Subscript[n, maximum])}];

ListPlot[xrplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[xrplot]}, PlotLabel -> "Estimated root as a function of number \
of iterations"];

Etplot = Table[Et[i], {i, 1, (Subscript[n, maximum])}];

ListPlot[Etplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[Etplot]}, PlotLabel -> "Absolute true error as a function of \
number of iterations"];

ϵtplot = 
        Table[ϵt[i], {i, 1, (Subscript[n, maximum])}];

ListPlot[ϵtplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[ϵtplot]}, PlotLabel -> "Absolute relative true error as a \
function of number of iterations"];

Eaplot = Table[Ea[i], {i, 1, (Subscript[n, maximum])}];

ListPlot[Eaplot,PlotJoined -> True, PlotRange -> All,AxesOrigin -> {1, Min[Eaplot]}, PlotLabel -> "Absolute approximate error as a function \
of number of iterations"];

ϵaplot = 
        Table[ϵa[i], {i, 1, (Subscript[n, maximum])}];

ListPlot[ϵaplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[ϵaplot]}, PlotLabel -> "Absolute relative approximate error as a \
function of number of iterations"];

sigdigplot = Table[sigdig[i], {i, 1, (Subscript[n, maximum])}];

<< Graphics`Graphics`

BarChart[sigdigplot];