%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution to CS&E Your Homework Assignment Project 10
% Multidimensional Integration: Partition and Conquer
%
% In this example, we estimate the area of a quarter circle
% of radius 0.5, bounded by the positive x and y axes and
% the curve sqrt(radius^2 - x^2).
%
% The estimate is the average value of sqrt(radius^2 - x^2) 
%  for x in [0,radius]
% times the length of the interval [0,radius]
%
% problem4.m Dianne P. O'Leary   09/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nn = 100000;
ni = 19;

radius = .5;
trueanswer = radius^2*pi/4;
disp(sprintf('Compute the area of the quarter circle of radius %f',...
   radius))
disp(' ')
disp(sprintf('True area: %5.8f',trueanswer));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the quasirandom Monte Carlo approximations
% for 10, 100, 1000, 10000, 100000 pts.
% We generate and save nn quasirandom numbers Xuqr in the
% interval [0,radius].
% This saves us some time and makes sure that each of the
% experiments uses the same sequence.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Xuqr = radius* quasirand(1,nn,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method 4.  Approximate the integral by the average
% function value times the length of the interval.
% Compute these Monte Carlo approximations
% for 10, 100, 1000, 10000, 100000 quasi-random points.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(sprintf('Method 4: quasirandom Monte Carlo'));
fvalsqr = sqrt(radius^2 - Xuqr(:,1).^2);
disp('         n   estimate        error        error*n')
for i=1:5,
  n = 10^i;
  answerqr(i) = radius * sum(fvalsqr(1:n))/n;
  error = answerqr(i) - trueanswer;
  disp( ...
   sprintf('%10.0f   %f   %e   %e', ...
        n, answerqr(i), error, error*n) )
end