锐客网

  1. 首页 > it技术 > >

CTRW的化学主方程推导(上)

算法

【CTRW的化学主方程推导(上)】连续时间随机行走(continuous time random walk)是一种常见的扩散模型，化学主方程也和扩散有关，所以思考CTRW模型下的化学主方程的推导。公式比较多，分两部分来写。
定义
m s m_{s} ms?：代表不同物种的个数
m r m_{r} mr?：代表不同反应的个数，物种与反应共同组成一个反应系统。
S j S_{j} Sj?：代表不同的物种， j j j的取值从1到 m s m_{s} ms?
n j n_{j} nj?：物种 S j S_{j} Sj?对应的粒子数
n = ( n 1 , ? ? , n m s ) T \mathbf{n}=\left(n_{1}, \cdots, n_{m_{s}}\right)^{\mathrm{T}} n=(n1?,?,nms??)T：粒子数目的状态向量
r i j ∈ N ( p i j ∈ N ) r_{i j} \in \mathbb{N}\left(p_{i j} \in \mathbb{N}\right) rij?∈N(pij?∈N)：经过反应 i i i后 n j n_{j} nj?的变化量。
s j ˉ = p i j ? r i j s_{\bar{j}}=p_{i j}-r_{i j} sjˉ??=pij??rij?：化学计量系数，代表物种 j j j经过反应 i i i后的粒子数目变化
根据以上定义，我们可以把反应 i i i对于状态空间的影响表示为 ∑ j r i j S j → ∑ j p i j S j \sum_{j} r_{i j} S_{j} \rightarrow \sum_{j} p_{i j} S_{j} ∑j?rij?Sj?→∑j?pij?Sj?
内部反应等待时间
单一反应 i i i可以由它的反应等待时间 τ i r \tau_{i}^{r} τir?来刻画（其中上标 r r r是为了与全局的等待时间区分，特指反应的等待时间)， τ i r \tau_{i}^{r} τir?的 PDF，我们记为ψ i r \psi_{i}^{r} ψir?，考虑 ψ i r \psi_{i}^{r} ψir?的表达式：
定义： h i ( n ) h_{i}(\mathbf{n}) hi?(n)为反应 i i i可能发生的反应事件数目，依赖于状态 n n n
那么会有 h i ( n ) = ∏ j n j ! r i j ! ( n j ? r i j ) ! ° h_{i}(\mathbf{n})=\prod_{j} \frac{n_{j} !}{r_{i j} !\left(n_{j}-r_{i j}\right) !} \circ hi?(n)=j∏?rij?!(nj??rij?)!nj?!?°
这 h i h_{i} hi?个可能事件中，任意单一反应事件可以由内部反应等待时间θ l ( i ) ( i i d , l = 1 , ? ? , h i ( n ) ) \theta_{l}^{(i)}\left(i i d, l=1, \cdots, h_{i}(\mathbf{n})\right) θl(i)?(iid,l=1,?,hi?(n))来刻画，其中θ l ( i ) \theta_{l}^{(i)} θl(i)?的 PDF记为 p i p_{i} pi?。我们可以导出：反应 i i i的等待时间τ i r = min ? ( θ l ( i ) ∣ l = 1 , ? ? , h i ( n ) ) \tau_{i}^{r}=\min \left(\theta_{l}^{(i)} \mid l=1, \cdots, h_{i}(\mathbf{n})\right) τir?=min(θl(i)?∣l=1,?,hi?(n))，对应的 PDF 为：
ψ i r ( t ; n ) = h i ( n ) p i ( t ) ∏ l = 1 , l ≠ i h i ( n ) ∫ t ∞ p l ( t l ) d t l = h i ( n ) p i ( t ) [ ∫ t ∞ p i ( t ′ ) d t ′ ] h i ( n ) ? 1 = ? ? ? t [ ∫ t ∞ p i ( t ′ ) d t ′ ] h t ( n ) \begin{aligned} \psi_{i}^{r}(t ; \mathbf{n}) &=h_{i}(\mathbf{n}) p_{i}(t) \prod_{l=1, l \neq i}^{h_{i}(\mathbf{n})} \int_{t}^{\infty} p_{l}\left(t_{l}\right) d t_{l} \\ &=h_{i}(\mathbf{n}) p_{i}(t)\left[\int_{t}^{\infty} p_{i}\left(t^{\prime}\right) d t^{\prime}\right]^{h_{i}(\mathbf{n})-1} \\ &=-\frac{\partial}{\partial t}\left[\int_{t}^{\infty} p_{i}\left(t^{\prime}\right) d t^{\prime}\right]^{h_{t}(\mathbf{n})} \end{aligned} ψir?(t; n)?=hi?(n)pi?(t)l=1,l?=i∏hi?(n)?∫t∞?pl?(tl?)dtl?=hi?(n)pi?(t)[∫t∞?pi?(t′)dt′]hi?(n)?1=??t??[∫t∞?pi?(t′)dt′]ht?(n)?
下一个发生的反应事件取决于哪个反应的等待时间最小，所以内部
反应事件的等待时间满足 τ r = min ? ( τ i r ∣ i = 1 , ? ? , m r ) \tau^{r}=\min \left(\tau_{i}^{r} \mid i=1, \cdots, m_{r}\right) τr=min(τir?∣i=1,?,mr?)，将其对应的PDF记为 ? i r \phi_{i}^{r} ?ir?。 ? i r \phi_{i}^{r} ?ir?同时也是反应等待时间 t t t后发生反应 i i i的联合密度函数。
如果不考虑全局的延迟，那么：
? i r ( t ; n ) = ψ i r ( t ; n ) ∏ l ≠ i ∫ t ∞ ψ l r ( t l ; n ) d t l = h i ( n ) p i ( t ) [ ∫ t ∞ p i ( t ′ ) d t ′ ] h i ( n ) ? 1 ∏ l ? i [ ∫ t ∞ p l ( t ′ ) d t ′ ] t ′ h l ( n ) = h i ( n ) p i ( t ) ∫ t ∞ p i ( t ′ ) d t ′ ∏ l = 1 m r [ ∫ t ∞ p l ( t ′ ) d t ′ ] h l ( n ) \begin{aligned} \phi_{i}^{r}(t ; \mathbf{n}) &=\psi_{i}^{r}(t ; \mathbf{n}) \prod_{l \neq i} \int_{t}^{\infty} \psi_{l}^{r}\left(t_{l} ; \mathbf{n}\right) d t_{l} \\ &=h_{i}(\mathbf{n}) p_{i}(t)\left[\int_{t}^{\infty} p_{i}\left(t^{\prime}\right) d t^{\prime}\right]^{h_{i}(\mathbf{n})-1} \prod_{l * i}\left[\int_{t}^{\infty} p_{l}\left(t^{\prime}\right) d t^{\prime}\right]_{t^{\prime}}^{h_{l}(\mathbf{n})} \\ &=\frac{h_{i}(\mathbf{n}) p_{i}(t)}{\int_{t}^{\infty} p_{i}\left(t^{\prime}\right) d t^{\prime}} \prod_{l=1}^{m_{r}}\left[\int_{t}^{\infty} p_{l}\left(t^{\prime}\right) d t^{\prime}\right]^{h_{l}(\mathbf{n})} \end{aligned} ?ir?(t; n)?=ψir?(t; n)l?=i∏?∫t∞?ψlr?(tl?; n)dtl?=hi?(n)pi?(t)[∫t∞?pi?(t′)dt′]hi?(n)?1l?i∏?[∫t∞?pl?(t′)dt′]t′hl?(n)?=∫t∞?pi?(t′)dt′hi?(n)pi?(t)?l=1∏mr??[∫t∞?pl?(t′)dt′]hl?(n)?
关于内部反应等待时间的推演到这里就可以了。
推广的化学主方程
经历了 k k k步反应后, 我们记CTRW的随机过程的粒子数目向量为 N k \mathbf{N}_{k} Nk?, 时间为 T k T_{k} Tk?，它们有以下递增关系：
N k + 1 = N k + s r k T k + 1 = T k + τ k \mathbf{N}_{k+1}=\mathbf{N}_{k}+\mathbf{s}_{r_{k}} \\ T_{k+1}=T_{k}+\tau_{k} Nk+1?=Nk?+srk??Tk+1?=Tk?+τk?
其中， s r k 1 = ( s r 1 1 , ? ? , s r k m s ) T \mathrm{s}_{r_{k} 1}=\left(s_{r_{1} 1}, \cdots, s_{r_{k} m_{s}}\right)^{\mathrm{T}} srk?1?=(sr1?1?,?,srk?ms??)T， r k ∈ ( 1 , ? ? , m r ) r_{k} \in\left(1, \cdots, m_{r}\right) rk?∈(1,?,mr?)， ( r k , τ k ) \left(r_{k}, \tau_{k}\right) (rk?,τk?)的联合分布由 ? i ( t ; n ) \phi_{i}(t ; \mathbf{n}) ?i?(t; n)给出，初始条件为N 0 = n 0 \mathbf{N}_{0}=\mathbf{n}_{0} N0?=n0?， T 0 = 0 T_{0}=0 T0?=0。
我们可以看到，因为 ( r k , τ k ) \left(r_{k}, \tau_{k}\right) (rk?,τk?)的联合分布依赖于当前的系统状态N k \mathbf{N}_{k} Nk?， 所以这样的递增关系定义了一个非同质的多维CTRW过程。所以我们可以利用CTRW的理论来推导对应的化学主方程。
由于我们更希望用时间 t t t来描述粒子数目的演化过程，所以我们考虑随时间增加步数的更新过程 K t K_{t} Kt?,并记 N ( t ) = N K\mathbf{N}(t)=\mathbf{N}_{K\text { }} N(t)=NK ?，显然过程K t K_{t} Kt?是 T k T_{k} Tk?的共轭，所以 K t = sup ? { k ∣ T k < t } K_{t}=\sup \left\{k \mid T_{k}P ( n , t ) = ? δ n , N K ? P(\mathbf{n}, t)=\left\langle\delta_{\mathbf{n}, \mathbf{N}_{K}}\right\rangle P(n,t)=?δn,NK???
其中 ? ? ? \langle·\rangle ???表示对等待时间过程 τ k \tau_{k} τk?的所有情况求平均， δ i , j \delta_{i, j} δi,j?为Kronecker delta记号。
由此，我们通过单位分解可得： P ( n , t ) = ∑ k = 0 ∞ ? δ n , N k δ k , K t ? P(\mathbf{n}, t)=\sum_{k=0}^{\infty}\left\langle\delta_{\mathbf{n}, \mathbf{N}_{k}} \delta_{k, K_{t}}\right\rangle P(n,t)=∑k=0∞??δn,Nk??δk,Kt???
与此同时，我们考虑示性函数，有 δ k , K t = 1 [ T k , T t + 1 ) ( t ) \delta_{k, K_{t}}=\mathbf{1}_{\left[T_{k}, T_{t+1}\right)}(t) δk,Kt??=1[Tk?,Tt+1?)?(t)，意即：
P ( n , t ) = ∫ 0 ∞ ∑ k = 0 ∞ ? δ n , N k δ ( t ′ ) 1 [ T k , T k + 1 ( t ? t ′ ) ? d t ′ = ∫ 0 ∞ ∑ k = 0 ∞ ? δ n , N k δ ( T k ? t ′ ) 1 [ 0 , τ k ) ( t ? t ′ ) ? d t ′ \begin{aligned} P(\mathbf{n}, t) &=\int_{0}^{\infty} \sum_{k=0}^{\infty}\left\langle\delta_{\mathbf{n}, \mathbf{N}_{k}} \delta\left(t^{\prime}\right) \mathbf{1}_{\left[T_{k}, T_{k+1}\right.}\left(t-t^{\prime}\right)\right\rangle d t^{\prime} \\ &=\int_{0}^{\infty} \sum_{k=0}^{\infty}\left\langle\delta_{\mathbf{n}, \mathbf{N}_{k}} \delta\left(T_{k}-t^{\prime}\right) \mathbf{1}_{\left[0, \tau_{k}\right)}\left(t-t^{\prime}\right)\right\rangle d t^{\prime} \end{aligned} P(n,t)?=∫0∞?k=0∑∞??δn,Nk??δ(t′)1[Tk?,Tk+1??(t?t′)?dt′=∫0∞?k=0∑∞??δn,Nk??δ(Tk??t′)1[0,τk?)?(t?t′)?dt′?
那么，就很容易能得到
P ( n , t ) = ∫ 0 t ∑ k = 0 ∞ R k ( n , t ′ ) ∑ i = 1 m r ∫ t ? t ′ ∞ ? i ( t ′ ′ ; n ) d t ′ ′ d t ′ P(\mathbf{n}, t)=\int_{0}^{t} \sum_{k=0}^{\infty} R_{k}\left(\mathbf{n}, t^{\prime}\right) \sum_{i=1}^{m_{r}} \int_{t-t^{\prime}}^{\infty} \phi_{i}\left(t^{\prime \prime} ; \mathbf{n}\right) d t^{\prime \prime} d t^{\prime} P(n,t)=∫0t?k=0∑∞?Rk?(n,t′)i=1∑mr??∫t?t′∞??i?(t′′; n)dt′′dt′
其中， R k ( n , t ) = ? δ n , N k δ ( T k ? t ) ? R_{k}(\mathbf{n}, t)=\left\langle\delta_{\mathbf{n}, \mathbf{N}_{k}} \delta\left(T_{k}-t\right)\right\rangle Rk?(n,t)=?δn,Nk??δ(Tk??t)?是经历 k k k步反应后在 t t t时刻状态为n \mathbf{n} n的联合密度， R ( n , t ) = ∑ k = 0 ∞ R k ( n , t ) R(\mathbf{n}, t)=\sum_{k=0}^{\infty} R_{k}(\mathbf{n}, t) R(n,t)=∑k=0∞?Rk?(n,t)是经历任意步反应后在 t t t时刻状态为 n \mathbf{n} n的概率密度。
以上的这些推导的作用在于：我们给出了一个直观的物理解释: 系统在更早的时刻t ′ t^{\prime} t′到达状态 n n n ，然后 在接下米的 t ? t ′ t-t^{\prime} t?t′吋间内没有发生反应，再对任意到达时间 t ′ t^{\prime} t′进行积分，这样就得到了在 t t t时刻状态为 n n n的概率。
欲知后事如何，请移步
CTRW的化学主方程推导(下)

    推荐阅读


    上一篇：xls与csv转换（多工作表合并与拆分）

    下一篇：修改ARM|修改ARM linux的系统时钟和硬件时钟