指数衰变无法正常工作的问题与解决方案

我们尝试编写一个程序来模拟指数衰变，但它并没有按预期工作。下面是代码示例：

from math import *
from random import random
from time import sleep

class atom():
    def __init__(self, pos):
        self.pos = pos
        self.isalive = True

class substance():
    def __init__(self, halflife):
        self.halflife = float(halflife)
        self.rate = float(log(2) / self.halflife)
        self.life = float(0)
        self.probdecay = float(0)
        self.struct = [atom(i) for i in range(100)]
        self.isalive = True
        print(self.rate)

    def sim(self):
        print(self.probdecay)
        for i in range(len(self.struct)):
            if self.struct[i].isalive:
                if random() < self.probdecay:
                    self.struct[i].isalive = False

        self.life += 1.0
        self.probdecay = 1.0 - pow(1.0 - self.rate, self.life)

    def getstr(self):
        s = ""
        for i in self.struct:
            if i.isalive:
                s += "o"
            else:
                s += " "
        return s

s = substance(5)

while s.isalive:
    s.sim()
    print(s.getstr())
    sleep(1)

运行以上代码，我们得到以下输出：

0.138629436112
0.0
oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
0.138629436112
oooooo o ooooooooooo oooooooo ooooooo o oooooooooooo oooooo oooooooo oooo ooo oooooooo  oooo   oooo 
0.258040751667
o  ooo o   oo oo ooo o    ooo o ooooo o oo  oooooo o   oooo  oo oooo oo   ooo ooo  ooo   ooo   oo   
0.360898143882
   o o     oo     oo o    ooo   ooo o   oo  o o oo     o o    o    o o    ooo o o   oo   o o   oo   
0.449496473814
   o        o        o    o o    o      o     o        o      o      o     o                   oo   
0.525812467226
   o        o                    o      o     o                      o                              
0.591548817506
            o                           o                                                           
0.648172174615

0.696945867656

0.738958091134

0.775146183762

显然，这并不是我们预期的指数衰变。

2、解决方案

问题在于我们混淆了粒子在某一步衰变的概率和粒子到某一步已经衰变的概率。粒子在任何一步衰变的概率都是相同的，而粒子到某一步已经衰变的概率是单调递增的。

为了解决这个问题，我们需要在初始化物质时一次性设置 probdecay，而不是在每次模拟步骤中逐渐增加它。修改后的代码如下：

def __init__(self, halflife):
    self.probdecay = 1 - 0.5 ** (1.0 / halflife)
    ...

现在，我们得到预期的指数衰变输出：

0.138629436112
0.393469340287
oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
0.393469340287
ooooooooooooooooooooooooooooooooooooooooo oooooooooooooooooooooooooooooooooooooooooooooooooo
0.393469340287
ooooooooooooooooooooooooooooo ooooooo ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
0.393469340287
ooooooooooooooooooooooooooooooooooooooooooo oooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
0.393469340287
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo oooooooooooooooooooooooooooooooooooooooooooo
0.393469340287
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo ooooooooooooooooooooooooooooooooooooooooo
0.393469340287
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo oooooooooooooooooooooooooooooooooooooooooo
0.393469340287
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo oooooooooooooooooooooooooooooooooooooooooo
0.393469340287
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo oooooooooooooooooooooooooooooooooooooooooo