clear
echo on
%This program demonstrates the technique called
%Monte Carlo Integration.  It integrates the
%function exp(xy) over the circle x^2+y^2<1.  We
%also illustrate the technique of imbedding the
%domain, in which we actually perform the integ-
%ration over the square domain |x|<1,|y|<1 and
%simply define the function to be zero outside
%of the circle but within the square.  The error
%is just the standard deviation of the function
%divided by the square root of the number of points.
pause
%We first set the sum equal to zero:
sum=0;
sumsq=0; %this bit is for calculating the error

n=200;

%Now for a bit of graphics.  Note that the graphics
%-really- slows down the process!.  If you use a large
%number of points, be sure to comment this out!
pause
xx=[-1:.01:1];
yplus=(1-xx.*xx).^.5;
yminus=-(1-xx.*xx).^.5;
plot(xx,yplus,xx,yminus)
drawnow
hold on
pause
%Now for the actual integration.  We keep
%track of both the sum of the integrand over
%all of the points sampled, and the sum of
%the square of the integrand.  The latter is
%used in calculating the error in the integ-
%ration process.  So:
pause
for i=1:n
  x=rand*2-1;
  y=rand*2-1;
  plot(x,y,'o')
  drawnow
  if x^2+y^2 < 1,
    f=exp(x*y);
    sum=sum+f;
    sumsq=sumsq+f^2.;
    echo off
  end
end
echo on
hold off
pause
%We now calculate the integral and the error.
%Note that the integral is just the average
%value times the area of integration (including
%the zero regions outside the imbedded domain).
integ=sum/n*4;
error=4*(sumsq/n-(sum/n)^2)^.5/n^.5;
[integ,error]
pause
%Note that the relative error depends both on the
%magnitude of the function, the degree of variability,
%and the number of points used in the evaluation.
%For 200 points, the error is about 4% for this 
%integrand.  The correct answer is:
correct_int=3.2077
%which is within error of the mean (usually, at least -
%remember that the value you get will be different
%every time you run the program!).
echo off