`郭忠智` `学号2014301020087`

Problem1.5

一、摘要

本次作业应用Python编写程序，解决了一个微分方程组——两种原子核A和B吸收和释放核子的相互转变，计算出每个时刻两种原子核数量并将其结果通过画图展现出来，可以利用数值运算可以方便地解决微分方程。

二、背景

利用计算机解决物理问题，在今天的物理研究里至关重要，利用编程语言，可以将科研工作者集中精力于物理问题，避免花过多时间在解决数学问题上。第一章通过介绍了一种简单的解决微分方程的方法——级数展开及数值逼近法。这种方法简单易懂，但是需要注意精度问题，取不同的时间间隔，会产生不同大小的误差，在具体的实际问题中，需要细心考虑所取步长对结果的影响。

三、正文

原问题所给出的微分方程组为：

$\frac{d{N_A}}{dt}=\frac{N_B}{\tau}-\frac{N_A}{\tau}$

$\frac{d{N_B}}{dt}=\frac{N_A}{\tau}-\frac{N_B}{\tau}$
（1）用微积分解该方程组，得到：

$N_A(t)=C_1+C_2e^{-2t/\tau}$

$N_B(t)=C_1-C_2e^{-2t/\tau}$
代入

$\tau=1s$ 、以及初值条件

时 、

$t=0时N_A=100、N_B=0$ ,解得：

$C_1=50$

$C_2=50$
所以

$N_A(t)=50+50e^{-2t/\tau}$

$N_B(t)=50-50e^{-2t/\tau}$
（2）利用数值逼近方法，写出的python程序代码如下：

# -*- coding: utf-8 -*-
"""
Spyder Editor
This is a temporary script file.
"""
# A decay problem with two types of nuclei A and B.
# Simulation of the decay.
# Program to problem_1.5 on the book Computational Physics.
import pylab as pl
class nuclei_decay:
    """
    A decay problem with two types of nuclei A and B.
    Simulation of the decay.
    Program to problem_1.5 on the book Computational Physics.
    """
    def __init__(self, number_of_A = 100, number_of_B = 0, time_constant = 1, time_of_duration = 5, time_step = 0.05):
        # unit of time is second
        self.N_A = [number_of_A]
        self.N_B = [number_of_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_A)
        print("Initial number of nuclei B ->", number_of_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_A[i] + (self.N_B[i] - self.N_A[i]) * self.dt / self.tau
            tmp_B = self.N_B[i] + (self.N_A[i] - self.N_B[i]) * self.dt / self.tau
            self.N_A.append(tmp_A)
            self.N_B.append(tmp_B)
            self.t.append(self.t[i] + self.dt)
    def show_results(self):
        plot1, = pl.plot(self.t, self.N_A, 'r')
        plot2, = pl.plot(self.t, self.N_B, 'g')
        pl.xlabel('time ($s$)')
        pl.ylabel('Number of Nuclei')
        pl.title("Number's change of nuclei A and nuclei B")
        pl.legend([plot1, plot2], ['nuclei A', 'nuclei B'], loc="best")
        pl.show()
    def store_results(self):
        myfile = open('nuclei_decay_data.txt', 'w')
        for i in range(len(self.t)):
            myfile.write(str(self.t[i]))
            myfile.write(str(self.N_A[i]))
        myfile.close()
        print myfile
a = nuclei_decay()
a.calculate()
a.show_results()
a.store_results()

（3）对比结果：改写show_results()函数，加入利用微积分所求得的结果的图形，改写的代码如下：

def show_results(self):
        import math
        L = []
        for x in range(501):
            L.append(x * 0.01)
        x = L
        y = []
        z = []
        for i in range(501):
            a = 50 + 50 * math.e ** (-2 * x[i])
            b = 50 - 50 * math.e ** (-2 * x[i])
            y.append(a)
            z.append(b)
        plot_1, = pl.plot(x, y, 'r:')
        plot_2, = pl.plot(x, z, 'g:')
        plot1, = pl.plot(self.t, self.N_A, 'r')
        plot2, = pl.plot(self.t, self.N_B, 'g')
        pl.xlabel('time ($s$)')
        pl.ylabel('Number of Nuclei')
        pl.title("Number's change of nuclei A and nuclei B")
        pl.legend([plot1, plot2, plot_1, plot_2], ['nuclei A', 'nuclei B', 'nuclei A for calculus', 'nuclei B for calculus'], loc="best")
        pl.show()

四、结论

运行（2）中所写代码，并将（1）（2）中所得结果作图，得到如下结果：
1 5
1 5

可以看到，利用数值逼近的方法可以得到所求结果的较好近似