Exercise_04 :

Consider a decay problem with two types of nuclei

冉峰 弘毅班 2014301020064

Abstract

• We consider a decay problem with two types of nuclei A and B, but suppose that the nucleiof A decay into ones of type B, while nuclei B decay into ones of type A, when the system reach a steady state the numberical are consistent in which $N_{A}$ and $N_{B}$ are constant.

Background

• We have that the corresponding rate equation of the decay between the nuclei A and B are:
  $$\frac{dN_{A}}{dt}=\frac{N_{B}}{\tau}-\frac{N_{A}}{\tau}$$
  $$\frac{dN_{B}}{dt}=\frac{N_{A}}{\tau}-\frac{N_{B}}{\tau}$$
  which for simplicity the problem, first we have assumed that two types of decay are characterized by the same time constant, $\tau$ and then for the different $\tau$.

The Main Body

Code 1

• code

import pylab as pl
class uranium_decay:
    def __init__(self, number_of_nuclei_A = 100,\
    number_of_nuclei_B = 0,time_constant_1 = 1,\
    time_constant_2 = 1, time_of_duration = 8,\
    time_step = 0.05):
        # unit of time is second
        self.n1_uranium = [number_of_nuclei_A ]
        self.n2_uranium = [number_of_nuclei_B ]
        self.t = [0]
        self.tau1 = time_constant_1
        self.tau2 = time_constant_2
        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 A ->", time_constant_1)
        print("Time constant B ->", time_constant_2)
        print("time step -> ", time_step)
        print("total time -> ", time_of_duration)
    def calculate(self):
        for i in range(self.nsteps):
            tmp1 = self.n1_uranium[i] - self.n1_uranium[i] /self.\
            tau1 * self.dt+self.n2_uranium[i] / self.tau2 * self.dt
            tmp2 = self.n2_uranium[i] - self.n2_uranium[i] /self.\
            tau2 *self.dt+self.n1_uranium[i] / self.tau1 * self.dt
            self.n1_uranium.append(tmp1)
            self.n2_uranium.append(tmp2)
            self.t.append(self.t[i] + self.dt)
    def show_results(self):
        pl.plot(self.t, self.n1_uranium)
        pl.plot(self.t, self.n2_uranium)
        pl.xlabel('time ($s$)')
        pl.ylabel('Number of Nuclei ')
        pl.legend(['$N_A$','$N_B$'])
        pl.show()
    def store_results(self):
        myfile = open('nuclei_decay_data.txt', 'w')
        for i in range(len(self.t)):
            print(self.t[i], self.n1_uranium[i], file = myfile)
            print(self.t[i], self.n2_uranium[i], file = myfile)
        myfile.close()
a = uranium_decay()
a.calculate()
a.show_results()
a.store_results()

Result

• running result for $N_{A0}=100$ ,$N_{B0}=0$，$\tau_A=\tau_B=1$
• running result for $N_{A0}=80$ ,$N_{B0}=0$,$\tau_A=\tau_B=1$
• running result for $N_{A0}=100$ ,$N_{B0}=50$,$\tau_A=\tau_B=1$
• running result for $N_{A0}=50$ ,$N_{B0}=100$,$\tau_A=\tau_B=1$
• running result for $N_{A0}=100$ ,$N_{B0}=0$,$\tau_A=1$,$\tau_B=2$
• running result for $N_{A0}=100$ ,$N_{B0}=0$,$\tau_A=2$,$\tau_B=1$
  

Code 2

We can see the results in one figure.

• code

