@computationalphysics-2014301020090 2016-10-08T13:31:05.000000Z 字数 6152 阅读 176

The 4th homework

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

Abstract

● The use of Python mapping function to study the equilibrium state of radioactive particles.
● Solve problem about two kinds of particles decay.
● The equilibrium state under different initial values is considered in the decay model.

background

To solve first-order differential problems,there is the Euler method which is based on Taylor expansion:

$f(x_0+△x)=f(x_0)+\frac{f'(x_0)}{1!} \ △x+\frac{f''(x_0)}{2!}△x^2+......+\frac{f^{(n)}(x_0)}{n!}△x^n+Rn(x)$
where

$f(x_0)$ is the value of our function at

$x=x_0$ .If we take

$△x$ to be small,then it is usually a good approximation to simply ignore the terms that involve second and higher powers of

$△x$ ,leaving us with:

$f(x_0+△x)≈f(x_0)+\frac{ df(x_0)}{dx} |_{x=x_0} △x$
For problem 1.5,it dives us two relevant

$N_A$ and

$N_B$ which are characterised by two first-order differential equations.We can raplace

$dN_A \over dt$ with former equation,so,

$N_A(T+△t)≈N_A(t)+({dN_B\over τ}-{dN_A \over τ})△t$
It's the same with

$N_B$

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

${dN_A\over dt}={N_B \over τ}-{N_A \over τ} \qquad\qquad {dN_B\over dt}={N_A \over τ}-{N_B \over τ}$
where for simplicity we have assumed that the two types of decay are characterized by the same time constant, . 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

$τ=1s$ . Show that your numerical results are consistent with the idea that the system reaches a steady state which

$N_AandN_B$ are constant. In such a steady state, the time derivatives

${dN_B\over dt}and{dN_A \over dt}$ should vanish.

● solution

Taylor expansion for $N_A$ ：

$N_A(T+△t)≈N_A(t)+({dN_B\over τ}-{dN_A \over τ})△t$
Similarly，

$N_B(T+△t)≈N_B(t)+({dN_A\over τ}-{dN_B \over τ})△t$

● Program

take $△t$ = 0.05 and the duration of time is 5s.
Codes are listed:

# -*- coding: utf-8 -*-
"""
Created on Fri Oct 07 2016
@author: 胡胜勇
"""
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： Problem 1.5.png
From the figure we can see finally and the system reaches a steady state which and are constant.

now we change the initial conditions：
● NA=50, NB=0
● NA=50, NB=20
● NA=80, NB=40
It's obvious that no matter how many $N_AandN_B$ about the initial conditions, the more one tend to decrease and the less one tend to increase and they will both reach the equilium number after a long time.And the equilium number is the average value of the initial number of $N_AandN_B$ .What's more,the more steps we do, the more precise the results are. Because the high order terms we dropped are getting smaller as the time step go smaller.

● Conclusion

Euler mothed is powerful to solve simple first order ordinary differential equations.draw figures by Python is a little difficult as well as interesting.It is essential to master the basic Python syntax，for me,there is still a long way to go.

addition

In order to celebrate the national day, we use Python to draw a five star red flag.
Codes are as follows:

#coding=utf-8
'''
Created on 2016-10-1
@author: 胡胜勇
'''
import turtle
import math
def draw_polygon(aTurtle, size=50, n=3):
    ''' 绘制正多边形
    args:
        aTurtle: turtle对象实例
        size: int类型，正多边形的边长
        n: int类型，是几边形        
    '''
    for i in xrange(n):
        aTurtle.forward(size)
        aTurtle.left(360.0/n)
def draw_n_angle(aTurtle, size=50, num=5, color=None):
    ''' 绘制正n角形，默认为黄色
    args:
        aTurtle: turtle对象实例
        size: int类型，正多角形的边长
        n: int类型，是几角形    
        color: str， 图形颜色，默认不填色
    '''
    if color:
        aTurtle.begin_fill()
        aTurtle.fillcolor(color)
    for i in xrange(num):
        aTurtle.forward(size)
        aTurtle.left(360.0/num)
        aTurtle.forward(size)
        aTurtle.right(2*360.0/num)
    if color:
        aTurtle.end_fill()
def draw_5_angle(aTurtle=None, start_pos=(0,0), end_pos=(0,10), radius=100, color=None):
    ''' 根据起始位置、结束位置和外接圆半径画五角星
    args:
        aTurtle: turtle对象实例
        start_pos: int的二元tuple，要画的五角星的外接圆圆心
        end_pos: int的二元tuple，圆心指向的位置坐标点
        radius: 五角星外接圆半径
        color: str， 图形颜色，默认不填色    
    '''
    aTurtle = aTurtle or turtle.Turtle()
    size = radius * math.sin(math.pi/5)/math.sin(math.pi*2/5)
    aTurtle.left(math.degrees(math.atan2(end_pos[1]-start_pos[1], end_pos[0]-start_pos[0])))
    aTurtle.penup()
    aTurtle.goto(start_pos)
    aTurtle.fd(radius)
    aTurtle.pendown()
    aTurtle.right(math.degrees(math.pi*9/10))
    draw_n_angle(aTurtle, size, 5, color)
def draw_5_star_flag(times=20.0):
    ''' 绘制五星红旗
    args:
        times: 五星红旗的规格为30*20， times为倍数，默认大小为10倍， 即300*200
    '''
    width, height = 30*times, 20*times
    # 初始化屏幕和海龟
    window = turtle.Screen()
    aTurtle = turtle.Turtle()
    aTurtle.hideturtle()
    aTurtle.speed(10)
    # 画红旗
    aTurtle.penup()
    aTurtle.goto(-width/2, height/2)
    aTurtle.pendown()
    aTurtle.begin_fill()
    aTurtle.fillcolor('red')
    aTurtle.fd(width)
    aTurtle.right(90)
    aTurtle.fd(height)
    aTurtle.right(90)
    aTurtle.fd(width)
    aTurtle.right(90)
    aTurtle.fd(height)
    aTurtle.right(90)    
    aTurtle.end_fill()
    # 画大星星
    draw_5_angle(aTurtle, start_pos=(-10*times, 5*times), end_pos=(-10*times, 8*times), radius=3*times, color='yellow')  
    # 画四个小星星
    stars_start_pos = [(-5, 8), (-3, 6), (-3, 3), (-5, 1)]
    for pos in stars_start_pos:
        draw_5_angle(aTurtle, start_pos=(pos[0]*times, pos[1]*times), end_pos=(-10*times, 5*times), radius=1*times, color='yellow')  
    # 点击关闭窗口
    window.exitonclick()
if __name__ == '__main__':
    draw_5_star_flag()

●result:###：用Python五星红旗.gif

Acknowledgments

Thanks to Yang Zongmeng,Luo Jiajia,Zhihu(知乎)