% ### EXintegrateMC1.m ### 10.08.14 (C. Bergevin) % This code demonstrates integration via a two simple Monte carlo methods % for various functions over a specified interval % ** Method 1 ** - Computes Np (uniformly distributed) (x,y) pairs over the relevant range, % than ascertains what fraction would contribute to the integral (i.e., which % ones are 'under' the curve). % ** Method 2 ** - From the uniformly-distributed x-values, the function is % evaluated and a mean such is subsequently determined (which then allows % the definite integral to be determined) % --> Both are compared to a trapz.m estimate (i.e., Riemann sum) clear % ----------------------------- % User Inputs N= 100; % number of 'data' points Np= 1000; % number of random samples to compute intLimts= [-3 8]; % limits to evaluate integral over % === % user specifies function to integrate % generic model parameters AA= 2; % model parameters BB= 2; CC= 0.5; % function %z= @(x) AA./(BB-exp(-CC*x.^2)); % Gaussian z= @(x) AA*sin(x)-3; % sinusoid %z= @(x) AA*exp(-x)-BB; % decaying exponential % ------------------------------------------------- x=linspace(intLimts(1),intLimts(2),N); % create x-values y= z(x); % determine associated y-values % *** figure(1); clf; % plot the function plot(x,y,'b-','LineWidth',2); hold on; grid on; xlabel('x'); ylabel('y'); % *** % Integrate via built-in Matlab function trapz.m a1 = trapz(x,y); msg = ['Integral calculated by trapz.m (i.e., Riemann sums)= ' num2str(a1,10)]; disp(msg); % *************** % Integrate via a brute-force Monte carlo routine % 1. Create a bounding box about the function if (min(y)>=0), yL=0; else yL= min(y); end % handle case if entire curve sits above x-axis if (max(y)<=0), yH=0; else yH= max(y); end % handle case if entire curve sits below x-axis yD= yH- yL; xL= min(x); xH= max(x); xD= xH- xL; % 2. Create uniform distribution of points inside that box; then for each % point determine whether it falls between the curve and x-axis (if so, % counter is updated +1) cntP= 0; cntN= 0; for nn=1:Np xT= xD*rand(1)- abs(xL); yT= yD*rand(1)- abs(yL); zT= z(xT); % evaluate the function at all the random x-values storeP(nn,:)= [xT yT]; % store away points (for visualization) % consider more genral case for instances above or below x-axis if (yT>=0) if (yT<=zT), cntP= cntP+1; end end if (yT<0) if (yT>=zT), cntN= cntN+1; end end zComp(nn)= zT; % compile all values of evaluated function end % estimate area as ratio determined via flagged counts areaMC= ((cntP-cntN)/Np)*yD*xD; % diff. here accounts for area being above or below x-axis msg = ['Integral calculated via Monte carlo Method 1 (area ratios) = ' num2str(areaMC)]; disp(msg); plot(storeP(:,1),storeP(:,2),'r+'); title('Method 1 (area ratios)') % visualize the points used % also report back estimate via average value msg = ['Integral calculated via Monte carlo Method 2 (average value) = ' num2str(mean(zComp)*diff(intLimts))]; disp(msg);