MATLAB - Radioactive Decay by Monte Carlo Method

% _________________________________________________________________________
% To simulate radio active decay by Monte Carlo method
% It gives nuclides number upto 10000 time steps
% By Mahesha MG, MIT
% Date: 15/02/2014
% _________________________________INPUT___________________________________
clc;
clear;
lambda=input('Enter decay constant in /s: ');
n=input('Enter total number of nuclides: ');
dt=0.0001/lambda;
nt=zeros(10000,1); %nt holds MC result
nta=zeros(10000,1); %nta holds analytic result
tmax=10000*dt;
t=(0:dt:tmax)';
nta=n*exp(-lambda*t); %Analytic solution
nt(1)=n;
% __________________________Monte Carlo Method_____________________________
for i=2:10001
    random=rand(n,1);
    count=length(random(random<=0.0001));
    n=n-count;
    if(n<=0)
        break;
    end
    nt(i)=n;
end
% _________________________________________________________________________
plot(t,nt,'r',t,nta,'b');
xlabel('Time');
ylabel('No. of nuclides');

Published with MATLAB® R2013b