Chapter 1 Problem 1.5

---

Name: Feng Chen
Student number： 2013301020145
March 19, 2016

Abstract

This is a programme to solve a specific first order ODE system via Euler method and store the datas in a txt file. Finally, I plotthe datas and compare this result to the analytical result. Also, I write a programme to read the datas in the .txt files to plot.

Introduction

This is a decay problem with two types of nuclei A ans B, which can decay into each other with the same time constant τ. The corresponding rate equations are:

$$\begin{equation} \frac{dN_A}{dt}=\frac{N_B}{\tau}-\frac{N_A}{\tau} \label{} \end{equation}$$
$$\begin{equation} \frac{dN_B}{dt}=\frac{N_A}{\tau}-\frac{N_B}{\tau} \label{} \end{equation}$$
Obviously, the number of these two types ($N_A$ and $N_B$) of particles will reach a steady state. In such a steady state, the time derivative should vanish.

Nummerical Approach -- Euler method

Firstly, we do the Taylor expansion for the particle numbers N at time t :

$$\begin{equation} N（t+\Delta t）=N(t)+\frac{dN(t)}{dt}\Delta t+\frac{1}{2}\frac{d^2N(t)}{dt^2}(\Delta t)^2+... \approx N(t)+\frac{dN(t)}{dt}\Delta t \label{} \end{equation}$$
When Δt is small enough, the error terms can be negligible.
Thus, the decay rate equations of particle A and B becomes:
$$\begin{equation} N_A(t+\Delta t)=N_A(t)+\frac{N_B(t)}{\tau}\Delta t-\frac{N_A(t)}{\tau}\Delta t \label{} \end{equation}$$ $$\begin{equation} N_B(t+\Delta t)=N_B(t)+\frac{N_A(t)}{\tau}\Delta t-\frac{N_B(t)}{\tau}\Delta t \label{} \end{equation}$$

Analytical Approcah

Firstly, we combine two rate equations:

$$\begin{equation} \tau\frac{d^2N_A}{dt^2}+2\frac{dN_A}{dt}=0 \label{} \end{equation}$$
Solving this equation, we get the function for particle A:
$$\begin{equation} N_A(t) = Ce^{-\frac{2t}{\tau}}+D \label{} \end{equation}$$ $$\begin{equation} N_B(t) = -Ce^{-\frac{2t}{\tau}}+D \label{} \end{equation}$$
where C and D are two constant that determine by the initial conditions. If the initial number of particles are given, we can get the complete expressions.
$$\begin{equation} N_A(t) = \frac{N_A(0)-N_B(0)}{2}e^{-\frac{2t}{\tau}}+\frac{N_A(0)+N_B(0)}{2} \label{} \end{equation}$$ $$\begin{equation} N_B(t) = -\frac{N_A(0)-N_B(0)}{2}e^{-\frac{2t}{\tau}}+\frac{N_A(0)+N_B(0)}{2} \end{equation}$$
In general, if the time constant of A is different with that of B, the results become:
$$\begin{equation} N_A(t) = \frac{N_A(0)\alpha-N_B(0)}{\alpha+1}e^{-\frac{(\alpha+1)t}{\tau_B}}+\frac{N_A(0)+N_B(0)}{\alpha+1} \label{} \end{equation}$$ $$\begin{equation} N_B(t) = -\frac{N_A(0)\alpha-N_B(0)}{\alpha+1}e^{-\frac{(\alpha+1)t}{\tau_B}}+\frac{N_A(0)+N_B(0)}{\alpha+1} \label{} \end{equation}$$
where α is the ratio of time constant B and time constant A:
$$\begin{equation} \alpha = \frac{\tau_B}{\tau_A} \label{} \end{equation}$$

Data analysis

Setting the initial number of A is 100 and that of B is zero and the same time constant 1s.

As the figures showed, when we make the time step smaller, the numerical result will be much more accurate. Also, due to the same time constant, the final number of both particles will be equal.

In general case, I set the different time constant for particel A and particle B. Sstting the initial number of A is 1000 and that of B is zero.

In this case the ratio of final number of B and A is equal to their time constant ratio. (For the particel number will reach the constant term in their expressions as time grows up.)

Finally, I add particle B at the initial time as well.

Read datas in the .txt file

In the previous programme, all the datas will be recorded in a .txt file and the format of datas are as follows:

As showed in the figure, the datas corresponding to the same time t are in the same row and each of them are seperated by a blank space.
Thus, we should transfer this data list into a 2-dimension list and each column is what we want.
The core codes are as follows:

def read_initial(initial_file):
    itxt= open(initial_file)
    initial_data=[]
    data =[]  
    for lines in itxt.readlines():
        lines=lines.replace("\n","").split(",")
        initial_data.append(lines)
    for i in range(len(initial_data[0][0].split(' '))):
        data.append([])
        for j in range(len(initial_data)):
            data[i].append(initial_data[j][0].split(' ')[i])
    itxt.close()
    return data     ### data[0] = x, data[1] = y1, data[2] = y2,...

After transferring the data into a 2D list, we can plot the number of particles versus time:

Conclusion

In this project, we can learn how to use Euler method to solve simple first order ordinary differential equation or equations. On the other hand, all the high order differential equation can be transferred into equation systems, which means Euler method is also a way to solve high order differential equation.

Programme codes

• Programme for problem 1.5
• Programme for data reading