clear
echo on
%Solution to problem 24
%
%For part a. you were asked to derive
%the weights and nodes for the three
%point gaussian quadrature scheme.  This
%rule is of the form:
%
%  gq3 = w(1)*f(x(1)) + w(2)*f(x(2)) + w(3)*f(x(3))
%
%From symmetry, we know that the node x(2)
%will be located at the midpoint x=0, and
%that the other two points will be such that
%x(1) = -x(3).  Also from symmetry, we have that
%w(1) = w(3).
pause

%There remain three undetermined values: w(1), w(2),
%and x(1).  These must be determined from the even
%integrals:
%
%  2 = 2*w(1) + w(2)    (since w(3) = w(1))
%  2/3 = 2*w(1)*x(1)^2  (since x(2) = 0)
%  2/5 = 2*w(1)*x(1)^4
pause

%Dividing the second equation by the third gives us
%that x(1)^2 = 3/5 or that x(1) = -(3/5)^.5
%Plugging this back into the second equation yields
%w(1) = 5/9.  Substituting this into the first equation
%gives us the last value, w(2) = 8/9.  Thus:
xstar=[-(3/5)^.5,0,(3/5)^.5]
wstar=[5/9,8/9,5/9]
pause

%OK, now we use this to compute the integral of sin(x)
%from 0 to pi.  Thus:
a=0;
b=pi;
w=(b-a)/2*wstar;
x=(a+b)/2+(b-a)/2*xstar;
f=sin(x);
igq=w*f';
answer=2;
ergq=igq-answer
pause

%So the error is pretty small.  Now we compare this to
%the three point (two panel) trapezoidal and Simpson's
%rules.  First we set up the evaluation points and weights:
xt=[a,(a+b)/2,b];
ft=sin(xt);
wtr=(b-a)/2*[0.5,1,0.5];
wsimp=(b-a)/6*[1,4,1];
pause

%And now for the quadrature results:
itr=wtr*ft';
isimp=wsimp*ft';
ertr=itr-answer
ersimp=isimp-answer
%Or, comparing the three methods directly,
format short e
[ergq,ertr,ersimp]

%So the gaussian quadrature is much more accurate than
%the other three point rules for this smooth function.
%This is because the gaussian rule is of degree 5 while
%the Simpson's rule is of degree 3 and the trapezoidal
%rule is of degree 1.
echo off