作业四

计算物理

---

1.摘要

运用Euler method,编写python程序，利用matplotlib绘图，模拟粒子衰变过程

2背景

通过以学习的模拟单个衰变的方法，考虑两个粒子A和B的衰变过程。假设A粒子衰变成B粒子的同时，B粒子也衰变成A粒子。模拟俩粒子组成的系统达到平衡的过程。

3正文

1.方程

已给出粒子数随时间改变速率的方程

$$\frac{dN_A}{dt}=\frac{N_B}{\tau}-\frac{N_A}{\tau}$$
$$\frac{dN_B}{dt}=\frac{N_A}{\tau}-\frac{N_B}{\tau}$$
运用泰勒近似得到的方程
$$N(t+\Delta t)\approx N(t)+\frac{dN}{dt}\Delta t$$
可得
$$N_A(t+\Delta t)\approx N_A(t)+(\frac{N_B}{\tau}-\frac{N_A}{\tau})\Delta t$$
$$N_B(t+\Delta t)\approx N_B(t)+(\frac{N_A}{\tau}-\frac{N_B}{\tau})\Delta t$$
同时解出微分方程
$$N_A = \frac{N_A+N_A e^-2t}{2}$$
$$N_A = \frac{N_B-N_B e^-2t}{2}$$

2程序代码

import pylab as pl
class uranium_decay_A_and_B:
    def __init__(self,number_of_nuclei_A = 100,number_of_nuclei_B = 0, time_constant=1,time_of_duration = 5,time_step=0.1):
        self.n_uranium_A = [number_of_nuclei_A]
        self.n_uranium_B = [number_of_nuclei_B]
        self.t = [0]
        self.tau = time_constant
        self.dt = time_step
        self.time = time_of_duration
        self.nsteps = int(time_of_duration // time_step + 1)
        print("Initial number of nuclei A ->", number_of_nuclei_A)
        print("Initial number of nuclei B ->", number_of_nuclei_B)
        print("Time constant ->", time_constant)
        print("time step -> ", time_step)
        print("total time -> ", time_of_duration)
    def calculate(self):
        for i in range(self.nsteps):
            tmp_A = self.n_uranium_A[i] + (self.n_uranium_B[i] - self.n_uranium_A[i])/self.tau*self.dt
            tmp_B = self.n_uranium_B[i] + (self.n_uranium_A[i] - self.n_uranium_B[i])/self.tau*self.dt
            self.n_uranium_A.append(tmp_A)
            self.n_uranium_B.append(tmp_B)
            self.t.append(self.t[i] + self.dt)
    def show_results(self):
        pl.plot(self.t,self.n_uranium_A,'r')
        pl.plot(self.t,self.n_uranium_B)
        pl.xlabel('time($s$)')
        pl.ylabel('Number of Nuclei')
        pl.show()
a = uranium_decay_A_and_B()
a.calculate()
a.show_results()
from math import *
class exact_results_check(uranium_decay_A_and_B):
    def show_results(self):
        self.et_A = []
        self.et_B = []
        for i in range(len(self.t)):
            temp_A = self.n_uranium_A[0]/2 + exp(-self.t[i])*self.n_uranium_A[0]/2
            temp_B = self.n_uranium_A[0]/2 - exp(-self.t[i])*self.n_uranium_A[0]/2
            self.et_A.append(temp_A)
            self.et_B.append(temp_B)
        pl.plot(self.t,self.et_A)
        pl.plot(self.t,self.et_B)
        pl.plot(self.t,self.n_uranium_A,'--')
        pl.plot(self.t,self.n_uranium_B,'--')
        pl.xlabel('time ($s$)')
        pl.ylabel('Number of Nuclei')
        pl.xlim(0, self.time)
        pl.show()
b = exact_results_check()
b.calculate()
b.show_results()

3结果

模拟结果如下

与准确结果比较如下，虚线表示近似模拟，实线是准确结果

4结论

Euler method近似计算的结果与准确结果有一定误差，从图形中很容易得出，且近似的值总比准确值小，系统趋于平衡的准确时间大于近似的计算的时间。

5

感谢李云飞同学的公式输入法总结