求不稳定物质的单位时间留存率

KeyTo9_Fans

以前我们讨论过有限状态的物质的衰变：

https://bbs.emath.ac.cn/thread-8767-1-1.html

现在我想研究无限状态的物质的衰变：

某种不稳定物质有无限多种状态，用0、1、2、3、……、n来表示，其中n→∞

其中，状态0是稳定状态，其余状态都不是稳定状态。

对于任意的i>0，处于状态i的物质，在经过1个单位时间后，有10%会变成状态(i+1)，有90%会变成状态(i-1)。

而处于状态0的物质，永远都会留在状态0。

我们假设一开始有100%的物质处于状态1。

而关于单位时间留存率，则需要用极限来定义：

假设经过了$t$（$t\rightarrow\infty$）个单位时间后，处于不稳定状态的物质所占的比例是$p_t$，

那么该物质的单位时间留存率就等于$p_t$的(1/t)次方。

问题1：对于上面所举的例子，该物质的单位时间留存率是多少？

问题2：把题目描述里的10%改成p，90%改成(1-p)，那么当p<0.5时，该物质的单位时间留存率好像是$2\sqrt{p(1-p)}$，当p≥0.5时，该物质的单位时间留存率好像是1。我们该如何证明这个有趣的结论呢？

KeyTo9_Fans

好像大致估算一下，经过了n个单位时间后，处于不稳定状态的物质所占的比例，就能得到第2问的结论。

估算方法如下：

某个粒子能经过n个单位时间还不稳定，最有可能的情况就是进行了(n/2)次(i+1)的转移和(n/2)次(i-1)的转移，

这个概率大概是  卡特兰(n/2) * p^(n/2) * (1-p)^(n/2) ≈ 2^n * (p(1-p))^(n/2)，

这个概率开n次方，等于$2\sqrt{p(1-p)}$，就是该物质的单位时间留存率。

程序模拟第1问的结果，得到的单位时间留存率是0.6，与$2\sqrt{0.1*0.9}$是一致的。

#####

上面这个估算是比较粗糙的，

不知道有没有办法通过矩阵的特征值得到严谨的推导过程呢？

mathe

为了方便起见，我们记$s=\sqrt{p(1-p)}$.
对于有0，1，2，..., n这些状态的情况，对应状态转移矩阵为
\(T_n=\begin{pmatrix}0&1-p&0&\cdots&0&0&0\\
p&0&1-p&\cdots&0&0&0\\
0&p&0&\cdots&0&0&0\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\
0&0&0&\cdots&p&0&1-p\\
0&0&0&\cdots&0&p&0\end{pmatrix}\)
为了求其特征方程，我们可以记$f_n(x)=\det(xE_n-T_n)$是$T_n$的特征多项式（$E_n$为n阶单位阵)
对行列式第一列展开得到$f_n(x)=x*f_{n-1}(x)-s^2 f_{n-2}(x)$而且$f_1(x)=x, f_2(x)=x^2-s^2$, 容易归纳证明n和$f_n$奇偶性相同，所以$f_n$的所有零点关于原点对称。
容易看出如果对于绝对值充分大的x有$|x|>s^2+1$而且$|f_2(x)|>|f_1(x)|$,就可以递推证明$|f_n(x)|$单调增，从而不为0.
也就是只要x满足$|x|>s^2+1$而且$|x^2-s^2|>|x|$,那么$f_n(x)$不等于0，由于可以给出了特征值的一个上界。
另外用数学归纳法可以证明，$f_n(x)$有n个不同的实数根,它的每两个实数根之间正好有一个$f_{n-1}(x)$的实数根，而且在$f_n(x)$的实数根处,$f_{n-1}(x)$的取值依次变号。

另外假设对于$T_n$有某个稳定的概率分布$(x_1,x_2,...,x_n)$,使得
\(T_n \begin{pmatrix}x_1\\x_2\\\vdots\\x_n\end{pmatrix}=\lambda_n \begin{pmatrix}x_1\\x_2\\\vdots\\x_n\end{pmatrix}\)
其中$\x_i \gt 0$而且显然$\lambda_n\ge 0$
我们记$x_0=x_{n+1}=0$,于是得到$p x_{i-1} -\lambda_n x_i+(1-p)x_{i+1}=0$对于一切$i=1,2,...,n$都成立。而且根据递推式，我们显然得出$x_{-1}=-\frac{1-p}{p} x_1\lt 0$
这个递推式特征方程为$px^2-\lambda_n x +(1-p)=0$,如果这个特征方程两个非零根$r_1,r_2$,那么我们有通项公式$x_h = a*r_1^h+b*r_2^h$
根据$x_0=0$得到$a+b=0$, 所以a=-b, 于是有$x_h=a(r_1^h-r_2^h) \gt 0$ 对于一切$h=1,2,...,n$成立，但是$x_{n+1}=0$。假设$r_1,r_2$是两个不同的实数根，显然不成立。
类似可以排除重根情况，所以我们得出$r_1,r_2$是共轭复数根，而且$r_1^{n+1}$是实数。可以得出$r_1$是$\sqrt{\frac{1-p}p} \exp(\frac{t\pi i}{n+1})$,
由此得到$\frac{\lambda_n}p = 2\sqrt{\frac{1-p}p}\cos(\frac{t\pi i}{n+1})$,在t=1时取到最大值$\frac{\lambda_n}p = 2\sqrt{\frac{1-p}p}\cos(\frac{\pi }{n+1})$
也就是最大特征值为$\lambda_n = 2\sqrt{(1-p)p}\cos(\frac{\pi }{n+1})$,在n趋向无穷时逼近$2\sqrt{(1-p)p}$

