计算科学 - 在MATLAB中使用蒙特卡罗方法模拟连续衰减 - 吾爱随笔录

我正在尝试使用 MATLAB 通过使用蒙特卡罗方法来模拟原子衰变过程。过程如下：

假设原子 1 衰变为原子 2，原子 2 又衰变为稳定的类型 3 原子。1 和 2 的衰变常数和。 $\lambda_1$ $\lambda_2$

假设在，和。 $t = 0$ $N_1 = N_0$ $N_2 = N_3 = 0$

$\lambda_1=0.0001$ , , 进程为。 $\lambda_2=0.00005$ $t_{\rm MAX}$ $150000\rm s$

我对蒙特卡洛方法不是很熟悉，所以我只设法模拟了单次衰减过程。当谈到连续衰减时，我完全迷失了。

我在下面附上了我的单衰变脚本。希望能在这里得到帮助。

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;
        end    
    end
    N(j) = N(j-1)-nd; % number of atoms left at the step 2
end

plot(t, N)