Radioactive Decay from Atom 1 to Atom 2 to Atom 3

I've got a monte carlo simulation of a radioactive decay from unstable atom 1 to unstable atom 2, which then decays to stable atom 3. I have a plot of the number of nuclides of each atom vs. time.

I'm trying to write a script for the analytical solution for the number of particles. The equations for each are N_1, N_2, N_3. They are currently commented out. One I have wrote this script I want a plot of them all together, and another plot comparing them to the monte carlo method. I'm struggling to do so at the minute. Any help would be much appreciated.

N0 = 1e5;  % initial number of atoms
lambda1 = 0.0001; %decays/s
lambda2 = 0.00005;
dt = 100; % time step (in seconds)
M = 1500; % total number of time steps
t = 0:dt:(M-1)*dt; % time series
p1 = lambda1*dt;  % probability of having a decay in the time dt
p2 = lambda2*dt;
% Define the array that will be populated with the number of atoms at each
% time step
N1 = zeros(M,1); N2 = zeros(M,1); N3 = zeros(M, 1); 
N1(1) = N0;

%N_1=N0*exp(-lambda1*(t));
%N_2=((lambda1*N0)/(lambda1-lambda2))*(exp(-lambda2*(t)))-(exp(-lambda1*(t)));
%N_3=((N0)/(lambda1-lambda2))*(lambda2*exp(-lambda1*(M*(t/dt))-lambda1*exp(-lambda2*(M*(t))+(lambda1-lambda2))));

tic
for j=2:M
    n_decaysN1toN2 = nnz(rand(N1(j-1),1) <= p1);
    n_decaysN2toN3 = nnz(rand(N2(j-1),1) <= p2);
      N1(j) = N1(j-1) - n_decaysN1toN2; % number of type 1 atoms that are left
      N2(j) = N2(j-1) + n_decaysN1toN2 - n_decaysN2toN3;
      N3(j) = N3(j-1) + n_decaysN2toN3;
end
toc

figure(1)
plot(t, [N1, N2, N3]), box off
legend({'# atoms type 1', '# atoms type 2', '#atoms type 3'})
title('title')
xlabel('time')
ylabel('number of particles')
legend boxoff

[Alexander Matthew]

Hi says,

can i ask you something?

n_decaysN1toN2 = nnz(rand(N1(j-1),1) <= p1);

n_decaysN2toN3 = nnz(rand(N2(j-1),1) <= p2);

Can you explain what that's mean and why you choose those code? Thanks

Answer 1

[David Goodmanson]

Hi says,

This is pretty close. In your expression for N_2, there is a misplaced right parenthesis. The factor in front is supposed to multiply both exponentials, and a correct expression is

N_2 = N0*(lambda1/(lambda1-lambda2))*(exp(-lambda2*t)-exp(-lambda1*t));

Visually it wasn't helping that there were unneeded parentheses around t and especially around the second exp, so they are gone. Here N0 is up front to make it consistent with the N_1 expression and to emphasize that everything in the whole problem is proportional to N0.

For N_3 you can work out an expression (that does not involve using dt), but it's a lot easier to take advantage of the fact that the particles are intraconverted, but they don't disappear. The total number of particles is always N0, and you already know N_1 and N_2, ...

Answer 2

[says]

Thank you, David Goodmanson! I'm trying now to create a plot that compares the analytical answers (N_1, N_2, N_3) with with the Monte Carlo Simulation (N1,N2,N3). I'm not sure how to do so though.

N0 = 1e5;  % initial number of atoms
lambda1 = 0.0001; %decays/s
lambda2 = 0.00005;
dt = 100; % time step (in seconds)
M = 1500; % total number of time steps
t = 0:dt:(M-1)*dt; % time series
p1 = lambda1*dt;  % probability of having a decay in the time dt
p2 = lambda2*dt;
% Define the array that will be populated with the number of atoms at each
% time step
N1 = zeros(M,1); N2 = zeros(M,1); N3 = zeros(M, 1); 
N1(1) = N0;

N_1 = zeros(M,1); N_2 = zeros(M,1); N_3 = zeros(M,1);
N_1(2)=N0;

N_1 = N0*exp(-lambda1*(t));
N_2 = N0*(lambda1/(lambda1-lambda2))*(exp(-lambda2*t)-exp(-lambda1*t));
N_3 = ((N0)/(lambda1-lambda2))*(lambda2*exp(-lambda1*(t))-lambda1*exp(-lambda2*(t)+(lambda1-lambda2)));

for k=2:M
N_1(k) = N_1(k-1);
N_2(k) = N_2(k-1) + n_decaysN1toN2 - n_decaysN2toN3;
N_3(k) = N_3(k-1) + n_decaysN2toN3;
N_1 = N0*exp(-lambda1*(t));
N_2 = N0*(lambda1/(lambda1-lambda2))*(exp(-lambda2*t)-exp(-lambda1*t));
N_3 = ((N0)/(lambda1-lambda2))*(lambda2*exp(-lambda1*(t))-lambda1*exp(-lambda2*(t)+(lambda1-lambda2)));
end

figure(2)
plot(t, [N_1, N_2, N_3]), box off
legend({'# atoms of type 1', '# atoms of type 2', '#atoms of type 3'})
title('Analytical Calculation of decay')
xlabel('time')
ylabel('number of particles')
legend boxoff

tic
for j=2:M
    n_decaysN1toN2 = nnz(rand(N1(j-1),1) <= p1);
    n_decaysN2toN3 = nnz(rand(N2(j-1),1) <= p2);
      N1(j) = N1(j-1) - n_decaysN1toN2; % number of type 1 atoms that are left
      N2(j) = N2(j-1) + n_decaysN1toN2 - n_decaysN2toN3;
      N3(j) = N3(j-1) + n_decaysN2toN3;
end
toc

figure(1)
plot(t, [N1, N2, N3]), box off
legend({'# atoms of type 1', '# atoms of type 2', '# atoms of type 3'})
title('Monte Carlo Calculation of Decay')
xlabel('time (s)')
ylabel('number of particles')
legend boxoff