% Skript som estimerer et integral ved Monte Carlo-metoden
% Selve integralet, og den analytiske løsninga av det, er
% spesifisert i starten av skriptet - i tillegg til
% maksimalt antall trekk.
% I tillegg til å estimere integralet, plotter skriptet
% estimatet for økende antall trekk. På den måten får vi 
% sett hvordan estimatet nærmer seg det faktiske integralet.

% Integrasjonsgrenser
a=-1; 
b=3;

% Integrand
funk=@(x) x*sin(x^2);

% Kjent analytisk svar:
Fasit=(cos(1)-cos(9))/2;

% Maksimalt antall trekk:
Nmax=input('Gi maksimalt antall trekk: ');

% Initierer summen
MCsum=0;

% Lager klart plott for estimatene - sammen med eksakt svar
plot([1 Nmax],Fasit*[1 1],'k--')
hold on

% Vi utfører trekkene og legger til summen for hver gang
for n=1:Nmax
  x=a+(b-a)*rand;
  f=funk(x);
  MCsum=MCsum+f;
  Estimat=(b-a)*MCsum/n;          % Integral-estimat
  plot(n,Estimat,'b.')            % Plotter estimat
end
hold off

% Skriver svaret til skjerm
Estimat

