Chapter 1 problem 1.5: The decay of two kinds of particles

homework python

Abstract

• First try to use python to draw figures.
• Solve problem about two kinds of particles decay.

Content

• Problem 1.5
  Consider again a decay problem with two types of nuclei A and B, but now suppose that nuclei of type A decay into ones of type B, while nucleiof type B decay into ones of type A. Strictly speaking, this is not a "decay" process, since it is possible for the type B nuclei to turn back into type A nuclei. A better analogy would be a reasonance in which a system can tunnel or move back and forth between two states A and B which have equal energies. The corresponding rate equations are
  

  $$\frac{\text{ d}N_A}{\text{d}t}=\frac{N_B}{\tau}-\frac{N_A}{\tau} \tag{1}$$
  $$\frac{\text{ d}N_B}{\text{d}t}=\frac{N_A}{\tau}-\frac{N_B}{\tau}\tag{2}$$
  where for simplicity we have assumed that the two types of decay are characterized by the same time constant, ${\tau}$. Solve this system of equations for the numbers of nuclei, ${N_A}$ and ${N_B}$, as functions of time .Consider different initial conditions, such as ${N_A}$ = 100, ${N_B}$ = 0, etc., and take ${\tau}$ = 1 s. Show that your numerical results are consistent with the idea that the system reaches a steady state which ${N_A}$ and ${N_B}$ are constant. In such a steady state, the time derivatives d${N_A}$/dt and d${N_B}$/dt should vanish.

• Ideas

  1. Nummerical Approach
  Taylor expansion for ${N_A}$
  

  $$N_A(t +△t)= N_A(t)+\frac{\text{ d}N_A}{\text{d}t}△t+\frac{1}{2!}\frac{\text{ d}^2N_A}{\text{d}t^2}(△t)^2+...\tag{3}$$
  Take O(1) cut, we have:
  $$N_A(t +△t) \approx N_A(t)+\frac{\text{ d}N_A}{\text{d}t}△t\tag{4}$$
  We can also have ${N_B}$ in the same way:
  $$N_B(t +△t) \approx N_B(t)+\frac{\text{ d}N_B}{\text{d}t}△t\tag{5}$$
  Δt is small enough to make the error terms negligible.
  So the equations (1) and (2) become:
  $$N_A(t +△t)=N_A+\frac{N_B}{\tau}△t-\frac{N_A}{\tau}△t\tag{6}$$
  $$N_B(t +△t)=N_B+\frac{N_A}{\tau}△t-\frac{N_B}{\tau}△t\tag{7}$$
  2.Program
  In this program I take △t = 0.05 and the duration of time is 5s.
  Codes are listed:

import pylab as pl
class uranium_decay:
     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.nuclei_a = [number_of_a]
         self.nuclei_b = [number_of_b]
         self.t = [0]
         self.tau = time_constant
         self.dt = time_step
         print("Initial number of nuclei ->", number_of_a)
         print("Time constant ->", time_constant)
         print("time step -> ", time_step)
         print("total time -> ", time_of_duration)
     def calculate(self):
        for i in range(100):
           tmpa = self.nuclei_a[i] + (self.nuclei_b[i] - self.nuclei_a[i]) / self.tau * self.dt
           self.nuclei_a.append(tmpa)
           tmpb = self.nuclei_b[i] + (self.nuclei_a[i] - self.nuclei_b[i]) / self.tau * self.dt
           self.nuclei_b.append(tmpb)
           self.t.append(self.t[i] + self.dt)
     def store(nuclei_a,nuclei_b,t):
           data = open('data.txt','w')
           for i in range(100):
               data.write(str(t[i]))
               data.write(' ')
               data.write(str(nuclei_a[i]))
               data.write(' ')
               data.write(str(nuclei_b[i]))
               data.write('\n')
               data.close
     def show_results(self):
       plot1, = pl.plot(self.t,self.nuclei_a,'g')
       plot2, = pl.plot(self.t,self.nuclei_b,'r')
       pl.xlabel('Time(s)')
       pl.ylabel('Number of Nuclei')
       pl.legend([plot1, plot2], ['NA', 'NB'], loc="best")
       pl.xlim(0, 5)
       pl.ylim(0, 100)
       pl.savefig('chapter1_1.5.png',dpi=144)
       pl.show()
b = uranium_decay()
b.calculate()
b.show_results()

• Result
  Figure:The number of nuclei of type A is marked by the green line and The number of nuclei of type B is marked by the red line.
  And here is the data I saved.
  From the figure we can see ${N_A}={N_B}$ finally and the system reaches a steady state which ${N_A}$ and ${N_B}$ are constant.

• Conclusion
  From this project I learned how to use Euler mothed to solve simple first order ordinary differential equation or equations and how to use python draw figures.