clear
echo on
format compact
format short e
%This program illustrates the use of mapping to
%apply gaussian quadrature to an arbitrary domain,
%and the utility of removing a singularity from
%an integrand using analytic techniques.  we will
%use the function cos(x)log(x) as an example.  We
%shall store the function under the name 'f1.m'.
%This function has an integrable singularity at
%x=0.  We want to integrate the function over the
%domain [0,1].  First, let's plot the integrand.

y=[.0001:.01:1];
plot(y,f1(y))
xlabel('x')
ylabel('cos(x)log(x)')
%Note that the integrand has a nasty functional
%dependence around x=0.
pause

%Let's see what this does to the accuracy of the
%gaussian quadrature scheme.  We first calculate
%the exact integral analytically (it turns out that
%we can get it as an infinite series).  The integral
%is -.946083.  So:
answer=-9.460830703672338e-01;

%Let's compare this to the quadrature rules for 1
%to 4 points.  We shall do this using the attached
%function 'gauss.m' which returns the four integrals.
%Thus:
a=0;
b=1;
er1=gauss('f1',a,b)-answer
semilogy(abs(er1))
xlabel('number of points in quadrature')
ylabel('error')
pause

%As you can see, the error decreases only very slowly
%as the number of points increases.  What happens if
%we instead remove the singularity analytically?  In
%this case we simply subtract off the singularity
%and integrate it analytically.  Thus we solve the
%problem: log(x) + (cos(x)-1)log(x).  The integral
%of log(x) over the domain is just -1.  First we plot
%the new integrand (cos(x)-1)log(x).  We have saved
%this function under the name 'f2.m':

plot(y,f2(y))
xlabel('x')
ylabel('(cos(x)-1)log(x)')
%This certainly is more well behaved near zero!
pause

%Now for the error of the gaussian quadrature routine
%on this new integrand:
er2=gauss('f2',a,b)-1-answer
%Note that we had to subtract off the analytic result
%for the integral of log(x).
n=[1:4];
semilogy(n,abs(er1),n,abs(er2))
xlabel('number of points in quadrature')
ylabel('error')
pause

%As you can see, this is a big improvement.  We can
%still do better, however, because while the function
%(cos(x)-1)log(x) is well behaved at the origin, its
%higher derivatives are still singular!  In particular,
%its second derivative is singular.  This is why the
%4 point quadrature error bends up rather than down.
%Just continuing on with higher order quadrature rules
%(increasing the number of points uniformly) would yield
%little benefit.  We -could- get an accurate result by
%using an adaptive algorithm on this function and increasing
%the number of points in the region of the singularity,
%but there is another way to do this!
pause
%It turns out that the integral of the original integrand
%can be simplified by integration by parts.  Doing this,
%we find that the integral of cos(x)log(x) over [0,1] is
%the same as -sin(x)/x.  This integrand is well behaved 
%everywhere.  Let's try integrating this one.  We save it
%under the name 'f3.m':
pause
%First we look at the integrand again:
plot(y,f3(y))
xlabel('x')
ylabel('-sin(x)/x')
pause
%And now for the integration error:
er3=gauss('f3',a,b)-answer
semilogy(n,abs(er1),n,abs(er2),n,abs(er3))
xlabel('number of points in quadrature')
ylabel('error')
pause

%We now have even higher accuracy!  The error in the four
%point rule is on the order of 5e-11 for the last integrand,
%as compared to 3e-2 for the original integrand and 2e-5 if
%we removed the singularity analytically. The moral of
%the story is that if you are having trouble integrating
%a function, either subtract off the asymptotic form
%that is singular, or transform the integrand into a
%form that is better behaved.
pause
%As a last example, let's see what the adaptive algorithm
%'quad' does with this function.  We can't even apply it to
%the first function because quad uses the Simpson's rule, 
%which requires evaluation at the end points.  Thus we shall
%integrate the second function instead.  We shall call it
%f2plot because we have added a bit of graphics to it.  We
%also have to prevent it from crashing at x=0.  Thus:
pause
plot(y,f2(y))
xlabel('x')
ylabel('(cos(x)-1)log(x)')
hold on
error=quad('f2plot',a,b)-1-answer
hold off
pause
%As you can see, the singular second derivative at x=0 really
%forced a large number of evaluations!  The accuracy of the
%end result isn't much better than the 4 point Gaussian rule
%either.

echo off


_______________________


function igq=gauss(f,a,b)
%This function calculates the integral of the function
%specified by the string variable 'f' for the 1, 2, 3,
%and 4 point gaussian quadrature rules.  The four values
%are returned as an array.

%We define the midpoint:
mid=(b+a)/2;

%We begin with the midpoint rule (n=1)
xstar=0;
wstar=2;
x=mid+xstar*(b-a)/2;
w=wstar*(b-a)/2;
igq(1)=sum(w.*feval(f,x));

%Now the two-point rule
xstar=[-1/3^.5,1/3^.5];
wstar=[1,1];
x=mid+xstar*(b-a)/2;
w=wstar*(b-a)/2;
igq(2)=sum(w.*feval(f,x));

%And the three point rule
xstar=[-0.6^.5,0,0.6^.5];
wstar=[5/9,8/9,5/9];
x=mid+xstar*(b-a)/2;
w=wstar*(b-a)/2;
igq(3)=sum(w.*feval(f,x));

%Finally the four point rule
xstar=[-0.86113631159405,-0.33998104358486,0.33998104358486,0.86113631159405];
wstar=[0.34785484513745,0.65214515486255,0.65214515486255,0.34785484513745];
x=mid+xstar*(b-a)/2;
w=wstar*(b-a)/2;
igq(4)=sum(w.*feval(f,x));


_____________________


function y=f1(x)
y=cos(x).*log(x);


_____________________


function y=f2(x)
y=(cos(x)-1).*log(x);


_____________________


function y=f3(x)
y=-sin(x)./x;


_____________________


function y=f2plot(x)
n=length(x);
for i=1:n,
 if x(i)==0
  y(i)=0;
 else
  y(i)=(cos(x(i))-1).*log(x(i));
 end
end
plot(x,y,'o')
drawnow