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.  
xx=[-1:.01:1];
yplus=(1-xx.*xx).^.5;
yminus=-(1-xx.*xx).^.5;
figure(1)
plot(xx,yplus,xx,yminus)
set(gca,'Fontsize',18)
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 integration
%process.  So:
pause
for i=1:n
  x=rand*2-1;
  y=rand*2-1;
  figure(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