% Sauer Computer Problem 9.1.3 (c)
%
% Gör Monte Carlo typ 2 (slumpa fram punkter inom en kvadrat och
% räkna hur många som hamnar mellan parablerna).
%
% Joar Bagge 2016-12-05

%%
% Illustration av Monte Carlo typ 2

clear

P1 = @(x) x.^2 - x + 1/2;
P2 = @(x) -x.^2 + x + 1/2;

figure(4), clf
fplot(P1, [0,1]), hold on
fplot(P2, [0,1])

n = 100; % antal punkter som slumpas fram
X = rand(1,n); % slumpade x-koordinater mellan 0 och 1
Y = rand(1,n); % slumpade y-koordinater mellan 0 och 1

% Följande blir en vektor med n element, där det är en etta om
% motsvarande (x,y)-punkt ligger mellan parablerna, annars en
% nolla.
in = Y < P2(X) & Y > P1(X);
out = ~in; % negerar den boolska vektorn
scatter(X(in), Y(in), 'm')
scatter(X(out), Y(out), 'k')
axis([0 1 0 1])

% Dela antalet punkter inuti arean (mellan parablerna) med det
% totala antalet punkter. Multiplicera med kvadratens area, men
% den är 1 i detta fall.
A = sum(in) / n

%%
% Test för flera olika n

ns = 10.^(2:6);
nmax = max(ns);
Xall = rand(1,nmax);
Yall = rand(1,nmax);
As = [];
for n=ns
  X = Xall(1:n);
  Y = Yall(1:n);
  in = Y < P2(X) & Y > P1(X);
  A = sum(in) / n;
  As = [As, A];
end

exact = 1/3; % exakt värde på arean
figure(5)
semilogx(ns, As, '*-') % uppskattningar
hold on
semilogx(ns, 0*ns+exact, 'k--') % exakt svar
title('Integraluppskattning')

figure(6)
loglog(ns, abs(As-exact), '*-') % fel i uppskattningarna
hold on
loglog(ns, 1./sqrt(ns), 'k--') % referenslinje
title('Fel')
% Kör koden flera gånger så blir det fler kurvor i plottarna!
