中子穿墙问题的MonteCarlo求解方法

2024-04-11 20:07:25

中子穿墙问题的MonteCarlo求解方法

一.问题的提出

如下图所示，代表一个中子穿过用于屏蔽中子的铅墙的示意图。铅墙的高度远大于左右的厚度。设中子垂直由左端进入铅墙，铅墙中前行一个单位距离后与一个铅原子碰撞。此时会改变方向，但是此时改变的方向是任意的，并且能继续运行一个单位后与另一个铅原子碰撞。这样下去，如果中子在铅墙里消耗掉所有的能量或者从左端溢出即被视为中子被铅墙挡住。但是，如果中子由右端溢出，那么我们就视作中子溢出。假设穿墙的厚度为5个单位，能量为可以运行7个单位，每运行一个单位后视为中子会发生碰撞，求中子溢出的概率？

二.问题的分析

我们关心的是每一次碰撞后中子在 x x x轴方向前进了多少。所以行进方向是正/负对结果影响不大，我们只考虑正值情况。用计算机抽样 θ i ∼ U ( 0 , π ) i = 1 , 2..7 \theta_i \sim U(0,\pi)\quad i = 1,2..7 θi∼U(0,π)i=1,2..7,看经过7次碰撞后，有多少能超过墙的右端。

S t e p 1 : Step1: Step1:随机产生 θ i ∼ U ( 0 , π ) i = 1 , 2..7 , \theta_i \sim U(0,\pi) \quad i = 1,2..7, θi∼U(0,π)i=1,2..7,令初值 C o u n t = 0 。 Count = 0。 Count=0。

S t e p 2 : Step2: Step2:求 S = ∑ i j c o s ( θ i ) j = j + 1 S = \sum_{i }^j cos(\theta_i)\quad j = j+1 S=∑ijcos(θi)j=j+1.

S t e p 3 : Step3: Step3:判断如果 S ≥ 5 S \geq 5 S≥5或者 S ≤ 0 S \leq 0 S≤0。

S t e p 4 : Step4: Step4:如果 S ≥ 5 , C o u n t = C o u n t + 1 S \geq 5,Count = Count+1 S≥5,Count=Count+1。否则转向 S t e p 1 Step1 Step1,直到到达总的实验次数。

三.问题的求解

求解用 m a t l a b matlab matlab求解，易如反掌。

N = 1e6;%设置总实验次数
countSuccess = 0;
countLeft = 0;
countLose = 0;
for i = 1:N
    Theta = unifrnd(0,pi,[1,7]);
    k = 1;
    S = 0;
    while k<=7
        S = S + cos(Theta(k));
        if S >= 5
            countSuccess  = countSuccess+1;
            break;
        elseif S<=0
            countLeft = countLeft + 1;
            break;
        end
        k = k + 1;
    end
    if k == 8 %此时正好能量消耗完
       countLose = countLose + 1;
    end
end
Psuccess = countSuccess/N;
Pleft = countLeft/N;
Plose = countLose/N;
disp(['中子成功从右端溢出的概率：',num2str(Psuccess)]);
disp(['中子成功从左端溢出的概率：',num2str(Pleft)]);
disp(['中子能量耗尽的概率：',num2str(Plose)]);

结果为:

中子成功从右端溢出的概率：0.002714
中子成功从左端溢出的概率：0.79054
中子能量耗尽的概率：0.20675

说明中子最后能从此铅墙溢出的概率为 0.27 % 0.27\% 0.27%。

码农公寓