- Email: [email protected]
Error analysis of large-time approximations for decay chains
Applied Mathematics and Computation 218 (2011) 3311–3319
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Error analysis of large-time approximations for decay chains Frank Massey ⇑, Jeffrey Prentis University of Michigan-Dearborn, Dearborn, MI 48128, United States
a r t i c l e
i n f o
a b s t r a c t Consider a radioactive decay chain X1 ? ? Xn? and let Nn(t) be the amount of Xn at time t. This paper establishes error bounds for large-time approximations to Nn(t) that include and generalize the transient equilibrium approximations and other known approximations. The error bounds allow one to find the range of t for which these approximations can be used with a given degree of precision. Ó 2011 Elsevier Inc. All rights reserved.
Keywords: Radioactive decay chain Transient equilibrium Bateman’s formula Large-time approximation Error bound Exponential convolution Reduced chain
1. Introduction Consider a radioactive decay chain X1 ? X2 ? ? Xn ? where a fraction bj of the nuclide Xj decays into Xj+1 with decay constant kj and half-life T = ln(2)/kj. The remaining fraction 1 bj of Xj decays with the same decay constant kj into other nuclides outside the chain. If Xn is stable then kn = 0. Let Nj(t) be that portion of the amount of Xj present at time t that has been produced by decays following the chain. The Nj(t) satisfy the radioactive decay equations:
dN 1 =dt ¼ k1 N1 ; dN j =dt ¼ bj1 kj1 Nj1 kj Nj
ð1Þ
for j P 2;
see Segrè [1, p. 172]. If Nj(0) = 0 for j P 2 then Bateman’s formula [2] for Nn(t) is
Nn ðtÞ ¼ N1 ð0Þb1;n k1;n En ðt; k1 ; . . . ; kn Þ;
ð2Þ
where
En ðt; k1 ; . . . ; kn Þ ¼
n X
C j ekj t ;
ð3Þ
j¼1
bi;j ¼ bi biþ1 bj1
and ki;j ¼ ki kiþ1 kj1
and C j ¼
n Y
ðki kj Þ1
i¼1 i–j
and bi,i = ki,i = 1 and E1(t ; k) = ekt. This formula is only defined if all the kj are distinct. Formulas for En(t ; k1, . . . , kn) when some of the kj are equal are more complicated; see Mathai [3, Section 2] and Cetnar [4, Section 3]. There are several useful large-time approximations for Nn(t). These involve the smallest of k1, . . . , kn which we denote by kp. Note that in (3) the ekj t for j – p go to zero faster than ekp t as t ? 1 which implies
⇑ Corresponding author. Address: Department of Mathematics & Statistics, University of Michigan-Dearborn, Dearborn, MI 48128, United States. E-mail address: [email protected] (F. Massey). 0096-3003/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.08.072
3312
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
Nn ðtÞ N1 ð0Þaekp t
ð4Þ kp t
for large t with a = b1,nk1,nCp. This approximation holds in the sense that N n ðtÞ=½N 1 ð0Þae ! 1 as t ? 1. The approximation (4) is widely used, see e.g. Rutherford et al. [5, p. 3] and Segrè [1, p. 174]. However, we could find no mention of how accurate it is. The following theorem gives an error bound for this approximation. Theorem 1. Let kp and kq be the smallest and 2nd smallest of k1, . . . , kn respectively and a = b1,nk1,nCp. Assume kp < kq and kq is strictly less than the 3rd smallest of k1, . . . , kn. Then
jNn ðtÞ N1 ð0Þaekp t j 6 ceðkq kp Þt N1 ð0Þaekp t
ð5Þ
Q where c ¼ nj¼1 ðkj kp Þ=ðkj kq Þ. Theorem j–p;q 1 is the special case of Theorem 3 below when m = 1 in that theorem. Note that it follows from (5) that the relative error in (4) will be less than e if t > ln(c/e)/(kq kp), which tells how large t must be for (4) to hold with a given degree e of precision. In particular, the relative error in (4) will be small if t is large compared to the second largest half-life unless the second largest half-life is close to the largest. See Example 1 in Section 3 for a typical application of Theorem 1. Products like the one defining c in Theorem 1 appear in a number of the inequalities below, so it is useful to make the following definition. Definition 1. Let ap and aq be the smallest and second smallest of a1, . . . , ak. If b is different from a1, . . . , ak and aq is strictly smaller than the third smallest of the aj let
/ða1 ; . . . ; ak jbÞ ¼
vða1 ; . . . ; ak Þ ¼
k Y
aj =ðaj bÞ;
j¼1 k Y
ðaj ap Þ=ðaj aq Þ:
j¼1 j–p;q
In the case k = 2 let v(a1, a2) = 1. Note that v(a1, . . . , ak) 1 if aq is substantially smaller than the third smallest of the aj which occurs frequently in decay chains. Also note that c = v(k1, . . . , kn) in Theorem 1. Another useful approximation is
Nn ðtÞ aNk ðtÞ
ð6Þ
for large t; see Segrè [1, p. 174]. Here p 6 k < n and a = (bk,nkk/kn)/(kk+1, kk+2, . . . , knjkp) and (6) holds in the sense that Nn(t)/ [aNk(t)] ? 1 as t ? 1. One way to see that this is valid is to apply (4) to both Nn(t) and Nk(t) giving Q N n ðtÞ=½N 1 ð0Þb1;n k1;n C p ekp t ! 1 and N k ðtÞ=½N 1 ð0Þb1;k k1;k C p;k ekp t ! 1 as t ? 1 where C p;k ¼ ki¼1 ðki kp Þ1 . Dividing gives i–p
Nn(t)/[aNk(t)] ? 1 as t ? 1. The following theorem gives two error bounds for (6). The first, (7), is similar to (5) while (8) shows that the approximation (6) is good for a larger range of values of t than (7). However (8) is not always applicable. Theorem 2. Let kp and kq be the smallest and 2nd smallest of k1, . . . , kn respectively, p 6 k < n and a = (bk,nkk/kn)/ (kk+1, kk+2, . . . , knjkp). For (7) assume kp < kq and kq is strictly less than the 3rd smallest of k1, . . . , kn and let c = v(k1, . . . , kn). For P (8) assume kq is strictly less than all of kk+1, . . . , kn and let W ¼ ðk 1Þvðkp ; kq ; kkþ1 ; . . . ; kn Þ nj¼kþ1 1=ðkj kp Þ. Then
jNn ðtÞ aNk ðtÞj 6 ceðkq kp Þt aNk ðtÞ;
ð7Þ
jNn ðtÞ aNk ðtÞj 6 aNk ðtÞW=t:
ð8Þ
Theorem 2 is the special case of Theorem 3 when m = k and r(j) = j for j = 1, . . . , n in that theorem. Note that it follows from (8) that if q < p then the relative error in (6) will be less than e if t > W/e. Again, this tells how large t must be for (6) to hold with a given degree of precision. If fact, if kk+1, . . . , kn are all substantially larger than kq then W will be of the same order of magnitude as the largest of the half-lives of Xk+1, . . . , Xn. In this case the relative error in (6) will be small if t is large compared to the half-lives of Xk+1, . . . , Xn. This will be a larger range of values of t than the range implied by (7). For a typical application, see Example 1 in Section 3. The activity, Aj(t) of Xj is the rate of decay of Xj, i.e. Aj(t) = kjNj(t). It follows from (6) that An(t)/Ap(t) bp,n[kp+1/ (kp+1 kp)] [kn/(kn kp)] for large t. In other words the activity of Xn is approximately proportional to that of Xp for large t, a situation that is termed transient equilibrium; see Segrè [1, p. 174]. A generalization of (4) and (6) is obtained by approximating the original chain by a reduced chain that is obtained by deleting some nuclides with small half-lives, usually the nuclides whose half-lives are less than a certain threshold. This type of approximation is also widely used; see e.g. Benedict et al. [6, p. 39], Ball and Adams [7], Bell [8] and Thomas and Barber [9]. To be precise, assume the following.
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
9 Let 1 6 m 6 n 1 and rð1Þ < rð2Þ < < rðmÞ be such that the nuclide > > > > > X p with the longest half-life is included in X rð1Þ ; . . . ; X rðmÞ and let > > > > > rðm þ 1Þ < rðm þ 2Þ < < rðnÞ be such that X rðmþ1Þ ; X rðmþ2Þ ; ; X rðnÞ > > > > > are the remaining nuclides: Treat the decays X rðmþ1Þ ! X rðmþ1Þþ1 ; > > > > = X rðmþ2Þ ! X rðmþ2Þþ1 ; . . . ; X rðnÞ ! X rðnÞþ1 as being instantaneous and consider the reduced chain X rð1Þ ! X rð2Þ ! ! X rðmÞ ! with li ¼ krðiÞ as the decay constant for X rðiÞ : In this reduced chain brðiÞ;rðjÞ is the fraction of X rðiÞ that decays into X rðjÞ : Let Mi ðtÞ be the amount of X rðiÞ present in the reduced chain at time t and
li;j ¼ li liþ1 lj1 :
> > > > > > > > > > > > > > > > > > > ;
3313
ð9Þ
Computations using the full chain are sometimes subject to round-off errors when the half-lives differ greatly in magnitude. The reduced chain decreases the variation in magnitude of the half-lives which can avoid this problem; see Bell [8] and Thomas and Barber [9]. Usually one takes Xr(1), . . . , Xr(m) to be the nuclides with the m longest half-lives, but other choices are useful. If we assume M1(0) = N1(0)b1,r(1) and Mi(0) = 0 for i P 2 then (2) and (3) imply
Mm ðtÞ ¼ N1 ð0Þb1;rðmÞ l1;m Em ðt; l1 ; . . . ; lm Þ:
ð10Þ
The approximations (4) and (6) generalize to
Nn ðtÞ aM m ðtÞ
ð11Þ
for large t where a = (br(m),nlm/kn)/(lm+1, lm+2, . . . , lnjkp). Note that br(m),nlm/kn = 1 if r(m) = n with the understanding that kn/kn = 1 if kn = 0. If kp << kr(i) for i > m then /(lm+1, lm+2, . . . , lnjkp) 1 and (br(m),nlm/kn)Mm(t) provides a good approximation to Nn(t) for large t. The following theorem provides error bounds for the approximation (11). Theorem 3. Let kp and kq be the smallest and 2nd smallest of k1, . . . , kn respectively and a = (br(m),nlm/kn)/(lm+1, lm+2, . . . , lnjkp). Let r(1) < r(2) < < r(m) and r(m + 1) < r(m + 2) < < r(n) be such that all the integers 1, . . . , n are included in r(1), . . . , r(n) and p is included in r(1), . . . , r(m). Let li = kr(i) and li,j = lili+1 lj1 and Mm(t) = N1(0)b1,r(m)l1,mEm(t ; l1, . . . , lm). For (13) assume kp < kq and kq is strictly less than the 3rd smallest of k1, . . . , kn and let c = v(k1, . . . , kn). For (14) assume q is included in r(1), . . . , r(m) and kq is strictly less than all of lm+1, . . . , ln and let
W ¼ ðm 1Þvðkp ; kq ; lmþ1 ; . . . ; ln Þ
n X
1=ðlj kp Þ:
ð12Þ
j¼mþ1
Then
jNn ðtÞ aMm ðtÞj 6 ceðkq kp Þt aM m ðtÞ;
ð13Þ
jNn ðtÞ aMm ðtÞj 6 aM m ðtÞW=t:
ð14Þ
The Proof of Theorem 3 is a fairly straightforward application of Theorems 5 and 8 in Section 2.2, so the proof is postponed until after Theorem 8. Example 2 in Section 3 illustrates the application of Theorem 3 to a specific chain. Section 2 establishes similar error bounds for related approximations when X1 is continually being produced from an outside source. These error bounds are also the basis for the Proof of Theorem 3. Section 2 also considers the case when more than one nuclide is initially present. Section 3 illustrates the error bounds with applications to specific decay chains. A future paper is devoted to approximations to Nn(t) for small and intermediate values of t. 2. Theory 2.1. Properties of En(t) Most of the theory is based on the following convolution formula. t
t
En ðt; k1 ; . . . ; kn Þ ¼ ek1 ekn
ð15Þ
for t P 0; see Hamawi [10]. This formula is valid even if some of the kj are equal. Here ⁄ denotes convolution, i.e. Rt gðtÞ hðtÞ ¼ 0 gðsÞhðt sÞds for t P 0. This is defined for g(t) and h(t) belonging to various classes of functions, e.g. if they are measurable and bounded on any finite interval. Convolution is commutative and associative, so En(t ; k1, . . . , kn) is a symmetric function of k1, . . . , kn, i.e. rearranging k1, . . . , kn in En(t ; k1, . . . , kn) does not change its value. We will use the following functions related to En(t). Definition 2. t
t
An ðt; k1 ; . . . ; kn Þ ¼ k1;nþ1 En ðt; k1 ; . . . ; kn Þ ¼ ðk1 ek1 Þ ðkn ekn Þ; n X k1;j Ej ðt; k1 ; . . . ; kj Þ: Sn ðt; k1 ; . . . ; kn Þ ¼ j¼1
3314
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
If one assumes that all the bj are one and at t = 0 there is one unit of X1 present and none of the other Xj then An(t ; k1, . . . , kn) = An(t), the rate of decay (or activity) of Xn at time t, and Sn(t ; k1, . . . , kn) is the total amount of the nuclides X1, X2, . . . , Xn present. An(t ; k1, . . . , kn) is also a symmetric function of k1, . . . , kn. An(t ; k1, . . . , kn) is the probability density function of a sum of independent exponential random variables and Sn(t ; k1, . . . , kn) is the corresponding survival (or complementary distribution) function; see Ross [11, p. 284]. Keilson and Sumita [12], Boland et al. [13], Bon and Pa˘lta˘nea [14], Pa˘lta˘nea [15], [16] and Mi et al. [17] have investigated relationships between An(t) = An(t ; k1, . . . , kn), Sn(t) = Sn(t ; k1, . . . , kn), Am(t) = Am(t ; l1, . . . , lm) and Sm(t) = Sm(t ; l1, . . . , lm) when the parameters l1, . . . , lm bear certain relationships to k1, . . . , kn (e.g. if m 6 n and kj 6 lj for j = 1, . . . , n). They show in various circumstances that (i) Sm(t) 6 Sn(t), (ii) An(t)/Sn(t) 6 Am(t)/Sm(t) and (iii) An(t)/Am(t) is non-decreasing. The relationships (i), (ii) and (iii) are increasingly stronger in the sense that (iii) implies (ii) implies (i). The relationship (i) for the case when l1, . . . , lm is a subset of k1, . . . , kn is related to Theorem 5 (b) below in the case a = 0. Simonsen [18] showed that (1)n1En(t ; k1, . . . , kn) is the (n 1)st divided difference with respect to k of the function ekt. It is also possible to express Nn(t) in terms of the exponential moment functions Mk(k1, . . . , kn) defined by Minor [19] and Mathews et al. [20]; see Harr [21, p. 13]. For other derivations of (2) see Erofeev [22], Pressyanov [23] and Moral and Pacheco [24]. The following theorem contains some basic properties of En(t) and An(t) and uses the function / in Defintion 1. Parts (a) and (f) are known but we include their short proofs. For (a) see Kaplan [25, p. 340]. Part (f) expresses conservation of nuclei. Theorem 4. Assume a1, . . . , ak and b are real numbers and g1(t), . . ., gk(t) are functions that are measurable and bounded on any = max{a1 . . . , ak}, Ek(t) = Ek(t ; a1, . . . , ak), Ak(t) = Ak(t ; a1, . . . , ak) and Sk(t) = Sk(t; a1, . . . , ak). finite interval. Let a = min{a1, . . . , ak}, a For (c) assume b is different from all the aj. For (d) assume k P 2 and 1 6 j 6 k. Then the following are true for t > 0. (a) (b) (c) (d) (e) (f) (g)
[ebtg1(t)]⁄ ⁄[ebt gk(t)] = ebt[g1(t)⁄ ⁄gk(t)] ebtEk(t ; a1, . . . , ak) = Ek(t; a1 b, . . . , ak b) ebtaAt k(t ; a1, . . . , ak) = /(a1, . . . , akjb)Ak(t ; a1 b, . . . , ak b) dðe j Ek ðtÞÞ ¼ eaj t Ek1 ðt; a1 ; . . . ; aj1 ; ajþ1 ; . . . ; ak Þ dt ebtEk(t) is a non-decreasing function of t if b P a dðSk ðtÞÞ ¼ Ak ðtÞ dt tk1eat /(k 1)! 6 Ek(t ; a1, . . . , ak) 6 tk1eat/(k 1)!
Proof. (a) The case k = 2 follows immediately from the definition of convolution and the general case is proved by induction. (b) and (c) follow from (a), (15) and Definition 2. It suffices to prove (d) for j = k in which case (b) and (15) give Rt eak t Ek ðtÞ ¼ Ek ðt; a1 ak ; . . . ; ak1 ak ; 0Þ ¼ 0 Ek1 ðs; a1 ak ; . . . ; ak1 ak Þds from which one obtains 0 a t ðe k Ek ðtÞÞ ¼ Ek1 ðt; a1 ak ; . . . ; ak1 ak Þ. Using (b) again gives (d). If one chooses j so that aj = a in (d) then one obtains (eatEk(t))0 P 0 which implies (e) for b = a. Multiplying by e(ba)t proves the general case. To prove (f), note that (d) implies E0j ¼ Ej1 aj Ej where E0 = 0. Multiplying by a1 aj1, summing on j, and using Definition 2 proves (f). To prove (g) note that ; . . . ; a Þ 6 Ek ðt; a1 ; . . . ; ak Þ 6 Ek ðt; a; . . . ; aÞ since ea t 6 eaj t 6 eat for all j. (g) follows from this and the known formula Ek ðt; a Ek(t ; a, . . . , a) = tk1eat/(k 1)! which is easily proved by induction. h 2.2. Decay chains with an input source We consider large-time approximations to Nn(t) when new X1 is being produced at a rate f(t) from a source Y outside the chain so one has Y ? X1 ? X2 ? ? Xn? and the first equation in (1) becomes dN1/dt = f(t) k1N1. We assume all the Nj(0) are zero in which case Nn(t) = b1,nk1,nEn(t; k1, . . . , kn)⁄f(t). This formula where f(t) is constant in time is derived in Hultqvist [26, pp. 48–50]; see also Skrable et al. [27]. The general case can be proved using Laplace transforms. We devote most of our attention to the case where f(t) = Ek(t ; a1, . . . , ak) with min{a1, . . . , ak} < min{k1, . . . , kn}. In this case we show
Nn ðtÞ aEk ðt; a1 ; . . . ; ak Þ
ð16Þ
for large t where a = b1,n/(k1, . . . , knja)/kn and a = min{a1, . . . , ak}; see Theorem 9 at the end of this section. Thus Nn(t) is approximately a constant times Ek(t ; a1, . . . , ak) for large t. Often approximations that are simpler than (16) are possible using Theorem 8 instead; see Example 3 in Section 3. The following theorem is the basis for the proof of the error bound (13) and also for the error bound (20) in Theorem 9 below. Theorem 5. Assume k1, . . . , kn, a1, . . . , ak and a are real numbers, kp is the smallest of k1, . . . , kn, a < kp, h = (k1 a) (kn a), Ek(t) = Ek(t ; a1, . . . , ak) and En(t) = En(t ; k1, . . . , kn). For (a) assume g(t) P 0 for t P 0 and eatg(t) is non-decreasing for t P 0. For (b) and (d) assume a = min{a1, . . . , ak}. For (c) assume kp is strictly smaller than all the other kj and let c1 = v(k1, . . . , kn, a). For (d) let m be the second smallest of k1, . . . , kn, a1, . . . , ak and c = v(k1, . . . , kn, a1, . . . , ak) and assume a < m and m is strictly smaller than the third smallest of k1, . . . , kn, a1, . . . , ak.Then the following are true for t > 0.
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
(a) (b) (c) (d)
3315
hEn(t) ⁄ g(t) 6 g(t) hEn(t) ⁄ Ek(t) 6 Ek(t) eat ð1 c1 eðkp aÞt Þ 6 hEn ðtÞ eat 6 eat Ek(t)(1 ce(va)t) 6 hEn(t) ⁄ Ek(t) 6 Ek(t)
Proof. (a) First suppose a = 0 so that g(t) is non-decreasing and k1, . . . , kn are strictly positive. Let An(t) = An(t ; k1, . . . , kn) = k1,n+1 Rt Rt En(t) as in Definition 2. Then An ðtÞ gðtÞ ¼ 0 An ðsÞgðt sÞds 6 gðtÞ 0 An ðsÞds 6 gðtÞ. The integral on the right is less than 1 since An(t) is a probability density function, as was noted in Section 2.1. Thus (a) is true if a = 0. To prove it in general we use Theorem 4 (a) and (b) to write hEn(t)⁄g(t) = eatz(t) where z(t) = hEn(t ; k1 a, . . . , kn a)⁄(eat g(t)). Using (a) in the case a = 0 one obtains z(t) 6 eatg(t) from which (a) follows. (b) follows from (a) and Theorem 4 (e). The right inequalities in (c) and (d) follow from (b). To prove the left inequality in (c), first note that we may assume p = n since En(t) is a symmetric function of k1, . . . , kn. We first do the case a = 0 so that k1, . . . , kn are strictly positive. Let An(t) = k1,n+1En(t). Then Rt R1 An ðtÞ 1 ¼ 0 An ðsÞds ¼ 1 t An ðsÞds since An(t) is a probability density function. By part (b) one has
An ðsÞ ¼ k1;nþ1 En1 ðs; k1 ; . . . ; kn1 Þ ekn s 6 /ðk1 ; . . . ; kn1 jkn Þkn ekn s R1 R1 Integrating gives t An ðsÞds 6 /ðk1 ; . . . ; kn1 jkn Þekp t . Combining with An ðtÞ 1 ¼ 1 t An ðsÞds gives the left inequality in (c) when a = 0. Note that when a = 0 one has v(k1, . . . , kn, a) = /(k1, . . . , kn1jkn). To prove the left inequality in (c) for general a we use Theorem 4 (a) and (b) to write hEn(t)⁄eat = eatz(t) where z(t) = hEn(t ; k1 a, . . . , kn a)⁄1. Using the left inequality in (c) when a = 0 one has zðtÞ P 1 /ðk1 a; . . . ; kn1 ajkn aÞeðkn aÞt ¼ 1 c1 eðkn aÞt . Combining with hEn(t)⁄eat = eat z(t) proves the left inequality in (c) for general a. To prove the left inequality in (d) we may assume a = ak since Ek(t) is a symmetric function of the a’s. Using (15) one has En(t)⁄Ek(t) = En+k1(t)⁄eat where En+k1(t) = En+k1(t ; k1, . . . , kn, a1, . . . , ak1). By (c) one has hfEn+k1(t)⁄eat P eat(1 ce(va)t) where f = (a1 a) (ak1 a). By (c) one has fEk(t) = fEk1(t ; a1, . . . , ak1)⁄eat 6 eat. Combining the above we get hfEn(t)⁄ Ek(t) P f(1 ce(va)t)Ek(t) from which the left inequality in (d) follows. h Remark. Let An(t) = An(t ; k1, . . . , kn) = k1,n+1En(t) and / = /(k1, . . . , knja) and note that En(t)⁄Ek(t) = En+k(t) where En+k(t) = En+k(t ; k1, . . . , kn, a1, . . . , ak). Part (d) of Theorem 5 implies jEn+k(t) h1Ek(t)j 6 ce(va)th1Ek(t) and jAn(t)⁄Ek(t) / Ek(t)j 6 ce(va)t/Ek(t). So the relative error of the approximations En+k(t) h1Ek(t) and An(t)⁄Ek(t) /Ek(t) is no more than ce(va)t, which implies these approximations are good for large t. Similar interpretations can be made for (19) in Theorem 8 below. Theorem 8 below is the basis for the proof of the error bound (14) and also for the error bound (21) in Theorem 9 below. Theorem 7 establishes an analogous error bound when (16) is generalized by replacing Ek(t) by a more general class of functions. These theorems use the elasticity of a function g(t) defined by E(g(t)) = tg0 (t)/g(t) where prime denotes differentiation. We start out with the following result connecting convolution and elasticity. Theorem 6. Let Ek(t) = Ek(t ; a1, . . . , ak) where a1, . . . , ak are non-negative and let g(t) and h(t) be continuously differentiable functions satisfying g(t) > 0 and h(t) > 0. For part (a) assume E(g(t)) 6 r and E(h(t)) 6 w for t > 0 for some numbers r and w, while for part (b) assume E(g(t)) P r and E(h(t)) P w for t > 0. Then the following are true for t > 0. (a) E(g(t) ⁄ h(t)) 6 r + w + 1 (b) E(g(t) ⁄ h(t)) P r + w + 1 (c) E(Ek(t)) 6 k 1. Proof. To prove (a) let u(t) = g(t)⁄h(t). It is not hard to see that uðtÞ ¼ t
tu0 ðtÞ ¼ t
Z
1
gðtsÞhðtð1 sÞÞds þ t
0
Z
1
tsg 0 ðtsÞhðtð1 sÞÞds þ t
0
Z
1
R1 0
gðtsÞhðtð1 sÞÞ so 0
gðtsÞtð1 sÞh ðtð1 sÞÞds:
0
Using the hypotheses one has
tu0 ðtÞ 6 uðtÞ þ rt
Z 0
1
gðtsÞhðtð1 sÞÞds þ wt
Z
1
gðtsÞhðtð1 sÞÞds ¼ ðr þ w þ 1ÞuðtÞ: 0
which proves (a). The proof of (b) is the same, but with the inequalities reversed. (c) follows from (a) since E(eat) 6 0 if a P 0. h The following theorem establishes an error bound for generalizations of (16) when Ek(t) is replaced by a more general class of functions. Theorem 7. Assume k1, . . . , kn are real numbers, En(t) = En(t ; k1, . . . , kn) and f(t) satisfies f(0) = 0, f(t) > 0 for t > 0, and there are at numbers a, b and r such that a < min{k1, . . . , kn}, a + b < min{k1, . . . , kn}, E(eatf(t)) 6 r for t > 0, and both eatf(t) and ebt dðe dtf ðtÞÞ are Pn non-decreasing for t P 0. Let h = (k1 a) (kn a) and W ¼ r/ðk1 a; . . . ; kn a j bÞ j¼1 1=ðkj aÞ. Then for t > 0 one has
3316
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
f ðtÞð1 W=tÞ 6 hEn ðtÞ f ðtÞ 6 f ðtÞ
ð17Þ
Proof. The right inequality in (17) follows from Theorem 5 (a), so it remains to prove the left inequality. We first consider the case a = 0, In this case the assumptions become (i) f(0) = 0, f(t) > 0 for t > 0 and f(t) is non-decreasing for t P 0, (ii) E(f(t)) 6 r for t > 0, (iii) k1, . . . , kn are all strictly positive, (iv) b < min{k1, . . . , kn} and (v) ebtf0 (t) is non-decreasing for t P 0. In this case we need to show
f ðtÞð1 r/ðk1 ; . . . ; kn jbÞU=tÞ 6 An ðtÞ f ðtÞ
ð18Þ
Pn
where An(t) = k1,n+1En(t) and U ¼ j¼1 1=kj . Let Sn(t) = Sn(t; k1, . . . , kn). Integrating by parts and using Theorem 4 (f), f(0) = 0 and Sn(0) = 1 one obtains An(t)⁄f(t) = f(t) Sn(t)⁄f0 (t). Theorem 5 (a) gives Aj(t ; k1, . . . , kj)⁄f0 (t) 6 /(k1, . . . , kjjb)f 0 (t). Dividing by kj, summing on j and using the fact that /(k1, . . . , kjjb) 6 /(k1, . . . , knjb) one obtains Sn(t)⁄f0 (t) 6 /(k1, . . . , knjb)Uf 0 (t). Using the fact that E(f(t)) 6 r one obtains Sn(t)⁄f0 (t) 6 r/(k1, . . . , knjb)Uf(t)/t. Combining with the previous gives (18). To prove the left inequality in (17) when a is not necessarily zero we use Theorem 4 (a) and (b) to obtain hEn(t)⁄f(t) = hat [En(t ; k1 a, . . . ,km a)⁄(eatf(t)]. Applying (18) to En(t ; k1 a, . . . , kn a)⁄(eatf(t)) one has n(t)⁄f(t) = he
eat f ðtÞð1 W=tÞ 6 hEn ðt; k1 a; . . . ; kn aÞ ðeat f ðtÞÞ Multiplying by eat and using the above formula for hEn(t)⁄f(t) gives the left inequality in (17). h The following theorem is the basis for the proofs of the error bounds (14) and (21). It can also be used for other similar approximations; Example 3 in Section 3 has a typical instance. Theorem 8. Assume k1, . . . , kn and a1, . . . , ak are real numbers and let En(t) = En(t ; k1, . . . , kn) and Ek(t) = Ek(t ; a1, . . . , ak). Assume k P 2 and let a and m be the smallest and second smallest of a1, . . . , ak. Assume m < min{k1, . . . , kn} and let h = (k1 a) (kn a) P and W ¼ ðk 1Þvða; m; k1 ; . . . ; kn Þ nj¼1 1=ðkj aÞ. Then for t > 0 one has
Ek ðtÞð1 W=tÞ 6 hEn ðtÞ Ek ðtÞ 6 Ek ðtÞ
ð19Þ
Proof. We may assume a = a1 and m = a2 since Ek(t) is a symmetric function of the a’s. Apply Theorem 7 to f(t) = Ek(t) with a = a1 and b = a2 a1. By Theorem 4 (b) and Theorem 6 (b) one has E(eatEk(t)) = E(Ek(t;0,a2 a, . . . , ak a)) 6 k 1 for t > 0. By Theorem 4 (b) and (d) one has ebt dðeat Ek ðtÞÞ=dt ¼ eða2 aÞt Ek1 ðt; a2 a; . . . ; ak aÞ. By Theorem 4 (e) one sees that eatEk(t) and ebtd(eatEk(t))/dt are non-decreasing. Therefore, all the hypotheses of Theorem 7 are satisfied and (19) follows from Theorem 7. h Here is the Proof of Theorem 3 which uses Theorems 5 and 8. Proof of Theorem 3. Let Em(t) = Em(t ; l1, . . . , lm), Enm(t) = Enm(t ; lm+1, . . . , ln) and En(t) = En(t ; k1, . . . , kn). Since l1, . . . , ln is a rearrangement of k1, . . . , kn one has En(t) = En(t; l1, . . . , ln) = Enm(t)⁄Em(t) and
k1;n ¼ k1;nþ1 =kn ¼ l1;nþ1 =kn ¼ l1;mþ1 /ðlmþ1 ; lmþ2 ; . . . ; ln jkp Þh=kn where h = (lm+1 kp) (ln kp). Theorem 5 (d) implies jhEnm ðtÞ Em ðtÞ Em ðtÞj 6 c2 eðkq kp Þt Em ðtÞ where c2 = v(l1, . . . , ln) and Theorem 8 implies jhEnm(t)⁄Em(t) Em(t)j 6 Em(t)W/t. Since l1, . . . , ln is a rearrangement of k1, . . . , kn one has v(l1, . . . , ln) = v(k1, . . . , kn). So c2 = c. We multiply the prceeding inequalities by N1(0)b1,nl1,m+1/(lm+1, lm+2, . . . , lnjkp)/kn and use the fact that (N1(0)b1,nl1,m+1/(lm+1,lm+2, . . . , lnjkp)/kn)Em(t) = aMm(t) and (N1(0)b1,nl1,m+1/(lm+1,lm+2, . . . , lnjkp)/kn) hEnm(t)⁄Em(t) = Nn(t) to get (13) and (14). h The following theorem has two error bounds for the approximation (16). (20) is similar to (13) while (21) is similar to (14). As with (13) and (14), the second error bound (21) shows that the relative error of the approximation (16) is small for a larger range of values of t than does the first, but is not applicable in some cases where (20) is. Theorem 9. Assume k1, . . . , kn are non-negative and a1, . . . , ak are arbitrary real numbers with a = min{a1, . . . , ak} and a < min{k1, . . . , kn}. Let Ek(t) = Ek(t; a1, . . . , ak) and Nn(t) = b1,nk1,nEn(t;k1, . . . , kn)⁄Ek(t) and a = b1,n/(k1, . . . , knja)/kn. For (20) let m be the second smallest of k1, . . . , kn, a1, . . . , ak, c = v(k1, . . . , kn, a1, . . . , ak) and assume a < m and m is strictly smaller than the third smallest of k1, . . . , kn, a1, . . . , ak. For (21) assume k P 2 and let m be the second smallest of a1, . . . , ak. Also assume P m < min{k1, . . . , kn} and let W ¼ ðk 1Þvða; m; k1 ; . . . ; kn Þ nj¼1 1=ðkj aÞ. Then for t P 0 one has
jNn ðtÞ aEk ðtÞj 6 ceðmaÞt aEk ðtÞ
ð20Þ
jNn ðtÞ aEk ðtÞj 6 aEk ðtÞW=t
ð21Þ
Proof. Let En(t) = En(t ; k1, . . . , kn). Theorem 5 (d) implies jhEn(t)⁄Ek(t) Ek(t)j 6 ce(ma)t Ek(t) where h = (k1 a) (kn a) and Theorem 8 implies jhEn(t)⁄Ek(t)Ek(t)j 6 Ek(t)W/t. Multiplying by a and using the fact that Nn(t) = haEn(t)⁄Ek(t) proves (20) and (21). h
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
3317
2.3. More than one initial nuclide In this section we consider the case where more than one Nj(0) is non-zero in which case (2) becomes
Nn ðtÞ ¼
n X
Nj;n ðtÞ;
ð22Þ
j¼1
Nj;n ðtÞ ¼ Nj ð0Þbj;n kj;n Enþ1j ðt; kj ; . . . ; kn Þ: We approximate each Nj,n(t) using the reduced chain (9). In addition to (9) we also assume k1, . . . , kn are all distinct and Xr(1), . . ., Xr(m) are the nuclides with the m largest half-lives, i.e. kr(i) < kr(k) for i 6 m and k > m. In fact, one can often choose m so kr(i) << kr(k) for i 6 m and k > m which makes applying the error bound (25) below simpler; see Example 4 in the next section. Using (11) a single term Nj,n(t) in (22) with j 6 r(m) can be approximated for large t by
Nj;n ðtÞ n/j Mj;m ðtÞ M j;m ðtÞ ¼ Nj ð0Þbj;rðmÞ lsðjÞ;m Emþ1sðjÞ ðt; lsðjÞ ; . . . ; lm Þ
ð23Þ
/j ¼ /ðlxðjÞ ; lxðjÞþ1 ; . . . ; ln jkpðjÞ Þ if j 6 rðnÞ and /j ¼ 1 if j > rðnÞ and n = br(m),nlm/kn. Here s(j) is such that Xr(s(j)) ? Xr(s(j)+1) ? ? Xr(m)? is the portion of the reduced chain that only includes nuclides from Xj, . . . , Xn, i.e. s(j) is the smallest i such that j 6 r(i). (s(j) is not defined if j > r(m).) Also x(j) is such that Xr(x(j)), Xr(x(j)+1), . . . , Xr(n) are the nuclides not in the reduced chain that are included in Xj, . . . , Xn, i.e. x(j) is the smallest i such that i > m and j 6 r(i). (x(j) is not defined if j > r(n).) p(j) is such that kp(j) is the smallest of kj, . . . , kn. If j >r(m) then Nj,n(t) 0 for t on the time scale of the reduced chain. Summing (23) one has
Nn ðtÞ nM m ðtÞ M m ðtÞ ¼
rX ðmÞ
ð24Þ
/j M j;m ðtÞ
j¼1
for large t. Note that Mm(t) is the amount of Xr(m) in the reduced chain at time t if the initial conditions are PrðiÞ Mi ð0Þ ¼ j¼rði1Þþ1 /j N j ð0Þbj;rðiÞ for i = 1, . . . , m where we set r(0) = 0. The following theorem gives an error bound for the approximation (24). Theorem 10. Assume the above and N1(0) > 0. For j < n let q(j) be such that kq(j) is the second smallest of kj, . . . , kn. If j 6 r(m 1) P W j ¼ ðm sðjÞÞvðkpðjÞ ; kqðjÞ ; lxðjÞ ; lxðjÞþ1 ; . . . ln Þ nk¼xðjÞ ½lk kpðjÞ 1 with Wj = 0 if j > r(n). Let W = max{Wj:1 6 j
let
6 r(m 1)} if m P 2 and W = 0 if m = 1. For r(m 1) < j 6 r(m) let cj = v(kj, . . . , kn) if j 6 n 1 with cn = 0 if r(m) = n. Let c = max{cj: r(m 1) < j 6 r(m)} and d = kq(r(m1)+1). Let g1 ðtÞ ¼ maxfW=t; ceðdlm Þt g. Let c = kp(r(m)+1) and P Q PrðmÞ l ¼ maxfl1 ; . . . ; lm g. Let V ¼ nj¼rðmÞþ1 N j ð0Þbj;n kpðjÞ wj where wj ¼ n i¼j ki =ðki kpðjÞ Þ and let RðtÞ ¼ j¼1 i–pðjÞ
/j N j ð0Þbj;n lsðjÞ;mþ1 tmsðjÞ =ðm sðjÞÞ! Let g2 ðtÞ ¼ VeðclÞt =RðtÞ if r(m) < n and g2(t) = 0 if r(m) = n. Then for t > 0 one has
jNn ðtÞ nM m ðtÞj 6 ½g1 ðtÞ þ g2 ðtÞnMm ðtÞ:
ð25Þ
Proof. Using the triangle inequality one has jNn(t) nMm(t)j 6 P1 + P2 where
P1 ¼
rX ðmÞ
jN j;n ðtÞ n/j M j;m ðtÞj
j¼1
and
P2 ¼
n X
N j;n ðtÞ:
j¼rðmÞþ1
By (14) one has
jNj;n ðtÞ n/j Mj;m ðtÞj 6 W j n/j M j;m ðtÞ=t 6 g1 ðtÞn/j Mj;m ðtÞ for 1 6 j 6 r(m 1). By (13) one has
jNj;n ðtÞ n/j Mj;m ðtÞj 6 cj eðkqðjÞ kpðjÞ Þt n/j Mj;m ðtÞ 6 g1 ðtÞn/j Mj;m ðtÞ for r(m 1) < j 6 r(m) since p(j) = r(m), kp(j) = kr(m) = lm and kq(j) P d. Summing these two previous inequalities it follows that P1 6 g1(t)nMm(t).
3318
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
For j > r(m) one can write kj;n Enþ1j ðt; kj ; . . . ; kn Þ ¼ ðkpðjÞ =kn Þwj hEnj ðtÞ ekpðjÞ t where Enj(t) = Enj(t ; kj, . . . , kp(j)1, Q kp(j)+1, . . . , kn) and h ¼ n i¼j ðki kpðjÞ Þ. By Theorem 5 (b) one has hEnj ðtÞ ekpðjÞ t 6 ekpðjÞ t 6 ect , so n X
P2 ¼
j¼rðmÞþ1
i–pðjÞ
Nj;n ðtÞ 6
n X
Nj ð0Þbj;n kpðjÞ wj ect =kn ¼ Vect =kn :
j¼rðmÞþ1
Theorem 4 (g) implies Emiþ1 ðt; li ; . . . ; lm Þ P t mi el t =ðm iÞ! for i = 1, . . . , m. Combining with (23) it follows that n/j M j;m ðtÞ P /j N j ð0Þbj;n lsðjÞ;mþ1 tmsðjÞ el t =½kn ðm sðjÞÞ! for j = 1, . . . , r(m). Therefore nM m ðtÞ P RðtÞel t =kn . So P2 6 Vect nM m ðtÞel t =RðtÞ ¼ g2 ðtÞnM m ðtÞ. Combining this with P1 6 g1(t)nMm(t) gives (25). h Remark. g1(t) comes from the error bounds (13) and (14) while g2(t) comes from the bounds in Theorem 4 (g) and Theorem 9(b) This results in g1(t) and g2(t) being a bit complicated. Both g1(t) and g2(t) approach zero as t goes to infinity, so (25) shows that the relative error in the approximation (24) is small for large t.
3. Examples We present four specific examples to illustrate the error bounds in Theorems 1, 2, 3, 8 and 10. This first example is an application of Theorems 1 and 2. Example 1. The chain 223Fr ? 223Ra ? 219Rn ? 215Po? has half-lives T1 = 21.8 min, T2 = 11.4 days, T3 = 3.96 s and T4 = 1.78 103 s and branching fractions b1 = 0.99994 and b2 = b3 = 1; see Firestone et al. [28]. In this case p = 2 and q = 1 in Theorems 1 and 2. Assume that the initial amounts of 223Ra, 219Rn and 215Po are zero. If one approximates the amount N4(t) of 215Po at time t using (4) and (6) with k = 2 then one obtains N 4 ðtÞ N 1 ð0Þb1;4 k1;4 C 2 ek2 t ¼ N 1 ð0Þb1 k1 k2 k3 ððk1 k2 Þðk3 k2 Þðk4 k2 ÞÞ1 ek2 t and N4(t) k2k3((k3 k2)(k4 k2))1N2(t ; k1, k2). If a 2% error is acceptable and one applies (5) and (7), then these approximations are good for t > t0 where t0 = ln(50c)/(k1 k2) where c = (k3 k2)(k4 k2)/(k3 k1)(k4 k1). However, c 1 and k1 k2 k1, so t0 ln(50)/k1 < 2.07 hours. On the other hand, if we apply (8), then N4(t) k2k3((k3 k2)(k4 k2))1N2(t ; k1, k2) with error less than 2% if t > t1 where t1 = 50 c[1/(k3 k2) + 1/ (k4 k2)] 50/k3 < 4.8 min. So the approximation N4(t) k2k3((k3 k2)(k4 k2))1N2(t ; k1, k2) has error less than 2% for t > 4.8 min which is better than the range given by (7). One final note, since b1 1 and kj/(kj k2) 1 for j = 1, 3, 4, the above approximations can be simplified to N 4 ðtÞ N 1 ð0Þðk2 =k4 Þek2 t and N4(t) (k2/k4)N2(t ; k1, k2) which have error less than 2% for t > 2.07 hours and t > 4.8 min respectively. The following example is a typical application of Theorem 3. Example 2. The chain 238U ? 234Th ? 234 mPa ? 234U? has half-lives T1 = 4.47 109 years, T2 = 24.1 days, T3 = 1.17 min and T4 = 2.46 105 years and branching fractions b1 = b2 = 1 and b3 = 0.9984; see Firestone et al. [28]. Here 234 m Pa denotes the metastable isomer of 234Pa. In this case p = 1. Assume that the initial amounts of 234Th, 234 m Pa and 234U are zero. Suppose one approximates the amount N4(t) of 234U at time t using (11) with the reduced chain 238U ? 234U ?. In the notation of (9) one has m = 2, r(1) = 1, r(2) = 4 and the deleted nuclides are 234Th and 234Pa with r(3) = 2 and r(4) = 3. Also, M1(0) = N1(0) and M2(t) = N1(0)b3k1E2(t ; k1, k4) and the approximation (11) is N4(t) /(k2, k3jk1)N1(0)b3k1E2(t ; k1, k4). If a 2% error is acceptable and one applies (14), then this approximation is good for t > t1 where t1 = 50v(k1, . . . , k4)[1/(k2 k1) + 1/ (k3 k1)] 50/k2 < 5 years. Also, /(k2,k3jk1)b3 1, so N4(t) N1(0)k1E2(t ; k1,k4) with error less than 2% for t > 5 years. This approximation is useful, not just for understanding the behavior of N4(t) for t > 5 years, but also because the calculation of N4(t) using (2) may encounter problems with overflow for large t. The following example is a typical application of Theorem 8. Example 3. Consider the series 92Kr ? 92Rb ? 92Sr ? 92 Y? with half-lives T1 = 3.0 s, T2 = 5.3 s, T3 = 2.71 h, T4 = 3.53 h and all bj = 1. Benedict et al. [6, p. 39] considers the case where 92Kr is produced at a constant rate P and the initial amounts of 92 Kr, 92Rb, 92Sr and 92Y are all zero. Let N4(t) be the amount of 92Y at time t. One has N4(t) = k1,4 En(t ; k1, . . . , k4)⁄P and it follows from (15) that N4(t) = Pk1k2k3E2(t ; k1,k2)⁄E3(t ; k3,k4,0). By Theorem 8, N 4 ðtÞ Pk3 E3 ðt; k3 ; k4 ; 0Þ ¼ Pðk4 Þ1 ½1 þ ðk4 k3 Þ1 ðk3 ek4 t k4 ek3 t Þ, an approximation that appears in Benedict et al. [6, p. 40] without an error bound. Note that in Theorem 8, a = 0, m = k4, h = k1k2, v(a,m,k1,k2) 1 and W = 2v(a,m,k1,k2)[1/(k1 a) + 1/(k2 a)] 2(T1 + T2)/(ln2). If an error of 2% error is acceptable then the approximation above is valid with error less than 2% for t > 50W 20 min. The following example is a typical application of Theorem 10. Example 4. Consider the chain 223Fr ? 223Ra ? 219Rn ? 215Po ? in Example 1. Suppose the initial amount of 223Fr, 223Ra, Rn and 215Po are N1(0) = 103 gm, N2(0) = 1.0 gm, N3(0) = 106 gm and N4(0) = 109 gm respectively. Suppose one approximates the amount N4(t) of 215Po at time t using (24) using the reduced chain 223Fr ? 223 Ra ?. In the notation of (9) one
219
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
3319
has m = 2, r(1) = 1, r(2) = 2 and the deleted nuclides are 219Rn and 215Po with r(3) = 3 and r(4) = 4. Also, M1(0) = /1N1(0), M2(0) = /2N2(0) where /1 = /(k3, k4jk1) and /2 = /(k3, k4jk2). So M 2 ðtÞ ¼ /1 N 1 ð0Þb1 k1 E2 ðt; k1 ; k2 Þ þ /2 N 2 ð0Þek2 t . One has n = b2,4k2/k4 = k2/k4 and the approximation (24) is N 4 ðtÞ ðk2 =k4 Þ½/1 N 1 ð0Þb1 k1 E2 ðt; k1 ; k2 Þ þ /N 2 ð0Þek2 t . Suppose an error of 2% error is acceptable and one applies Theorem 10. One has W = W1 = v(k1, . . . , k4)[1/(k3 k2) + 1/(k4 k2)] 1/k3. Also c = c2 = v(k2, . . . , k4) 1, d = k3 and d lm d. So ceðdlm Þt ek3 t and g1 ðtÞ maxf1=ðk3 tÞ; ek3 t g ¼ 1=ðk3 tÞ. The values of j in ¼ k1 ; c ¼ k3 and c l k3 . Also V = N3(0)b3kp(3)k4(k4 k3)1 + the range r(m) < j 6 n are j = 3 and 4. One has l kp(4)N4(0) = N3(0)k3k4(k4 k3)1 + k4N4(0) k3N3(0) + k4N4(0) and R(t) > /2N2(0)b2,4l2,3 k2N2(0). So g2 ðtÞ K ½k3 N 3 ð0Þþ k4 N 4 ð0Þek3 t =½k2 N 2 ð0Þ 0:8ek3 t . So g1 ðtÞ þ g2 ðtÞ 6 1=ðk3 tÞ þ 0:8ek3 t . So the error in the approximation above for N4(t) will be less than 2% if 1=ðk3 tÞ þ 0:8ek3 t < 0:02. Solving numerically, one sees that this occurs it t > 50/k3 4.8 min. Furthermore, b1 1, /1 1 and /2 1 so the above approximation can be simplified to N 4 ðtÞ ðk2 =k4 Þ½N 1 ð0Þk1 E2 ðt; k1 ; k2 Þ þ N 2 ð0Þek2 t . 4. Conclusions In this paper we have established error bounds for several frequently-used, large-time approximations for a radioactive decay chain X1? X2 ? ? Xn?. These approximations allow one to approximate the amount Nn(t) of Xn in the original chain by the amount of some nucleus in a simpler chain for sufficiently large times. The error bounds in this paper allow users of these approximations to determine how large time must be in order that the approximation hold with a given degree of precision. References [1] E. Segrè, Nuclei and Particles, second ed., Benjamin/Cummings, Reading, Massachusetts, 1977. [2] H. Bateman, Solution of a system of differential equations occurring in the theory of radioactive transformation, Proc. Cambridge. Philos. Soc. 15 (1910) 423–427. [3] A. Mathai, Storage capacity of dam with gamma type inputs, Ann. Inst. Statist. Math. 34 (Part A) (1982) 591–597. [4] J. Cetnar, General solution of Bateman equations for nuclear transmutations, Ann. Nucl. Energy 33 (2006) 640–645. [5] E. Rutherford, J. Chadwick, C.D. Ellis, Radiations from Radioactive Substances, Cambridge University Press, Cambridge, 1951. [6] M. Benedict, T.H. Pigford, H.W. Levi, Nuclear Chemical Engineering, second ed., McGraw-Hill, New York, 1981. [7] S.J. Ball, R.K. Adams, ‘‘MATEXP,’’ A General Purpose Digital Computer Program for Solving Ordinary Differential Equations by the Matrix Exponential Method, Oak Ridge National Laboratory, Oak Ridge, TN, 1967. [8] M.J. Bell, ORIGEN–The ORNL Isotope Generation and Depletion Code, Oak Ridge National Laboratory, Oak Ridge, TN, 1973. [9] G.F. Thomas, D.H. Barber, Numerical solution of Bateman systems, Ann. Nucl. Energy 20 (1993) 407–414. [10] J.M. Hamawi, A useful recurrence formula for the equations of radioactive decay, Nucl. Technol. 11 (1971) 84–88. [11] S.M. Ross, Introduction to Probability Models, eighth ed., Academic Press, Boston, 2003. [12] J. Keilson, U. Sumita, Uniform stochastic ordering and related inequalities, Canad. J. Stat. 10 (1982) 181–198. [13] P.J. Boland, E. El-Neweihi, F. Proschan, Schur properties of convolutions of exponential and geometric random variables, J. Multivariate Anal. 48 (1994) 157–167. [14] J.-L. Bon, E. Pa˘lta˘nea, Ordering properties of convolutions of exponential random variables, Lifetime Data Anal. 5 (1999) 185–192. [15] E. Pa˘lta˘nea, On the approximation of the convolutions of exponential type, in RoGer 2000-Brasßov. Proceedings of the 4thRomanian-German Seminar on Approximation Theory and its Applications (H. Gonska, D. Kacsó and L. Beutel, Eds.), Schriftenreihe des Fachbereichs Mathematik der GerhardMercator-Universität, No. 485, Gerhard-Mercator-Universität, Fachbereich 11/Mathematik, Duisburg, 2000, pp. 107-109. [16] E. Pa˘lta˘nea, On the comparison in hazard rate ordering of fail-safe systems, J. Stat. Planning Infer. 138 (2008) 1993–1997. [17] J. Mi, W. Shi, Y.Y. Zhao, Some properties of convolutions of Pascal and Erlang random variables, Stat. Probab. Lett. 78 (2008) 2378–2387. [18] W. Simonsen, On divided differences and osculatory interpolation, Skandinavisk Aktuarietidsskrift 31 (1948) 157–163. [19] B.M. Minor, Exponential characteristic spatial quadrature for discrete ordinates neutral particle transport in two-dimensional cartesian coordinates, Ph.D. Thesis, Air Force Institute of Technology, Dayton, OH, 1993. [20] K.A. Mathews, G. Sjoden, B. Minor, Exponential characteristic spatial quadrature for discrete ordinates radiation transport in slab geometry, Nuclear Sci. Eng. 118 (1994) 24–37. [21] L.J. Harr, Precise calculation of complex radioactive decay chains, Master’s Thesis, Air Force Institute of Technology, Dayton, OH, 2007. [22] B.V. Erofeev, Zhur. Fiz. Kim. 24 (1950) 721. [23] D.S. Pressyanov, Short solution of the radioactive decay equations, Am. J. Phys 70 (2002) 444–445. [24] L. Moral, A.F. Pacheco, Algebraic approach to the radioactive decay equations, Am. J. Phys 71 (2003) 684–686. [25] W. Kaplan, Operational Methods for Linear Systems, Addison-Wesley, Reading, Massachusets, 1962. [26] B. Hultqvist, Studies on Naturally Occurring Ionizing Radiations, Almqvist & Wiksells Boktrycheri AB, Stockholm, 1956. [27] K. Skrable, C. French, G. Chabot, A. Major, A general equation for the kinetics of linear first order phenomena and suggested applications, Health Phys. 27 (1974) 155–157. [28] R.B. Firestone, V.S. Shirley, C.M. Baglin, S.Y. Chu, J. Zipkin, Table of Isotopes, eighth ed., Wiley, New York, 1996.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Error analysis of large-time approximations for decay chains Frank Massey ⇑, Jeffrey Prentis University of Michigan-Dearborn, Dearborn, MI 48128, United States
a r t i c l e
i n f o
a b s t r a c t Consider a radioactive decay chain X1 ? ? Xn? and let Nn(t) be the amount of Xn at time t. This paper establishes error bounds for large-time approximations to Nn(t) that include and generalize the transient equilibrium approximations and other known approximations. The error bounds allow one to find the range of t for which these approximations can be used with a given degree of precision. Ó 2011 Elsevier Inc. All rights reserved.
Keywords: Radioactive decay chain Transient equilibrium Bateman’s formula Large-time approximation Error bound Exponential convolution Reduced chain
1. Introduction Consider a radioactive decay chain X1 ? X2 ? ? Xn ? where a fraction bj of the nuclide Xj decays into Xj+1 with decay constant kj and half-life T = ln(2)/kj. The remaining fraction 1 bj of Xj decays with the same decay constant kj into other nuclides outside the chain. If Xn is stable then kn = 0. Let Nj(t) be that portion of the amount of Xj present at time t that has been produced by decays following the chain. The Nj(t) satisfy the radioactive decay equations:
dN 1 =dt ¼ k1 N1 ; dN j =dt ¼ bj1 kj1 Nj1 kj Nj
ð1Þ
for j P 2;
see Segrè [1, p. 172]. If Nj(0) = 0 for j P 2 then Bateman’s formula [2] for Nn(t) is
Nn ðtÞ ¼ N1 ð0Þb1;n k1;n En ðt; k1 ; . . . ; kn Þ;
ð2Þ
where
En ðt; k1 ; . . . ; kn Þ ¼
n X
C j ekj t ;
ð3Þ
j¼1
bi;j ¼ bi biþ1 bj1
and ki;j ¼ ki kiþ1 kj1
and C j ¼
n Y
ðki kj Þ1
i¼1 i–j
and bi,i = ki,i = 1 and E1(t ; k) = ekt. This formula is only defined if all the kj are distinct. Formulas for En(t ; k1, . . . , kn) when some of the kj are equal are more complicated; see Mathai [3, Section 2] and Cetnar [4, Section 3]. There are several useful large-time approximations for Nn(t). These involve the smallest of k1, . . . , kn which we denote by kp. Note that in (3) the ekj t for j – p go to zero faster than ekp t as t ? 1 which implies
⇑ Corresponding author. Address: Department of Mathematics & Statistics, University of Michigan-Dearborn, Dearborn, MI 48128, United States. E-mail address: [email protected] (F. Massey). 0096-3003/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.08.072
3312
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
Nn ðtÞ N1 ð0Þaekp t
ð4Þ kp t
for large t with a = b1,nk1,nCp. This approximation holds in the sense that N n ðtÞ=½N 1 ð0Þae ! 1 as t ? 1. The approximation (4) is widely used, see e.g. Rutherford et al. [5, p. 3] and Segrè [1, p. 174]. However, we could find no mention of how accurate it is. The following theorem gives an error bound for this approximation. Theorem 1. Let kp and kq be the smallest and 2nd smallest of k1, . . . , kn respectively and a = b1,nk1,nCp. Assume kp < kq and kq is strictly less than the 3rd smallest of k1, . . . , kn. Then
jNn ðtÞ N1 ð0Þaekp t j 6 ceðkq kp Þt N1 ð0Þaekp t
ð5Þ
Q where c ¼ nj¼1 ðkj kp Þ=ðkj kq Þ. Theorem j–p;q 1 is the special case of Theorem 3 below when m = 1 in that theorem. Note that it follows from (5) that the relative error in (4) will be less than e if t > ln(c/e)/(kq kp), which tells how large t must be for (4) to hold with a given degree e of precision. In particular, the relative error in (4) will be small if t is large compared to the second largest half-life unless the second largest half-life is close to the largest. See Example 1 in Section 3 for a typical application of Theorem 1. Products like the one defining c in Theorem 1 appear in a number of the inequalities below, so it is useful to make the following definition. Definition 1. Let ap and aq be the smallest and second smallest of a1, . . . , ak. If b is different from a1, . . . , ak and aq is strictly smaller than the third smallest of the aj let
/ða1 ; . . . ; ak jbÞ ¼
vða1 ; . . . ; ak Þ ¼
k Y
aj =ðaj bÞ;
j¼1 k Y
ðaj ap Þ=ðaj aq Þ:
j¼1 j–p;q
In the case k = 2 let v(a1, a2) = 1. Note that v(a1, . . . , ak) 1 if aq is substantially smaller than the third smallest of the aj which occurs frequently in decay chains. Also note that c = v(k1, . . . , kn) in Theorem 1. Another useful approximation is
Nn ðtÞ aNk ðtÞ
ð6Þ
for large t; see Segrè [1, p. 174]. Here p 6 k < n and a = (bk,nkk/kn)/(kk+1, kk+2, . . . , knjkp) and (6) holds in the sense that Nn(t)/ [aNk(t)] ? 1 as t ? 1. One way to see that this is valid is to apply (4) to both Nn(t) and Nk(t) giving Q N n ðtÞ=½N 1 ð0Þb1;n k1;n C p ekp t ! 1 and N k ðtÞ=½N 1 ð0Þb1;k k1;k C p;k ekp t ! 1 as t ? 1 where C p;k ¼ ki¼1 ðki kp Þ1 . Dividing gives i–p
Nn(t)/[aNk(t)] ? 1 as t ? 1. The following theorem gives two error bounds for (6). The first, (7), is similar to (5) while (8) shows that the approximation (6) is good for a larger range of values of t than (7). However (8) is not always applicable. Theorem 2. Let kp and kq be the smallest and 2nd smallest of k1, . . . , kn respectively, p 6 k < n and a = (bk,nkk/kn)/ (kk+1, kk+2, . . . , knjkp). For (7) assume kp < kq and kq is strictly less than the 3rd smallest of k1, . . . , kn and let c = v(k1, . . . , kn). For P (8) assume kq is strictly less than all of kk+1, . . . , kn and let W ¼ ðk 1Þvðkp ; kq ; kkþ1 ; . . . ; kn Þ nj¼kþ1 1=ðkj kp Þ. Then
jNn ðtÞ aNk ðtÞj 6 ceðkq kp Þt aNk ðtÞ;
ð7Þ
jNn ðtÞ aNk ðtÞj 6 aNk ðtÞW=t:
ð8Þ
Theorem 2 is the special case of Theorem 3 when m = k and r(j) = j for j = 1, . . . , n in that theorem. Note that it follows from (8) that if q < p then the relative error in (6) will be less than e if t > W/e. Again, this tells how large t must be for (6) to hold with a given degree of precision. If fact, if kk+1, . . . , kn are all substantially larger than kq then W will be of the same order of magnitude as the largest of the half-lives of Xk+1, . . . , Xn. In this case the relative error in (6) will be small if t is large compared to the half-lives of Xk+1, . . . , Xn. This will be a larger range of values of t than the range implied by (7). For a typical application, see Example 1 in Section 3. The activity, Aj(t) of Xj is the rate of decay of Xj, i.e. Aj(t) = kjNj(t). It follows from (6) that An(t)/Ap(t) bp,n[kp+1/ (kp+1 kp)] [kn/(kn kp)] for large t. In other words the activity of Xn is approximately proportional to that of Xp for large t, a situation that is termed transient equilibrium; see Segrè [1, p. 174]. A generalization of (4) and (6) is obtained by approximating the original chain by a reduced chain that is obtained by deleting some nuclides with small half-lives, usually the nuclides whose half-lives are less than a certain threshold. This type of approximation is also widely used; see e.g. Benedict et al. [6, p. 39], Ball and Adams [7], Bell [8] and Thomas and Barber [9]. To be precise, assume the following.
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
9 Let 1 6 m 6 n 1 and rð1Þ < rð2Þ < < rðmÞ be such that the nuclide > > > > > X p with the longest half-life is included in X rð1Þ ; . . . ; X rðmÞ and let > > > > > rðm þ 1Þ < rðm þ 2Þ < < rðnÞ be such that X rðmþ1Þ ; X rðmþ2Þ ; ; X rðnÞ > > > > > are the remaining nuclides: Treat the decays X rðmþ1Þ ! X rðmþ1Þþ1 ; > > > > = X rðmþ2Þ ! X rðmþ2Þþ1 ; . . . ; X rðnÞ ! X rðnÞþ1 as being instantaneous and consider the reduced chain X rð1Þ ! X rð2Þ ! ! X rðmÞ ! with li ¼ krðiÞ as the decay constant for X rðiÞ : In this reduced chain brðiÞ;rðjÞ is the fraction of X rðiÞ that decays into X rðjÞ : Let Mi ðtÞ be the amount of X rðiÞ present in the reduced chain at time t and
li;j ¼ li liþ1 lj1 :
> > > > > > > > > > > > > > > > > > > ;
3313
ð9Þ
Computations using the full chain are sometimes subject to round-off errors when the half-lives differ greatly in magnitude. The reduced chain decreases the variation in magnitude of the half-lives which can avoid this problem; see Bell [8] and Thomas and Barber [9]. Usually one takes Xr(1), . . . , Xr(m) to be the nuclides with the m longest half-lives, but other choices are useful. If we assume M1(0) = N1(0)b1,r(1) and Mi(0) = 0 for i P 2 then (2) and (3) imply
Mm ðtÞ ¼ N1 ð0Þb1;rðmÞ l1;m Em ðt; l1 ; . . . ; lm Þ:
ð10Þ
The approximations (4) and (6) generalize to
Nn ðtÞ aM m ðtÞ
ð11Þ
for large t where a = (br(m),nlm/kn)/(lm+1, lm+2, . . . , lnjkp). Note that br(m),nlm/kn = 1 if r(m) = n with the understanding that kn/kn = 1 if kn = 0. If kp << kr(i) for i > m then /(lm+1, lm+2, . . . , lnjkp) 1 and (br(m),nlm/kn)Mm(t) provides a good approximation to Nn(t) for large t. The following theorem provides error bounds for the approximation (11). Theorem 3. Let kp and kq be the smallest and 2nd smallest of k1, . . . , kn respectively and a = (br(m),nlm/kn)/(lm+1, lm+2, . . . , lnjkp). Let r(1) < r(2) < < r(m) and r(m + 1) < r(m + 2) < < r(n) be such that all the integers 1, . . . , n are included in r(1), . . . , r(n) and p is included in r(1), . . . , r(m). Let li = kr(i) and li,j = lili+1 lj1 and Mm(t) = N1(0)b1,r(m)l1,mEm(t ; l1, . . . , lm). For (13) assume kp < kq and kq is strictly less than the 3rd smallest of k1, . . . , kn and let c = v(k1, . . . , kn). For (14) assume q is included in r(1), . . . , r(m) and kq is strictly less than all of lm+1, . . . , ln and let
W ¼ ðm 1Þvðkp ; kq ; lmþ1 ; . . . ; ln Þ
n X
1=ðlj kp Þ:
ð12Þ
j¼mþ1
Then
jNn ðtÞ aMm ðtÞj 6 ceðkq kp Þt aM m ðtÞ;
ð13Þ
jNn ðtÞ aMm ðtÞj 6 aM m ðtÞW=t:
ð14Þ
The Proof of Theorem 3 is a fairly straightforward application of Theorems 5 and 8 in Section 2.2, so the proof is postponed until after Theorem 8. Example 2 in Section 3 illustrates the application of Theorem 3 to a specific chain. Section 2 establishes similar error bounds for related approximations when X1 is continually being produced from an outside source. These error bounds are also the basis for the Proof of Theorem 3. Section 2 also considers the case when more than one nuclide is initially present. Section 3 illustrates the error bounds with applications to specific decay chains. A future paper is devoted to approximations to Nn(t) for small and intermediate values of t. 2. Theory 2.1. Properties of En(t) Most of the theory is based on the following convolution formula. t
t
En ðt; k1 ; . . . ; kn Þ ¼ ek1 ekn
ð15Þ
for t P 0; see Hamawi [10]. This formula is valid even if some of the kj are equal. Here ⁄ denotes convolution, i.e. Rt gðtÞ hðtÞ ¼ 0 gðsÞhðt sÞds for t P 0. This is defined for g(t) and h(t) belonging to various classes of functions, e.g. if they are measurable and bounded on any finite interval. Convolution is commutative and associative, so En(t ; k1, . . . , kn) is a symmetric function of k1, . . . , kn, i.e. rearranging k1, . . . , kn in En(t ; k1, . . . , kn) does not change its value. We will use the following functions related to En(t). Definition 2. t
t
An ðt; k1 ; . . . ; kn Þ ¼ k1;nþ1 En ðt; k1 ; . . . ; kn Þ ¼ ðk1 ek1 Þ ðkn ekn Þ; n X k1;j Ej ðt; k1 ; . . . ; kj Þ: Sn ðt; k1 ; . . . ; kn Þ ¼ j¼1
3314
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
If one assumes that all the bj are one and at t = 0 there is one unit of X1 present and none of the other Xj then An(t ; k1, . . . , kn) = An(t), the rate of decay (or activity) of Xn at time t, and Sn(t ; k1, . . . , kn) is the total amount of the nuclides X1, X2, . . . , Xn present. An(t ; k1, . . . , kn) is also a symmetric function of k1, . . . , kn. An(t ; k1, . . . , kn) is the probability density function of a sum of independent exponential random variables and Sn(t ; k1, . . . , kn) is the corresponding survival (or complementary distribution) function; see Ross [11, p. 284]. Keilson and Sumita [12], Boland et al. [13], Bon and Pa˘lta˘nea [14], Pa˘lta˘nea [15], [16] and Mi et al. [17] have investigated relationships between An(t) = An(t ; k1, . . . , kn), Sn(t) = Sn(t ; k1, . . . , kn), Am(t) = Am(t ; l1, . . . , lm) and Sm(t) = Sm(t ; l1, . . . , lm) when the parameters l1, . . . , lm bear certain relationships to k1, . . . , kn (e.g. if m 6 n and kj 6 lj for j = 1, . . . , n). They show in various circumstances that (i) Sm(t) 6 Sn(t), (ii) An(t)/Sn(t) 6 Am(t)/Sm(t) and (iii) An(t)/Am(t) is non-decreasing. The relationships (i), (ii) and (iii) are increasingly stronger in the sense that (iii) implies (ii) implies (i). The relationship (i) for the case when l1, . . . , lm is a subset of k1, . . . , kn is related to Theorem 5 (b) below in the case a = 0. Simonsen [18] showed that (1)n1En(t ; k1, . . . , kn) is the (n 1)st divided difference with respect to k of the function ekt. It is also possible to express Nn(t) in terms of the exponential moment functions Mk(k1, . . . , kn) defined by Minor [19] and Mathews et al. [20]; see Harr [21, p. 13]. For other derivations of (2) see Erofeev [22], Pressyanov [23] and Moral and Pacheco [24]. The following theorem contains some basic properties of En(t) and An(t) and uses the function / in Defintion 1. Parts (a) and (f) are known but we include their short proofs. For (a) see Kaplan [25, p. 340]. Part (f) expresses conservation of nuclei. Theorem 4. Assume a1, . . . , ak and b are real numbers and g1(t), . . ., gk(t) are functions that are measurable and bounded on any = max{a1 . . . , ak}, Ek(t) = Ek(t ; a1, . . . , ak), Ak(t) = Ak(t ; a1, . . . , ak) and Sk(t) = Sk(t; a1, . . . , ak). finite interval. Let a = min{a1, . . . , ak}, a For (c) assume b is different from all the aj. For (d) assume k P 2 and 1 6 j 6 k. Then the following are true for t > 0. (a) (b) (c) (d) (e) (f) (g)
[ebtg1(t)]⁄ ⁄[ebt gk(t)] = ebt[g1(t)⁄ ⁄gk(t)] ebtEk(t ; a1, . . . , ak) = Ek(t; a1 b, . . . , ak b) ebtaAt k(t ; a1, . . . , ak) = /(a1, . . . , akjb)Ak(t ; a1 b, . . . , ak b) dðe j Ek ðtÞÞ ¼ eaj t Ek1 ðt; a1 ; . . . ; aj1 ; ajþ1 ; . . . ; ak Þ dt ebtEk(t) is a non-decreasing function of t if b P a dðSk ðtÞÞ ¼ Ak ðtÞ dt tk1eat /(k 1)! 6 Ek(t ; a1, . . . , ak) 6 tk1eat/(k 1)!
Proof. (a) The case k = 2 follows immediately from the definition of convolution and the general case is proved by induction. (b) and (c) follow from (a), (15) and Definition 2. It suffices to prove (d) for j = k in which case (b) and (15) give Rt eak t Ek ðtÞ ¼ Ek ðt; a1 ak ; . . . ; ak1 ak ; 0Þ ¼ 0 Ek1 ðs; a1 ak ; . . . ; ak1 ak Þds from which one obtains 0 a t ðe k Ek ðtÞÞ ¼ Ek1 ðt; a1 ak ; . . . ; ak1 ak Þ. Using (b) again gives (d). If one chooses j so that aj = a in (d) then one obtains (eatEk(t))0 P 0 which implies (e) for b = a. Multiplying by e(ba)t proves the general case. To prove (f), note that (d) implies E0j ¼ Ej1 aj Ej where E0 = 0. Multiplying by a1 aj1, summing on j, and using Definition 2 proves (f). To prove (g) note that ; . . . ; a Þ 6 Ek ðt; a1 ; . . . ; ak Þ 6 Ek ðt; a; . . . ; aÞ since ea t 6 eaj t 6 eat for all j. (g) follows from this and the known formula Ek ðt; a Ek(t ; a, . . . , a) = tk1eat/(k 1)! which is easily proved by induction. h 2.2. Decay chains with an input source We consider large-time approximations to Nn(t) when new X1 is being produced at a rate f(t) from a source Y outside the chain so one has Y ? X1 ? X2 ? ? Xn? and the first equation in (1) becomes dN1/dt = f(t) k1N1. We assume all the Nj(0) are zero in which case Nn(t) = b1,nk1,nEn(t; k1, . . . , kn)⁄f(t). This formula where f(t) is constant in time is derived in Hultqvist [26, pp. 48–50]; see also Skrable et al. [27]. The general case can be proved using Laplace transforms. We devote most of our attention to the case where f(t) = Ek(t ; a1, . . . , ak) with min{a1, . . . , ak} < min{k1, . . . , kn}. In this case we show
Nn ðtÞ aEk ðt; a1 ; . . . ; ak Þ
ð16Þ
for large t where a = b1,n/(k1, . . . , knja)/kn and a = min{a1, . . . , ak}; see Theorem 9 at the end of this section. Thus Nn(t) is approximately a constant times Ek(t ; a1, . . . , ak) for large t. Often approximations that are simpler than (16) are possible using Theorem 8 instead; see Example 3 in Section 3. The following theorem is the basis for the proof of the error bound (13) and also for the error bound (20) in Theorem 9 below. Theorem 5. Assume k1, . . . , kn, a1, . . . , ak and a are real numbers, kp is the smallest of k1, . . . , kn, a < kp, h = (k1 a) (kn a), Ek(t) = Ek(t ; a1, . . . , ak) and En(t) = En(t ; k1, . . . , kn). For (a) assume g(t) P 0 for t P 0 and eatg(t) is non-decreasing for t P 0. For (b) and (d) assume a = min{a1, . . . , ak}. For (c) assume kp is strictly smaller than all the other kj and let c1 = v(k1, . . . , kn, a). For (d) let m be the second smallest of k1, . . . , kn, a1, . . . , ak and c = v(k1, . . . , kn, a1, . . . , ak) and assume a < m and m is strictly smaller than the third smallest of k1, . . . , kn, a1, . . . , ak.Then the following are true for t > 0.
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
(a) (b) (c) (d)
3315
hEn(t) ⁄ g(t) 6 g(t) hEn(t) ⁄ Ek(t) 6 Ek(t) eat ð1 c1 eðkp aÞt Þ 6 hEn ðtÞ eat 6 eat Ek(t)(1 ce(va)t) 6 hEn(t) ⁄ Ek(t) 6 Ek(t)
Proof. (a) First suppose a = 0 so that g(t) is non-decreasing and k1, . . . , kn are strictly positive. Let An(t) = An(t ; k1, . . . , kn) = k1,n+1 Rt Rt En(t) as in Definition 2. Then An ðtÞ gðtÞ ¼ 0 An ðsÞgðt sÞds 6 gðtÞ 0 An ðsÞds 6 gðtÞ. The integral on the right is less than 1 since An(t) is a probability density function, as was noted in Section 2.1. Thus (a) is true if a = 0. To prove it in general we use Theorem 4 (a) and (b) to write hEn(t)⁄g(t) = eatz(t) where z(t) = hEn(t ; k1 a, . . . , kn a)⁄(eat g(t)). Using (a) in the case a = 0 one obtains z(t) 6 eatg(t) from which (a) follows. (b) follows from (a) and Theorem 4 (e). The right inequalities in (c) and (d) follow from (b). To prove the left inequality in (c), first note that we may assume p = n since En(t) is a symmetric function of k1, . . . , kn. We first do the case a = 0 so that k1, . . . , kn are strictly positive. Let An(t) = k1,n+1En(t). Then Rt R1 An ðtÞ 1 ¼ 0 An ðsÞds ¼ 1 t An ðsÞds since An(t) is a probability density function. By part (b) one has
An ðsÞ ¼ k1;nþ1 En1 ðs; k1 ; . . . ; kn1 Þ ekn s 6 /ðk1 ; . . . ; kn1 jkn Þkn ekn s R1 R1 Integrating gives t An ðsÞds 6 /ðk1 ; . . . ; kn1 jkn Þekp t . Combining with An ðtÞ 1 ¼ 1 t An ðsÞds gives the left inequality in (c) when a = 0. Note that when a = 0 one has v(k1, . . . , kn, a) = /(k1, . . . , kn1jkn). To prove the left inequality in (c) for general a we use Theorem 4 (a) and (b) to write hEn(t)⁄eat = eatz(t) where z(t) = hEn(t ; k1 a, . . . , kn a)⁄1. Using the left inequality in (c) when a = 0 one has zðtÞ P 1 /ðk1 a; . . . ; kn1 ajkn aÞeðkn aÞt ¼ 1 c1 eðkn aÞt . Combining with hEn(t)⁄eat = eat z(t) proves the left inequality in (c) for general a. To prove the left inequality in (d) we may assume a = ak since Ek(t) is a symmetric function of the a’s. Using (15) one has En(t)⁄Ek(t) = En+k1(t)⁄eat where En+k1(t) = En+k1(t ; k1, . . . , kn, a1, . . . , ak1). By (c) one has hfEn+k1(t)⁄eat P eat(1 ce(va)t) where f = (a1 a) (ak1 a). By (c) one has fEk(t) = fEk1(t ; a1, . . . , ak1)⁄eat 6 eat. Combining the above we get hfEn(t)⁄ Ek(t) P f(1 ce(va)t)Ek(t) from which the left inequality in (d) follows. h Remark. Let An(t) = An(t ; k1, . . . , kn) = k1,n+1En(t) and / = /(k1, . . . , knja) and note that En(t)⁄Ek(t) = En+k(t) where En+k(t) = En+k(t ; k1, . . . , kn, a1, . . . , ak). Part (d) of Theorem 5 implies jEn+k(t) h1Ek(t)j 6 ce(va)th1Ek(t) and jAn(t)⁄Ek(t) / Ek(t)j 6 ce(va)t/Ek(t). So the relative error of the approximations En+k(t) h1Ek(t) and An(t)⁄Ek(t) /Ek(t) is no more than ce(va)t, which implies these approximations are good for large t. Similar interpretations can be made for (19) in Theorem 8 below. Theorem 8 below is the basis for the proof of the error bound (14) and also for the error bound (21) in Theorem 9 below. Theorem 7 establishes an analogous error bound when (16) is generalized by replacing Ek(t) by a more general class of functions. These theorems use the elasticity of a function g(t) defined by E(g(t)) = tg0 (t)/g(t) where prime denotes differentiation. We start out with the following result connecting convolution and elasticity. Theorem 6. Let Ek(t) = Ek(t ; a1, . . . , ak) where a1, . . . , ak are non-negative and let g(t) and h(t) be continuously differentiable functions satisfying g(t) > 0 and h(t) > 0. For part (a) assume E(g(t)) 6 r and E(h(t)) 6 w for t > 0 for some numbers r and w, while for part (b) assume E(g(t)) P r and E(h(t)) P w for t > 0. Then the following are true for t > 0. (a) E(g(t) ⁄ h(t)) 6 r + w + 1 (b) E(g(t) ⁄ h(t)) P r + w + 1 (c) E(Ek(t)) 6 k 1. Proof. To prove (a) let u(t) = g(t)⁄h(t). It is not hard to see that uðtÞ ¼ t
tu0 ðtÞ ¼ t
Z
1
gðtsÞhðtð1 sÞÞds þ t
0
Z
1
tsg 0 ðtsÞhðtð1 sÞÞds þ t
0
Z
1
R1 0
gðtsÞhðtð1 sÞÞ so 0
gðtsÞtð1 sÞh ðtð1 sÞÞds:
0
Using the hypotheses one has
tu0 ðtÞ 6 uðtÞ þ rt
Z 0
1
gðtsÞhðtð1 sÞÞds þ wt
Z
1
gðtsÞhðtð1 sÞÞds ¼ ðr þ w þ 1ÞuðtÞ: 0
which proves (a). The proof of (b) is the same, but with the inequalities reversed. (c) follows from (a) since E(eat) 6 0 if a P 0. h The following theorem establishes an error bound for generalizations of (16) when Ek(t) is replaced by a more general class of functions. Theorem 7. Assume k1, . . . , kn are real numbers, En(t) = En(t ; k1, . . . , kn) and f(t) satisfies f(0) = 0, f(t) > 0 for t > 0, and there are at numbers a, b and r such that a < min{k1, . . . , kn}, a + b < min{k1, . . . , kn}, E(eatf(t)) 6 r for t > 0, and both eatf(t) and ebt dðe dtf ðtÞÞ are Pn non-decreasing for t P 0. Let h = (k1 a) (kn a) and W ¼ r/ðk1 a; . . . ; kn a j bÞ j¼1 1=ðkj aÞ. Then for t > 0 one has
3316
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
f ðtÞð1 W=tÞ 6 hEn ðtÞ f ðtÞ 6 f ðtÞ
ð17Þ
Proof. The right inequality in (17) follows from Theorem 5 (a), so it remains to prove the left inequality. We first consider the case a = 0, In this case the assumptions become (i) f(0) = 0, f(t) > 0 for t > 0 and f(t) is non-decreasing for t P 0, (ii) E(f(t)) 6 r for t > 0, (iii) k1, . . . , kn are all strictly positive, (iv) b < min{k1, . . . , kn} and (v) ebtf0 (t) is non-decreasing for t P 0. In this case we need to show
f ðtÞð1 r/ðk1 ; . . . ; kn jbÞU=tÞ 6 An ðtÞ f ðtÞ
ð18Þ
Pn
where An(t) = k1,n+1En(t) and U ¼ j¼1 1=kj . Let Sn(t) = Sn(t; k1, . . . , kn). Integrating by parts and using Theorem 4 (f), f(0) = 0 and Sn(0) = 1 one obtains An(t)⁄f(t) = f(t) Sn(t)⁄f0 (t). Theorem 5 (a) gives Aj(t ; k1, . . . , kj)⁄f0 (t) 6 /(k1, . . . , kjjb)f 0 (t). Dividing by kj, summing on j and using the fact that /(k1, . . . , kjjb) 6 /(k1, . . . , knjb) one obtains Sn(t)⁄f0 (t) 6 /(k1, . . . , knjb)Uf 0 (t). Using the fact that E(f(t)) 6 r one obtains Sn(t)⁄f0 (t) 6 r/(k1, . . . , knjb)Uf(t)/t. Combining with the previous gives (18). To prove the left inequality in (17) when a is not necessarily zero we use Theorem 4 (a) and (b) to obtain hEn(t)⁄f(t) = hat [En(t ; k1 a, . . . ,km a)⁄(eatf(t)]. Applying (18) to En(t ; k1 a, . . . , kn a)⁄(eatf(t)) one has n(t)⁄f(t) = he
eat f ðtÞð1 W=tÞ 6 hEn ðt; k1 a; . . . ; kn aÞ ðeat f ðtÞÞ Multiplying by eat and using the above formula for hEn(t)⁄f(t) gives the left inequality in (17). h The following theorem is the basis for the proofs of the error bounds (14) and (21). It can also be used for other similar approximations; Example 3 in Section 3 has a typical instance. Theorem 8. Assume k1, . . . , kn and a1, . . . , ak are real numbers and let En(t) = En(t ; k1, . . . , kn) and Ek(t) = Ek(t ; a1, . . . , ak). Assume k P 2 and let a and m be the smallest and second smallest of a1, . . . , ak. Assume m < min{k1, . . . , kn} and let h = (k1 a) (kn a) P and W ¼ ðk 1Þvða; m; k1 ; . . . ; kn Þ nj¼1 1=ðkj aÞ. Then for t > 0 one has
Ek ðtÞð1 W=tÞ 6 hEn ðtÞ Ek ðtÞ 6 Ek ðtÞ
ð19Þ
Proof. We may assume a = a1 and m = a2 since Ek(t) is a symmetric function of the a’s. Apply Theorem 7 to f(t) = Ek(t) with a = a1 and b = a2 a1. By Theorem 4 (b) and Theorem 6 (b) one has E(eatEk(t)) = E(Ek(t;0,a2 a, . . . , ak a)) 6 k 1 for t > 0. By Theorem 4 (b) and (d) one has ebt dðeat Ek ðtÞÞ=dt ¼ eða2 aÞt Ek1 ðt; a2 a; . . . ; ak aÞ. By Theorem 4 (e) one sees that eatEk(t) and ebtd(eatEk(t))/dt are non-decreasing. Therefore, all the hypotheses of Theorem 7 are satisfied and (19) follows from Theorem 7. h Here is the Proof of Theorem 3 which uses Theorems 5 and 8. Proof of Theorem 3. Let Em(t) = Em(t ; l1, . . . , lm), Enm(t) = Enm(t ; lm+1, . . . , ln) and En(t) = En(t ; k1, . . . , kn). Since l1, . . . , ln is a rearrangement of k1, . . . , kn one has En(t) = En(t; l1, . . . , ln) = Enm(t)⁄Em(t) and
k1;n ¼ k1;nþ1 =kn ¼ l1;nþ1 =kn ¼ l1;mþ1 /ðlmþ1 ; lmþ2 ; . . . ; ln jkp Þh=kn where h = (lm+1 kp) (ln kp). Theorem 5 (d) implies jhEnm ðtÞ Em ðtÞ Em ðtÞj 6 c2 eðkq kp Þt Em ðtÞ where c2 = v(l1, . . . , ln) and Theorem 8 implies jhEnm(t)⁄Em(t) Em(t)j 6 Em(t)W/t. Since l1, . . . , ln is a rearrangement of k1, . . . , kn one has v(l1, . . . , ln) = v(k1, . . . , kn). So c2 = c. We multiply the prceeding inequalities by N1(0)b1,nl1,m+1/(lm+1, lm+2, . . . , lnjkp)/kn and use the fact that (N1(0)b1,nl1,m+1/(lm+1,lm+2, . . . , lnjkp)/kn)Em(t) = aMm(t) and (N1(0)b1,nl1,m+1/(lm+1,lm+2, . . . , lnjkp)/kn) hEnm(t)⁄Em(t) = Nn(t) to get (13) and (14). h The following theorem has two error bounds for the approximation (16). (20) is similar to (13) while (21) is similar to (14). As with (13) and (14), the second error bound (21) shows that the relative error of the approximation (16) is small for a larger range of values of t than does the first, but is not applicable in some cases where (20) is. Theorem 9. Assume k1, . . . , kn are non-negative and a1, . . . , ak are arbitrary real numbers with a = min{a1, . . . , ak} and a < min{k1, . . . , kn}. Let Ek(t) = Ek(t; a1, . . . , ak) and Nn(t) = b1,nk1,nEn(t;k1, . . . , kn)⁄Ek(t) and a = b1,n/(k1, . . . , knja)/kn. For (20) let m be the second smallest of k1, . . . , kn, a1, . . . , ak, c = v(k1, . . . , kn, a1, . . . , ak) and assume a < m and m is strictly smaller than the third smallest of k1, . . . , kn, a1, . . . , ak. For (21) assume k P 2 and let m be the second smallest of a1, . . . , ak. Also assume P m < min{k1, . . . , kn} and let W ¼ ðk 1Þvða; m; k1 ; . . . ; kn Þ nj¼1 1=ðkj aÞ. Then for t P 0 one has
jNn ðtÞ aEk ðtÞj 6 ceðmaÞt aEk ðtÞ
ð20Þ
jNn ðtÞ aEk ðtÞj 6 aEk ðtÞW=t
ð21Þ
Proof. Let En(t) = En(t ; k1, . . . , kn). Theorem 5 (d) implies jhEn(t)⁄Ek(t) Ek(t)j 6 ce(ma)t Ek(t) where h = (k1 a) (kn a) and Theorem 8 implies jhEn(t)⁄Ek(t)Ek(t)j 6 Ek(t)W/t. Multiplying by a and using the fact that Nn(t) = haEn(t)⁄Ek(t) proves (20) and (21). h
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
3317
2.3. More than one initial nuclide In this section we consider the case where more than one Nj(0) is non-zero in which case (2) becomes
Nn ðtÞ ¼
n X
Nj;n ðtÞ;
ð22Þ
j¼1
Nj;n ðtÞ ¼ Nj ð0Þbj;n kj;n Enþ1j ðt; kj ; . . . ; kn Þ: We approximate each Nj,n(t) using the reduced chain (9). In addition to (9) we also assume k1, . . . , kn are all distinct and Xr(1), . . ., Xr(m) are the nuclides with the m largest half-lives, i.e. kr(i) < kr(k) for i 6 m and k > m. In fact, one can often choose m so kr(i) << kr(k) for i 6 m and k > m which makes applying the error bound (25) below simpler; see Example 4 in the next section. Using (11) a single term Nj,n(t) in (22) with j 6 r(m) can be approximated for large t by
Nj;n ðtÞ n/j Mj;m ðtÞ M j;m ðtÞ ¼ Nj ð0Þbj;rðmÞ lsðjÞ;m Emþ1sðjÞ ðt; lsðjÞ ; . . . ; lm Þ
ð23Þ
/j ¼ /ðlxðjÞ ; lxðjÞþ1 ; . . . ; ln jkpðjÞ Þ if j 6 rðnÞ and /j ¼ 1 if j > rðnÞ and n = br(m),nlm/kn. Here s(j) is such that Xr(s(j)) ? Xr(s(j)+1) ? ? Xr(m)? is the portion of the reduced chain that only includes nuclides from Xj, . . . , Xn, i.e. s(j) is the smallest i such that j 6 r(i). (s(j) is not defined if j > r(m).) Also x(j) is such that Xr(x(j)), Xr(x(j)+1), . . . , Xr(n) are the nuclides not in the reduced chain that are included in Xj, . . . , Xn, i.e. x(j) is the smallest i such that i > m and j 6 r(i). (x(j) is not defined if j > r(n).) p(j) is such that kp(j) is the smallest of kj, . . . , kn. If j >r(m) then Nj,n(t) 0 for t on the time scale of the reduced chain. Summing (23) one has
Nn ðtÞ nM m ðtÞ M m ðtÞ ¼
rX ðmÞ
ð24Þ
/j M j;m ðtÞ
j¼1
for large t. Note that Mm(t) is the amount of Xr(m) in the reduced chain at time t if the initial conditions are PrðiÞ Mi ð0Þ ¼ j¼rði1Þþ1 /j N j ð0Þbj;rðiÞ for i = 1, . . . , m where we set r(0) = 0. The following theorem gives an error bound for the approximation (24). Theorem 10. Assume the above and N1(0) > 0. For j < n let q(j) be such that kq(j) is the second smallest of kj, . . . , kn. If j 6 r(m 1) P W j ¼ ðm sðjÞÞvðkpðjÞ ; kqðjÞ ; lxðjÞ ; lxðjÞþ1 ; . . . ln Þ nk¼xðjÞ ½lk kpðjÞ 1 with Wj = 0 if j > r(n). Let W = max{Wj:1 6 j
let
6 r(m 1)} if m P 2 and W = 0 if m = 1. For r(m 1) < j 6 r(m) let cj = v(kj, . . . , kn) if j 6 n 1 with cn = 0 if r(m) = n. Let c = max{cj: r(m 1) < j 6 r(m)} and d = kq(r(m1)+1). Let g1 ðtÞ ¼ maxfW=t; ceðdlm Þt g. Let c = kp(r(m)+1) and P Q PrðmÞ l ¼ maxfl1 ; . . . ; lm g. Let V ¼ nj¼rðmÞþ1 N j ð0Þbj;n kpðjÞ wj where wj ¼ n i¼j ki =ðki kpðjÞ Þ and let RðtÞ ¼ j¼1 i–pðjÞ
/j N j ð0Þbj;n lsðjÞ;mþ1 tmsðjÞ =ðm sðjÞÞ! Let g2 ðtÞ ¼ VeðclÞt =RðtÞ if r(m) < n and g2(t) = 0 if r(m) = n. Then for t > 0 one has
jNn ðtÞ nM m ðtÞj 6 ½g1 ðtÞ þ g2 ðtÞnMm ðtÞ:
ð25Þ
Proof. Using the triangle inequality one has jNn(t) nMm(t)j 6 P1 + P2 where
P1 ¼
rX ðmÞ
jN j;n ðtÞ n/j M j;m ðtÞj
j¼1
and
P2 ¼
n X
N j;n ðtÞ:
j¼rðmÞþ1
By (14) one has
jNj;n ðtÞ n/j Mj;m ðtÞj 6 W j n/j M j;m ðtÞ=t 6 g1 ðtÞn/j Mj;m ðtÞ for 1 6 j 6 r(m 1). By (13) one has
jNj;n ðtÞ n/j Mj;m ðtÞj 6 cj eðkqðjÞ kpðjÞ Þt n/j Mj;m ðtÞ 6 g1 ðtÞn/j Mj;m ðtÞ for r(m 1) < j 6 r(m) since p(j) = r(m), kp(j) = kr(m) = lm and kq(j) P d. Summing these two previous inequalities it follows that P1 6 g1(t)nMm(t).
3318
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
For j > r(m) one can write kj;n Enþ1j ðt; kj ; . . . ; kn Þ ¼ ðkpðjÞ =kn Þwj hEnj ðtÞ ekpðjÞ t where Enj(t) = Enj(t ; kj, . . . , kp(j)1, Q kp(j)+1, . . . , kn) and h ¼ n i¼j ðki kpðjÞ Þ. By Theorem 5 (b) one has hEnj ðtÞ ekpðjÞ t 6 ekpðjÞ t 6 ect , so n X
P2 ¼
j¼rðmÞþ1
i–pðjÞ
Nj;n ðtÞ 6
n X
Nj ð0Þbj;n kpðjÞ wj ect =kn ¼ Vect =kn :
j¼rðmÞþ1
Theorem 4 (g) implies Emiþ1 ðt; li ; . . . ; lm Þ P t mi el t =ðm iÞ! for i = 1, . . . , m. Combining with (23) it follows that n/j M j;m ðtÞ P /j N j ð0Þbj;n lsðjÞ;mþ1 tmsðjÞ el t =½kn ðm sðjÞÞ! for j = 1, . . . , r(m). Therefore nM m ðtÞ P RðtÞel t =kn . So P2 6 Vect nM m ðtÞel t =RðtÞ ¼ g2 ðtÞnM m ðtÞ. Combining this with P1 6 g1(t)nMm(t) gives (25). h Remark. g1(t) comes from the error bounds (13) and (14) while g2(t) comes from the bounds in Theorem 4 (g) and Theorem 9(b) This results in g1(t) and g2(t) being a bit complicated. Both g1(t) and g2(t) approach zero as t goes to infinity, so (25) shows that the relative error in the approximation (24) is small for large t.
3. Examples We present four specific examples to illustrate the error bounds in Theorems 1, 2, 3, 8 and 10. This first example is an application of Theorems 1 and 2. Example 1. The chain 223Fr ? 223Ra ? 219Rn ? 215Po? has half-lives T1 = 21.8 min, T2 = 11.4 days, T3 = 3.96 s and T4 = 1.78 103 s and branching fractions b1 = 0.99994 and b2 = b3 = 1; see Firestone et al. [28]. In this case p = 2 and q = 1 in Theorems 1 and 2. Assume that the initial amounts of 223Ra, 219Rn and 215Po are zero. If one approximates the amount N4(t) of 215Po at time t using (4) and (6) with k = 2 then one obtains N 4 ðtÞ N 1 ð0Þb1;4 k1;4 C 2 ek2 t ¼ N 1 ð0Þb1 k1 k2 k3 ððk1 k2 Þðk3 k2 Þðk4 k2 ÞÞ1 ek2 t and N4(t) k2k3((k3 k2)(k4 k2))1N2(t ; k1, k2). If a 2% error is acceptable and one applies (5) and (7), then these approximations are good for t > t0 where t0 = ln(50c)/(k1 k2) where c = (k3 k2)(k4 k2)/(k3 k1)(k4 k1). However, c 1 and k1 k2 k1, so t0 ln(50)/k1 < 2.07 hours. On the other hand, if we apply (8), then N4(t) k2k3((k3 k2)(k4 k2))1N2(t ; k1, k2) with error less than 2% if t > t1 where t1 = 50 c[1/(k3 k2) + 1/ (k4 k2)] 50/k3 < 4.8 min. So the approximation N4(t) k2k3((k3 k2)(k4 k2))1N2(t ; k1, k2) has error less than 2% for t > 4.8 min which is better than the range given by (7). One final note, since b1 1 and kj/(kj k2) 1 for j = 1, 3, 4, the above approximations can be simplified to N 4 ðtÞ N 1 ð0Þðk2 =k4 Þek2 t and N4(t) (k2/k4)N2(t ; k1, k2) which have error less than 2% for t > 2.07 hours and t > 4.8 min respectively. The following example is a typical application of Theorem 3. Example 2. The chain 238U ? 234Th ? 234 mPa ? 234U? has half-lives T1 = 4.47 109 years, T2 = 24.1 days, T3 = 1.17 min and T4 = 2.46 105 years and branching fractions b1 = b2 = 1 and b3 = 0.9984; see Firestone et al. [28]. Here 234 m Pa denotes the metastable isomer of 234Pa. In this case p = 1. Assume that the initial amounts of 234Th, 234 m Pa and 234U are zero. Suppose one approximates the amount N4(t) of 234U at time t using (11) with the reduced chain 238U ? 234U ?. In the notation of (9) one has m = 2, r(1) = 1, r(2) = 4 and the deleted nuclides are 234Th and 234Pa with r(3) = 2 and r(4) = 3. Also, M1(0) = N1(0) and M2(t) = N1(0)b3k1E2(t ; k1, k4) and the approximation (11) is N4(t) /(k2, k3jk1)N1(0)b3k1E2(t ; k1, k4). If a 2% error is acceptable and one applies (14), then this approximation is good for t > t1 where t1 = 50v(k1, . . . , k4)[1/(k2 k1) + 1/ (k3 k1)] 50/k2 < 5 years. Also, /(k2,k3jk1)b3 1, so N4(t) N1(0)k1E2(t ; k1,k4) with error less than 2% for t > 5 years. This approximation is useful, not just for understanding the behavior of N4(t) for t > 5 years, but also because the calculation of N4(t) using (2) may encounter problems with overflow for large t. The following example is a typical application of Theorem 8. Example 3. Consider the series 92Kr ? 92Rb ? 92Sr ? 92 Y? with half-lives T1 = 3.0 s, T2 = 5.3 s, T3 = 2.71 h, T4 = 3.53 h and all bj = 1. Benedict et al. [6, p. 39] considers the case where 92Kr is produced at a constant rate P and the initial amounts of 92 Kr, 92Rb, 92Sr and 92Y are all zero. Let N4(t) be the amount of 92Y at time t. One has N4(t) = k1,4 En(t ; k1, . . . , k4)⁄P and it follows from (15) that N4(t) = Pk1k2k3E2(t ; k1,k2)⁄E3(t ; k3,k4,0). By Theorem 8, N 4 ðtÞ Pk3 E3 ðt; k3 ; k4 ; 0Þ ¼ Pðk4 Þ1 ½1 þ ðk4 k3 Þ1 ðk3 ek4 t k4 ek3 t Þ, an approximation that appears in Benedict et al. [6, p. 40] without an error bound. Note that in Theorem 8, a = 0, m = k4, h = k1k2, v(a,m,k1,k2) 1 and W = 2v(a,m,k1,k2)[1/(k1 a) + 1/(k2 a)] 2(T1 + T2)/(ln2). If an error of 2% error is acceptable then the approximation above is valid with error less than 2% for t > 50W 20 min. The following example is a typical application of Theorem 10. Example 4. Consider the chain 223Fr ? 223Ra ? 219Rn ? 215Po ? in Example 1. Suppose the initial amount of 223Fr, 223Ra, Rn and 215Po are N1(0) = 103 gm, N2(0) = 1.0 gm, N3(0) = 106 gm and N4(0) = 109 gm respectively. Suppose one approximates the amount N4(t) of 215Po at time t using (24) using the reduced chain 223Fr ? 223 Ra ?. In the notation of (9) one
219
F. Massey, J. Prentis / Applied Mathematics and Computation 218 (2011) 3311–3319
3319
has m = 2, r(1) = 1, r(2) = 2 and the deleted nuclides are 219Rn and 215Po with r(3) = 3 and r(4) = 4. Also, M1(0) = /1N1(0), M2(0) = /2N2(0) where /1 = /(k3, k4jk1) and /2 = /(k3, k4jk2). So M 2 ðtÞ ¼ /1 N 1 ð0Þb1 k1 E2 ðt; k1 ; k2 Þ þ /2 N 2 ð0Þek2 t . One has n = b2,4k2/k4 = k2/k4 and the approximation (24) is N 4 ðtÞ ðk2 =k4 Þ½/1 N 1 ð0Þb1 k1 E2 ðt; k1 ; k2 Þ þ /N 2 ð0Þek2 t . Suppose an error of 2% error is acceptable and one applies Theorem 10. One has W = W1 = v(k1, . . . , k4)[1/(k3 k2) + 1/(k4 k2)] 1/k3. Also c = c2 = v(k2, . . . , k4) 1, d = k3 and d lm d. So ceðdlm Þt ek3 t and g1 ðtÞ maxf1=ðk3 tÞ; ek3 t g ¼ 1=ðk3 tÞ. The values of j in ¼ k1 ; c ¼ k3 and c l k3 . Also V = N3(0)b3kp(3)k4(k4 k3)1 + the range r(m) < j 6 n are j = 3 and 4. One has l kp(4)N4(0) = N3(0)k3k4(k4 k3)1 + k4N4(0) k3N3(0) + k4N4(0) and R(t) > /2N2(0)b2,4l2,3 k2N2(0). So g2 ðtÞ K ½k3 N 3 ð0Þþ k4 N 4 ð0Þek3 t =½k2 N 2 ð0Þ 0:8ek3 t . So g1 ðtÞ þ g2 ðtÞ 6 1=ðk3 tÞ þ 0:8ek3 t . So the error in the approximation above for N4(t) will be less than 2% if 1=ðk3 tÞ þ 0:8ek3 t < 0:02. Solving numerically, one sees that this occurs it t > 50/k3 4.8 min. Furthermore, b1 1, /1 1 and /2 1 so the above approximation can be simplified to N 4 ðtÞ ðk2 =k4 Þ½N 1 ð0Þk1 E2 ðt; k1 ; k2 Þ þ N 2 ð0Þek2 t . 4. Conclusions In this paper we have established error bounds for several frequently-used, large-time approximations for a radioactive decay chain X1? X2 ? ? Xn?. These approximations allow one to approximate the amount Nn(t) of Xn in the original chain by the amount of some nucleus in a simpler chain for sufficiently large times. The error bounds in this paper allow users of these approximations to determine how large time must be in order that the approximation hold with a given degree of precision. References [1] E. Segrè, Nuclei and Particles, second ed., Benjamin/Cummings, Reading, Massachusetts, 1977. [2] H. Bateman, Solution of a system of differential equations occurring in the theory of radioactive transformation, Proc. Cambridge. Philos. Soc. 15 (1910) 423–427. [3] A. Mathai, Storage capacity of dam with gamma type inputs, Ann. Inst. Statist. Math. 34 (Part A) (1982) 591–597. [4] J. Cetnar, General solution of Bateman equations for nuclear transmutations, Ann. Nucl. Energy 33 (2006) 640–645. [5] E. Rutherford, J. Chadwick, C.D. Ellis, Radiations from Radioactive Substances, Cambridge University Press, Cambridge, 1951. [6] M. Benedict, T.H. Pigford, H.W. Levi, Nuclear Chemical Engineering, second ed., McGraw-Hill, New York, 1981. [7] S.J. Ball, R.K. Adams, ‘‘MATEXP,’’ A General Purpose Digital Computer Program for Solving Ordinary Differential Equations by the Matrix Exponential Method, Oak Ridge National Laboratory, Oak Ridge, TN, 1967. [8] M.J. Bell, ORIGEN–The ORNL Isotope Generation and Depletion Code, Oak Ridge National Laboratory, Oak Ridge, TN, 1973. [9] G.F. Thomas, D.H. Barber, Numerical solution of Bateman systems, Ann. Nucl. Energy 20 (1993) 407–414. [10] J.M. Hamawi, A useful recurrence formula for the equations of radioactive decay, Nucl. Technol. 11 (1971) 84–88. [11] S.M. Ross, Introduction to Probability Models, eighth ed., Academic Press, Boston, 2003. [12] J. Keilson, U. Sumita, Uniform stochastic ordering and related inequalities, Canad. J. Stat. 10 (1982) 181–198. [13] P.J. Boland, E. El-Neweihi, F. Proschan, Schur properties of convolutions of exponential and geometric random variables, J. Multivariate Anal. 48 (1994) 157–167. [14] J.-L. Bon, E. Pa˘lta˘nea, Ordering properties of convolutions of exponential random variables, Lifetime Data Anal. 5 (1999) 185–192. [15] E. Pa˘lta˘nea, On the approximation of the convolutions of exponential type, in RoGer 2000-Brasßov. Proceedings of the 4thRomanian-German Seminar on Approximation Theory and its Applications (H. Gonska, D. Kacsó and L. Beutel, Eds.), Schriftenreihe des Fachbereichs Mathematik der GerhardMercator-Universität, No. 485, Gerhard-Mercator-Universität, Fachbereich 11/Mathematik, Duisburg, 2000, pp. 107-109. [16] E. Pa˘lta˘nea, On the comparison in hazard rate ordering of fail-safe systems, J. Stat. Planning Infer. 138 (2008) 1993–1997. [17] J. Mi, W. Shi, Y.Y. Zhao, Some properties of convolutions of Pascal and Erlang random variables, Stat. Probab. Lett. 78 (2008) 2378–2387. [18] W. Simonsen, On divided differences and osculatory interpolation, Skandinavisk Aktuarietidsskrift 31 (1948) 157–163. [19] B.M. Minor, Exponential characteristic spatial quadrature for discrete ordinates neutral particle transport in two-dimensional cartesian coordinates, Ph.D. Thesis, Air Force Institute of Technology, Dayton, OH, 1993. [20] K.A. Mathews, G. Sjoden, B. Minor, Exponential characteristic spatial quadrature for discrete ordinates radiation transport in slab geometry, Nuclear Sci. Eng. 118 (1994) 24–37. [21] L.J. Harr, Precise calculation of complex radioactive decay chains, Master’s Thesis, Air Force Institute of Technology, Dayton, OH, 2007. [22] B.V. Erofeev, Zhur. Fiz. Kim. 24 (1950) 721. [23] D.S. Pressyanov, Short solution of the radioactive decay equations, Am. J. Phys 70 (2002) 444–445. [24] L. Moral, A.F. Pacheco, Algebraic approach to the radioactive decay equations, Am. J. Phys 71 (2003) 684–686. [25] W. Kaplan, Operational Methods for Linear Systems, Addison-Wesley, Reading, Massachusets, 1962. [26] B. Hultqvist, Studies on Naturally Occurring Ionizing Radiations, Almqvist & Wiksells Boktrycheri AB, Stockholm, 1956. [27] K. Skrable, C. French, G. Chabot, A. Major, A general equation for the kinetics of linear first order phenomena and suggested applications, Health Phys. 27 (1974) 155–157. [28] R.B. Firestone, V.S. Shirley, C.M. Baglin, S.Y. Chu, J. Zipkin, Table of Isotopes, eighth ed., Wiley, New York, 1996.