import pylab as pl
class uranium_decay:
    def __init__(self, number_of_nuclei_A = 100,\
    number_of_nuclei_B = 0,time_constant_1 = 2,\
    time_constant_2 = 1, time_of_duration = 8,\
    time_step = 0.05,time_step1 = 0.10,number_of_nuclei_A1 = 100,number_of_nuclei_B1 = 0):
        # unit of time is second
        self.n1_uranium = [number_of_nuclei_A ]
        self.n2_uranium = [number_of_nuclei_B ]
        self.n11_uranium = [number_of_nuclei_A1 ]
        self.n21_uranium = [number_of_nuclei_B1 ]
        self.t = [0]
        self.tau1 = time_constant_1
        self.tau2 = time_constant_2
        self.dt = time_step
        self.dt1 = time_step1
        self.time = time_of_duration
        self.nsteps = int(time_of_duration // time_step + 1)
        self.nsteps1 = int(time_of_duration // time_step1 + 1)
        print("Initial number of nuclei A ->", number_of_nuclei_A)
        print("Initial number of nuclei B ->", number_of_nuclei_B)
        print("Initial number of nuclei A ->", number_of_nuclei_A1)
        print("Initial number of nuclei B ->", number_of_nuclei_B1)
        print("Time constant A ->", time_constant_1)
        print("Time constant B ->", time_constant_2)
        print("time step -> ", time_step)
        print("time step -> ", time_step1)
        print("total time -> ", time_of_duration)
    def calculate(self):
        for i in range(self.nsteps):
            tmp1 = self.n1_uranium[i] - self.n1_uranium[i] /self.\
            tau1 * self.dt+self.n2_uranium[i] / self.tau2 * self.dt
            tmp2 = self.n2_uranium[i] - self.n2_uranium[i] /self.\
            tau2 *self.dt+self.n1_uranium[i] / self.tau1 * self.dt
            tmp11 = self.n11_uranium[i] - self.n11_uranium[i] /\
            self.tau1 * self.dt1+self.n21_uranium[i] /\
            self.tau2 * self.dt1
            tmp21 = self.n21_uranium[i] - self.n21_uranium[i] /\
            self.tau2 *self.dt1+self.n11_uranium[i] /\
            self.tau1 * self.dt1
            self.n1_uranium.append(tmp1)
            self.n2_uranium.append(tmp2)
            self.n11_uranium.append(tmp11)
            self.n21_uranium.append(tmp21)
            self.t.append(self.t[i] + self.dt)
    def show_results(self):
        pl.plot(self.t, self.n1_uranium)
        pl.plot(self.t, self.n2_uranium)
        pl.plot(self.t, self.n11_uranium)
        pl.plot(self.t, self.n21_uranium)
        pl.xlabel('time ($s$)')
        pl.ylabel('Number of Nuclei ')
        pl.legend(['$N_A$','$N_B$'])
        pl.show()
    def store_results(self):
        myfile = open('nuclei_decay_data.txt', 'w')
        for i in range(len(self.t)):
            print(self.t[i], self.n1_uranium[i], file = myfile)
            print(self.t[i], self.n2_uranium[i], file = myfile)
            print(self.t[i], self.n11_uranium[i], file = myfile)
            print(self.t[i], self.n21_uranium[i], file = myfile)
        myfile.close()
a = uranium_decay()
a.calculate()
a.show_results()
a.store_results()

Result

• running result for $N_{A0}=100$ ,$N_{B0}=0$，$\tau_A=\tau_B=1$,$dt_1=0.05$,$dt_2=0.1$
• running result for $N_{A01}=100$, $N_{A02}=80$,$N_{B0}=0$，$\tau_A=\tau_B=1$,$dt_1=0.05$,$dt_2=0.05$
• running result for $N_{A01}=100$, $N_{A02}=80$,$N_{B0}=0$，$N_{B02}=50$,$\tau_A=\tau_B=1$,$dt_1=0.05$,$dt_2=0.05$
• running result for $N_{A01}=100$, $N_{A02}=80$,$N_{B0}=0$，$N_{B02}=50$,$\tau_A=1$,$\tau_B=2$,$dt_1=0.05$,$dt_2=0.05$
• running result for $N_{A01}=100$, $N_{A02}=80$,$N_{B0}=0$，$N_{B02}=50$,$\tau_A=2$,$\tau_B=1$,$dt_1=0.05$,$dt_2=0.05$
  
  

From the figure we can see that the system reaches a steady state in which $N_A$ and $N_B$ are constant.In such a steady state, the time derivatives $\frac{dN_{A}}{dt}$ and$\frac{dN_{B}}{dt}$ wanish.

Conclusion

• With the learning of the matplotlib, we can use it to solve the problem more easy and directly. we can easy to have a picture of the problem which can help us to knowthe answer more and help us to know how to consider the problem when we can not solve the problem easyly.

Thanks For

• Lesson plan of Cai Hao:Chapter 1