I'm trying to use MATLAB to simulate an atom decay process by using Monte Carlo approach. The process is as follow:

Suppose that atom 1 decay to atom 2, which, in turn, decay to stable atoms of type 3. The decay constants of 1 and 2 are $\lambda_1$ and $\lambda_2$, respectively.

Assume that at $t = 0$, $N_1 = N_0$ and $N_2 = N_3 = 0$.

$\lambda_1=0.0001$, $\lambda_2=0.00005$, $t_{\rm MAX}$ for the process is $150000\rm s$.

I'm not really familiar with Monte Carlo method, so I have only managed to simulate the single decay process.When it come to consecutive decay I got completely lost.

I attached my script for single decay below. Hope to get help here.

N0=input('Enter total number of nuclides: '); %Determine N0
lambda1 = 0.0001; %decays/s
dt = 100; %time step
M = 1500; % total number of time steps
t = 0:dt:(M-1)*dt; %time series
p = lambda1*dt; %probability of having a decay in the time dt

% Define the array that will be populated with the number of atoms at each time step N = zeros(M,1);

N(1) = N0;

for j=2:M % Loop over time
    nd = 0; % number of decays
    for i=1:N(j-1)
        %generate a single random number  
        r = rand();
        %binary value =1 if the condition is True -> Atom has decayed
        rb = r < dt*lambda1; 
        if rb == 1
           nd = nd+1;
    N(j) = N(j-1)-nd; % number of atoms left at the step 2

plot(t, N)

1 Answer 1


The decay of radioisotopes can be described by a set of coupled first order differential equations $$\frac{dN_{i}(t)}{dt} = -\lambda_{i}N_{i}(t) + \sum_{j \neq i}\lambda_{j\to i} N_{j}(t)$$ where $\lambda_{j\to i}$ represents the decay constant associated with the decay of isotope $j$ to isotope $i$, $\lambda_{i}$ is the total decay constant of isotope $i$ and $N_{i}(t)$ is the time-dependent concentration of isotope $i$.

This can be easily written in matrix form $$\frac{d\bar{N}(t)}{dt} = \bar{A}\bar{N}(t)$$ where $\bar{A}$ is typically very sparse. The numerical difficulties associated with this set of coupled differential equations are twofold:

  1. The system can be large (in our applications we track about 3500 isotopes)
  2. The system is very stiff due to the big differences between decay constants (some isotopes decay with half-lives in milliseconds, others in millions of years)

Luckily, there exist solvers than can handle these problems: solvers for stiff sets of ODEs like implicit Runge-Kutta methods (RADAU) or methods based on the matrix exponential method (CRAM). The application of RADAU to this problem is described in Advanced Method for Calculations of Core Burn-Up, Activation of Structural Materials, and Spallation Products Accumulation in Accelerator-Driven Systems (note that I'm co-author of this paper).

I see no reason nor advantage to use a Monte Carlo method to solve this problem. Do you intend to model the stochastic nature of radioactive decay? Why?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.