Computational Physics - Department of Physics

5.3 Gaussian Quadrature 127
∫ 100
1

exp(−x)

x

dx,

and ∫ 3

0

1

2 +x^2

dx.

A program example which uses the trapezoidal rule, Simpson’s rule and the Gauss-Legendre
method is included here. For the corresponding Fortran program, replace program1.cpp with
program1.f90. The Python program is listed as program1.py.

http://folk.uio.no/mhjensen/compphys/programs/chapter05/cpp/program1.cpp
#include
#include"lib.h"
using namespacestd;
// Here we define various functions called by the main program
// this function defines the function to integrate
doubleint_function(doublex);
// Main function begins here
intmain()
{
intn;
doublea, b;
cout <<"Read in the number of integration points"<< endl;
cin >> n;
cout <<"Read in integration limits"<< endl;
cin >> a >> b;
// reserve space in memory for vectors containing the mesh points
// weights and function values for the use of the gauss-legendre
// method
doublex =new double[n];
doublew =new double[n];
// set up the mesh points and weights
gauss_legendre(a, b,x,w, n);
// evaluate the integral with the Gauss-Legendre method
// Note that we initialize the sum
doubleint_gauss = 0.;
for(inti = 0; i < n; i++){
int_gauss+=w[i]int_function(x[i]);
}
// final output
cout <<"Trapez-rule = "<< trapezoidal_rule(a, b,n, int_function)
<< endl;
cout <<"Simpson's rule = "<< simpson(a, b,n, int_function)
<< endl;
cout <<"Gaussian quad = "<< int_gauss << endl;
delete[] x;
delete[] w;
return0;
}// end of main program
// this function defines the function to integrate
doubleint_function(doublex)
{
doublevalue = 4./(1.+xx);
returnvalue;
}// end of function to evaluate

To be noted in this program is that we can transfer the name of agiven function to integrate.
In Table 5.2 we show the results for the first integral using various mesh points, while Table
5.3 displays the corresponding results obtained with the second integral. We note here that,
since the area over where we integrate is rather large and theintegrand goes slowly to zero