clear
echo on
%In this example we demonstrate the
%use of the matlab routine quad8.m to
%perform a 2-D integral.  We shall
%integrate the function 1/exp(x^2+y^2)
%over the elliptical disk described by:

%  (x/a)^2 + (y/b)^2 < 1

%where a and b are the major and minor
%axes of the bounding ellipse, respectively.
pause
%First, let's plot up this function over its
%domain.  We shall take a to be 2 and b to be
%unity:
global a b
a=2;
b=1;
x=[-a:.05:a];
y=[-b:.05:b];
z=zeros(length(y),length(x));
for j=1:length(x)
 for i=1:length(y)
  if (x(j)^2/a^2+y(i)^2/b^2)<1
   z(i,j)=1/exp(x(j)^2+y(i)^2);
  end
  echo off
 end
end
echo on
figure(1)
mesh(x,y,z)
xlabel('x')
ylabel('y')
pause
%As you can see from the plot, the function is
%a maximum in the center and decreases toward
%the edges of the ellipse.  This is clearer
%in the contour plot:
figure(2)
contour(x,y,z)
xlabel('x')
ylabel('y')
hold on
plot(x,(1-x.^2/a^2).^.5,'w')
plot(x,-(1-x.^2/a^2).^.5,'w')
%The roughness around the edges is due to the
%finite discretization of the bounding ellipse.
pause
%Now for the integral itself.  We shall do this
%as two nested integrals.  The inner integral
%integrates over:

% -b*(1-(x/a)^2)^.5 < y < b*(1-(x/a)^2)^.5

%and the outer integral integrates over -a < x < a.
%Actually, because of symmetry, we can just integrate
%over the positive quadrant and multiply the result
%by four, this being four times faster!
pause
%Our program is complicated a bit by the
%requirement that any function fed to quad8.m 
%must be in vector form - it must return a vector
%of output values for a vector of input values.
%You can get around this if you write our own
%quadrature routine, but it is not too much of
%a problem.
pause
%In the main program we just integrate over the
%outer loop:
integ=4*quad8('inner',0,a)
%where the inner loop is given in the attached
%function inner.m, which returns the integral
%over y for an array of input x's.
pause
%Note that the points are unevenly spaced:  most
%of the points are located near the end of the
%ellipsoidal region, where the function value is
%small and isn't changing much!  This is a consequence
%of the mapping we are doing using the above integration
%method!  An optimized 2-D gaussian quadrature routine
%would probably give us just as accurate result with
%only a few function evaluations!
hold off
echo off



_____________________________



function out=inner(x)
global a b xpass
%We use the variable xpass to pass
%the values of x to the function we
%are integrating.
n=length(x);
out=zeros(size(x));
for i=1:n
 xpass=x(i);
 ylim=b*(1-xpass^2/a^2)^.5;
 out(i)=quad8('funin',0,ylim,1e-4);
end
drawnow



_____________________________



function out=funin(y)
global xpass
out=1./exp(y.^2+xpass^2);
%Comment this next bit of graphics out
%if you want the program to run faster.
for i=1:length(y)
 plot(xpass,y(i),'o')
end