Download Radioactive Decay-Modeling and Simulation-Lecture Handouts and more Lecture notes Mathematical Modeling and Simulation in PDF only on Docsity!
Handout#
Radioactive Decay Some nuclei are energetically unstable and can spontaneously transform into more stable forms by various processes known collectively as radioactive decay. The rate at which a particular radioactive
sample will decay depends on the identity of the sample. The number of atoms in a sample that decay depends on the total number of atoms in the sample. This fact yields a rate of decay called an
exponential decay.
Example:
In 1986 the Chernobyl nuclear power plant exploded, and scattered radioactive material over Europe.
The radioactive released contained two radioactive elements iodine‐ 131 (I 131 ) whose half‐life is 8 days and cesium (Cs 137 ) whose half‐life is 30 years. Develop a mathematical model for the radioactive decay and show how much of it will remain over time.
Mathematical Model:
The number of radioactive atoms decaying over a period of time is a dependent quantity. The time is an independent variable.
A quantity is said to be subject to exponential decay if it decreases at a rate proportional to its value. Therefore, if N(t) denotes the amount of a radioactive substance present at time t , then (the rate
dN(t) / dt is negative, since N is decreasing. The positive constant k is called the rate constant for the particular radioisotope. The equation that describes exponential decay is
Nt
dt
dNt
= − λ ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐(I)
Solution to the differential equation: By rearranging,
dt
Nt
dNt
Integrating on both side
ln N ( t )=−λ t + C
where C is the constant of integration
ln N ( t )=−λ t + C
t C C t
N t e e e
− λ + − λ
t
N t Ne
− λ
( )= 0 ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐(II)
where the final substitution, N 0 = eC , is obtained by evaluating the equation at t = 0, as N 0 is defined as being the quantity at t = 0.
Derivation of Half‐life:
The number of radioactive remaining at the end of one half‐life is
( )^0
y
y t =
The value of t that satisfies this is the half‐life.
y = ye − kt
0
0
= e −^ kt
Take reciprocal of both sides
kt
= e − kt = e
Take natural log on both sides
ln 2 = ln( e kt^ )= kt
The half‐life is given by
k
ln( 2 )
General Solution with MATLAB
1. Plot the graph of the exponential decay function e‐kx^ and find the effect of k on the function for
k=0.2, 0.1 and 0. >> fplot('[exp(-2x),exp(-x),exp(-0.5x)]',[0,5]) xlabel('x'); ylabel('y'); title('Plot of Exponential Decay Functions'); legend('y = exp(-2x)','y = exp(-x)','y = exp(-0.5x)',0);
The larger the absolute value of the rate constant k, the faster the function decays.
2. Develop a mathematical model for the radioactive decay of Iodine and Cesium and show how much of it will remain over time. The decay constant for Iodine is given by
k = ln2/8 = 0.693/8 =0.0866 per day If t is measure in days, the amount of Iodine left at time t would be
NI(t) = N 0 e‐0.866t The decay constant for Iodine is given by
k = ln2/30 = 0.693/30 =0.023 per year
(^00) 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
1
x
y
Plot of Exponential Decay Functions
y = exp(-2x) y = exp(-x) y = exp(-0.5x)
Let us find the slopes to show that the data is not linear >> change_in_y = diff(data(:,2)); change_in_x = diff(data(:,1));slopes = change_in_y./change_in_x
slopes = -0. -0. -0. -0.
The slope is different at every point, which means that we cannot use a linear function. Let's try to use
the exponential function to model our data. We calculate the ratio as:
k x x
y y x x
y y = = −
−
−
− ....
log( ) log( ) log( ) log( ) 2 1
2 1 1 0
1 0
>> change_in_log_y = diff(log(data(:,2))); ratios = change_in_log_y./change_in_x
ratios = -0. -0. -0. -0.
The ratios are very close to each other so we can use an exponential function to model our data. First we need to find the mean of our ratios. >> k = mean(ratios) k = -0.
We can use the equation for exponential decay to model the data N=N 0 e‐kt. Let N 0 =12 is the initial mass of I‐ 131 (on day 0), k is decay ratio and t is time (in days). Thus we will use the function N=12e‐0.0866t^ to find the mass of radioactive iodine at any given day t. n0=12; t = 0:0.01:20;n = n0.exp(k.t); plot(t,n,'b-',data(:,1),data(:,2),'ro') xlabel('Time (days)');ylabel('Mass of I-131 (μg)'); title('Mass of I-131 versus time');legend('Model','Data',0)
(^20 2 4 6 8 10 12 14 16 18 )
3
4
5
6
7
8
9
10
11
12
Time (days)
Mass of I-131 ( g)
Mass of I-131 versus time Model Data
5. It is experimentally proven that the next dose of Iodine can be given when there are less than 5 μ g of the I‐ 131 left in patient's thyroid gland. How much of the I‐ 131 dose will still be in the patient’s
thyroid gland after 10 days. Use the graph of the model to find the approximate value.
We can use the graph of the model and find the approximate value. From the graph, we approximate
that M=5 μ g at a time t=10 days.
6. Using the exponential delay equation, write the MATLAB code to determine the exact time t at
which the mass M=5 μ g.
We need solve the equation N 0 e kt^ = N for t.
Divide by N 0 and take the log of both sides
0
log
N
N
kt
Divide through by k
k
N
N
t
log
Given k and N 0 from our model, if we know a given N , we can calculate the corresponding t value with the above equation.
>> n0=12; t = 0:0.01:20; n = 5; k=-0. t=log(n/n0)/k t =
The next iodine treatment for the patient should be taken 10 days after the first injection.
7. Calculate the half‐life of I‐131.
The half life is the day at which only half of the initial amount remains (N=1/2 N 0 ). We can use the equation for radioactive decay
1 / 2 1 / 2 0 0
N ( t )= n = Nekt
Canceling out N 0 on both sides
1 / 2
1 kt
= e
Take log of both sides and solve for half life t1/
k
t
log( 1 / 2 )
>> t_hl=log(1/2)/k t_hl =
We calculated that the half life of I‐ 131 is 8 days.
