第四次作业

---

一 摘要

用Python语言编写解一个二元一阶齐次常微分方程组的程序，并解决实际问题。

二 背景介绍

用欧勒法（Euler method）来解常微分方程。对函数N(t)作泰勒展开：

$N(t)=N(0)+{\frac{\mathrm{d} N}{\mathrm{d} t}}\Delta t+\frac{\mathrm{d}^2 N}{\mathrm{d} t^2}\Delta t^2+\cdots$

取一阶近似：$N(t)=N(0)+{\frac{\mathrm{d} N}{\mathrm{d} t}}\Delta t+\frac{\mathrm{d}^2 N}{\mathrm{d} t^2}\Delta t^2$

对于题目方程：

$$\frac{\mathrm{d} N_{A}}{\mathrm{d} t}=\frac{N_{B}}{\tau }-\frac{N_{A}}{\tau }$$
$$\frac{\mathrm{d} N_{B}}{\mathrm{d} t}=\frac{N_{A}}{\tau }-\frac{N_{B}}{\tau }$$

由欧勒法有：

$$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$$
其中$\tau$是一个时间常数，我们可以用这组方程近似求出$\Delta t$时间后$N_{A},N_{B}$的值.

三 正文

import pylab as pl
class uranium_decay:
    def __init__(self, number_of_nuclei_A = 100, number_of_nuclei_B=0, time_constant = 1, time_of_duration = 5, time_step = 0.05):
        # unit of time is second
        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):
            tmpA = self.n_uranium_A[i] + (self.n_uranium_B[i] / self.tau -self.n_uranium_A[i] / self.tau)* self.dt
            tmpB = self.n_uranium_B[i] + (self.n_uranium_A[i] / self.tau -self.n_uranium_B[i] / self.tau)* self.dt
            self.n_uranium_A.append(tmpA)
            self.n_uranium_B.append(tmpB)
            self.t.append(self.t[i] + self.dt)
    def show_results(self):
        pl.plot(self.t, self.n_uranium_A)
        pl.xlabel('time ($s$)')
        pl.ylabel('Number of Nuclei A')
        pl.show()
        pl.plot(self.t, self.n_uranium_B)
        pl.xlabel('time ($s$)')
        pl.ylabel('Number of Nuclei B')
        pl.show()
a = uranium_decay()
a.calculate()
a.show_results()

最后$N_{A}和N_{B}$都趋近于50

四、结论

最后$N_{A}和N_{B}都趋近于$50
因为$N_{A}和N_{B}$都有相同的衰变常数，且A到B和B到A的过程是可逆的，所以最后会趋近于100的中值，即$N_{A}=N_{B}$.

五 感谢

感谢蔡老师的PPT 及杜威同学帮我解决了图片上传的问题。