The following python code allows you to run interactively either in a browser or using ipython notebook. It compares the trapezoidal rule and Gaussian quadrature with the exact result from symbolic python SYMPY up to 1000 integration points for the integral I = 2 = \int_0^{\infty} x^2 \exp{-x} dx. For the trapezoidal rule the results will vary strongly depending on how the infinity limit is approximated. Try to run the code below for different finite approximations to \infty .
xxxxxxxxxxfrom math import expimport numpy as npfrom sympy import Symbol, integrate, exp, oo# function for the trapezoidal ruledef TrapezoidalRule(a,b,f,n): h = (b-a)/float(n) s = 0 x = a for i in range(1,n,1): x = x+h s = s+ f(x) s = 0.5*(f(a)+f(b)) +s return h*s# function for the Gaussian quadrature with Laguerre polynomialsdef GaussLaguerreRule(n): s = 0 xgauleg, wgauleg = np.polynomial.laguerre.laggauss(n) for i in range(1,n,1): s = s+ xgauleg[i]*xgauleg[i]*wgauleg[i] return s# function to computedef function(x): return x*x*exp(-x)# Integration limits for the Trapezoidal rulea = 0.0; b = 10000.0# define x as a symbol to be used by sympyx = Symbol('x')# find result from sympyexact = integrate(function(x), (x, a, oo))# set up the arrays for plotting the relative errorn = np.zeros(40); Trapez = np.zeros(4); LagGauss = np.zeros(4);# find the relative error as function of integration pointsfor i in range(1, 3, 1): npts = 10**i n[i] = npts Trapez[i] = abs((TrapezoidalRule(a,b,function,npts)-exact)/exact) LagGauss[i] = abs((GaussLaguerreRule(npts)-exact)/exact)print "Integration points=", n[1], n[2]print "Trapezoidal relative error=", Trapez[1], Trapez[2]print "LagGuass relative error=", LagGauss[1], LagGauss[2]