另外由于$r_1=\sqrt{\frac{1-p}p}  \exp(\frac{t\pi i}{n+1}), r_2=\sqrt{\frac{1-p}p} \exp(\frac{-t\pi i}{n+1})$时, 我们有通项公式$x_u=2a(\frac{1-p}p)^{u/2} \sin(\frac{ut\pi }{n+1}) $, 我们得到了特征值$2\sqrt{(1-p)p}\cos(\frac{t\pi }{n+1})$对应的特征向量的第u个分量取值为$2a(\frac{1-p}p)^{u/2} \sin(\frac{ut\pi }{n+1}) $
可以看出$p\lt \frac{1}{2}$时(我这里的p应该和Fans的定义反了),主要权重会分布在后面,也就是n无穷大时,会一直向后发散.但是$p \gt \frac{1}{2}$时, 总体上越后面越小.
但是有一个问题是,实际上绝对值最大的特征值有两个而不是一个,两者互为相反数.而对应特征向量的区别是正特征值对应的特征向量所有分量大于0,而负特征值对应的特征向量每个分量绝对值和前者相同,但是符号交替变换.

对于初始向量\(\begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix}\), 如果反复使用矩阵$T_n$进行迭代计算,那么我们可以先将\(\begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix}\)表示为各特征向量的组合,也就是需要找出系数$a_1,a_2,...,a_n$
使得\(\begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix}=\sum_{k=1}^n a_k v_k\),其中\(v_t=\begin{pmatrix}(\frac{1-p}p)^{1/2} \sin(\frac{t\pi }{n+1})\\ (\frac{1-p}p)^{2/2} \sin(\frac{2t\pi }{n+1})\\ \cdots\\ (\frac{1-p}p)^{n/2} \sin(\frac{n t\pi }{n+1})\end{pmatrix}\)
也就是\(\begin{pmatrix}y_1(\frac{1-p}p)^{-1/2}\\y_2(\frac{1-p}p)^{-2/2}\\\vdots\\y_n(\frac{1-p}p)^{-n/2}\end{pmatrix}= \begin{pmatrix} \sin(\frac{\pi }{n+1})&\sin(\frac{\pi }{n+1})&\cdots&\sin(\frac{n\pi }{n+1})\\  \sin(\frac{2\pi }{n+1})&\sin(\frac{4\pi }{n+1})&\cdots&\sin(\frac{2n\pi }{n+1})\\ \vdots&\vdots&\vdots&\vdots\\ \sin(\frac{n\pi }{n+1})&\sin(\frac{2n\pi }{n+1})&\cdots&\sin(\frac{n^2\pi }{n+1})\end{pmatrix}\begin{pmatrix}a_1\\a_2\\\vdots\\a_n\end{pmatrix}\)
右边这个正好对应离散正弦变换,除了相差一个常数倍以外,正好是一个对称正交变换阵.所以我们可以解得
\(\begin{pmatrix}a_1\\a_2\\\vdots\\a_n\end{pmatrix}= \frac1{n+1}\begin{pmatrix} \sin(\frac{\pi }{n+1})&\sin(\frac{\pi }{n+1})&\cdots&\sin(\frac{n\pi }{n+1})\\  \sin(\frac{2\pi }{n+1})&\sin(\frac{4\pi }{n+1})&\cdots&\sin(\frac{2n\pi }{n+1})\\ \vdots&\vdots&\vdots&\vdots\\ \sin(\frac{n\pi }{n+1})&\sin(\frac{2n\pi }{n+1})&\cdots&\sin(\frac{n^2\pi }{n+1})\end{pmatrix}\begin{pmatrix}y_1(\frac{1-p}p)^{-1/2}\\y_2(\frac{1-p}p)^{-2/2}\\\vdots\\y_n(\frac{1-p}p)^{-n/2}\end{pmatrix}\)
我们比较关注其中的$a_1,a_n$,对应绝对值最大的两个特征值.而我们实际的初始向量为$y_1=1,y_2=y_3=...=y_n=0$,得到$a_1=a_n=a$
所以可以得到迭代充分大u次以后,向量近似为$a (2s\cos(\frac{\pi i}{n+1}))^u v_1+a (-2s\cos(\frac{\pi i}{n+1}))^u v_n = a (2s\cos(\frac{\pi i}{n+1}))^u (v_1 +(-1)^u v_n)$
于是比较有意思,u为偶数时,向量所有偶数分量抵消为0,只余下奇数分量,u为奇数时,余下向量奇数分量抵消为0,只余下偶数分量.
显然$\frac{1-p}p$是否小于1决定了留存率, 因为\(v_t=\begin{pmatrix}(\frac{1-p}p)^{1/2} \sin(\frac{t\pi }{n+1})\\ (\frac{1-p}p)^{2/2} \sin(\frac{2t\pi }{n+1})\\ \cdots\\ (\frac{1-p}p)^{n/2} \sin(\frac{n t\pi }{n+1})\end{pmatrix}\)