clear
echo on
format short e
format compact
%This program illustrates the accuracy of Gaussian
%quadrature routines.  It calculates the integral
%of exp(x) over the domain [-1,1] using quadrature
%rules with n = 1 to 5.  The resulting error is
%compared to a quadrature rule where the weights
%are optimized, but the nodes are evenly spaced.

%First for the Gaussian Quadrature

%We begin with the midpoint rule (n=1)
ngq(1)=1;
wgq1=2;
xgq1=0;
igq(1)=wgq1*exp(xgq1);
ergq(1)=igq(1)-(exp(1)-exp(-1))
pause

%Now the two-point rule
ngq(2)=2;
wgq2=[1,1];
xgq2=[-1/3^.5,1/3^.5];
igq(2)=wgq2*exp(xgq2)';
ergq(2)=igq(2)-(exp(1)-exp(-1))
pause

%And the three point rule
ngq(3)=3;
wgq3=[5/9,8/9,5/9];
xgq3=[-0.6^.5,0,0.6^.5];
igq(3)=wgq3*exp(xgq3)';
ergq(3)=igq(3)-(exp(1)-exp(-1))
pause

%And the four point rule
ngq(4)=4;
p=(30+480^.5)/70;
q=(30-480^.5)/70;
r=(1-3*q)/3/(p-q);
wgq4=[r,1-r,1-r,r]
xgq4=[-p^.5,-q^.5,q^.5,p^.5]
igq(4)=wgq4*exp(xgq4)';
ergq(4)=igq(4)-(exp(1)-exp(-1))
pause

%And finally the five point rule:
ngq(5)=5;
p=(70+1120^.5)/126;
q=(70-1120^.5)/126;
r1=(3-5*q)/15/p/(p-q);
r2=(5*p-3)/15/q/(p-q);
wgq5=[r1,r2,128/225,r2,r1]
xgq5=[-p^.5,-q^.5,0,q^.5,p^.5]
igq(5)=wgq5*exp(xgq5)';
ergq(5)=igq(5)-(exp(1)-exp(-1))
pause

%We may plot this up:
semilogy(ngq,abs(ergq))
xlabel('number of points')
ylabel('error')
%Note that the error falls off faster
%than exponentially in the number of
%points!
pause

%OK, now for non-optimally placed points.  In this
%case we will start with n=1, which is the same
%as the Gaussian Quadrature midpoint rule:
n(1)=1;
i(1)=2*exp(0);
er(1)=i(1)-(exp(1)-exp(-1))
pause

%Now we go to the two point rule with evenly
%spaced nodes, the trapezoidal rule:
n(2)=2;
x2=[-1,1];
w2=[1,1];
i(2)=w2*exp(x2)';
er(2)=i(2)-(exp(1)-exp(-1))
%It is interesting to note that the midpoint
%rule is actually -more- accurate than the two
%point Trapezoidal rule!
pause

%Now for n=3 (this is Simpson's rule)
n(3)=3;
x3=[-1,0,1];
w3=[1/3,4/3,1/3];
i(3)=w3*exp(x3)';
er(3)=i(3)-(exp(1)-exp(-1))
pause

%and for n=4
n(4)=4;
x4=[-1,-1/3,1/3,1];
w4=[1/4,3/4,3/4,1/4];
i(4)=w4*exp(x4)';
er(4)=i(4)-(exp(1)-exp(-1))
pause

%Finally, for n=5
n(5)=5;
x5=[-1,-.5,0,.5,1];
w5=[7/45,32/45,12/45,32/45,7/45];
i(5)=w5*exp(x5)';
er(5)=i(5)-(exp(1)-exp(-1))
pause

%And then we can plot this up together
%with the Gaussian quadrature result:
semilogy(ngq,abs(ergq),n,abs(er))
xlabel('number of points')
ylabel('error')
%As you can see, the optimally spaced nodes of
%Gaussian quadrature gives a much more accurate
%calculation of the integral than do non-optimally
%spaced nodes.
echo off