合并石子问题(2)

mathe

设$p_k(N)\to f(k-\log_d(N))$,让N趋向无穷，并且$x=k-\log_d(N)$我们可以得到积分方程
\(f(x)=\int_0^1 f(x-1-\log_d(t))f(x-1-\log_d(1-t))dt+2\sum_{s=1}^{\infty}\int_0^1 f(x-\log_d(t))f(x-s-\log_d(1-t))dt\)
要求解满足对于任意x都有\(\sum_{s=-\infty}^{+\infty}f(x+s)=1\)

mathe

可以看出,石子数目越大，曲线和指数函数偏离好像越大

KeyTo9_Fans

记$p_i(t)$为$t$时刻体积为$2^i$的石子个数所占的比例，

由于一开始，体积为$2^0$的石子所占的比例是$1$，其余体积的石子所占的比例是$0$，

所以我们可以将初始状态记为：

$p_0(0)=1$，$p_i(0)=0$（$i=1,2,3,...,99$）

以后假设每过$dt$的时间就合并$dt$比例的石子，

那么任意给定时刻$t$，都可以根据此刻各种体积的石子所占的比例$p_0(t)$、$p_1(t)$、$p_2(t)$、……

推算得到$(t+dt)$时刻各种体积的石子所占的比例$p_0(t)$、$p_1(t)$、$p_2(t)$、……如下：

$p_0(t+dt)=(p_0(t)-2p_0(t)dt)/(1-dt)$

$p_1(t+dt)=(p_1(t)+(p_0^2(t)-2p_1(t)\sum_{i=1}^99 p_i(t))dt)/(1-dt)$

$p_2(t+dt)=(p_2(t)+(p_1^2(t)-2p_2(t)\sum_{i=2}^99 p_i(t))dt)/(1-dt)$

$p_3(t+dt)=(p_3(t)+(p_2^2(t)-2p_3(t)\sum_{i=3}^99 p_i(t))dt)/(1-dt)$

$p_4(t+dt)=(p_4(t)+(p_3^2(t)-2p_4(t)\sum_{i=4}^99 p_i(t))dt)/(1-dt)$

……

$p_{99}(t+dt)=(p_{99}(t)+(p_{98}^2(t)-2p_{99}(t)\sum_{i=99}^{99} p_i(t))dt)/(1-dt)$

我令$dt=10^{-7}$，然后从$t=0$模拟到$t=30$，结果如下：

1. p0 0.00000000000009358
   
2. p1 0.00000000000029893
   
3. p2 0.00000000000091428
   
4. p3 0.00000000000278434
   
5. p4 0.00000000000847554
   
6. p5 0.00000000002579824
   
7. p6 0.00000000007852545
   
8. p7 0.00000000023901800
   
9. p8 0.00000000072752974
   
10. p9 0.00000000221447550
    
11. 10 0.00000000674048285
    
12. 11 0.00000002051687117
    
13. 12 0.00000006244982729
    
14. 13 0.00000019008651625
    
15. 14 0.00000057859044596
    
16. 15 0.00000176112762728
    
17. 16 0.00000536054889374
    
18. 17 0.00001631639542636
    
19. 18 0.00004966250436453
    
20. 19 0.00015114745896772
    
21. 20 0.00045991248557135
    
22. 21 0.00139846517444711
    
23. 22 0.00424347691642835
    
24. 23 0.01279492803884426
    
25. 24 0.03784517638466895
    
26. 25 0.10566980697559133
    
27. 26 0.24939136999902001
    
28. 27 0.37382343066513757
    
29. 28 0.19887733788537268
    
30. 29 0.01523737438557692
    
31. 30 0.00003361129914276
    
32. 31 0.00000000007286201
    
33. 32 0.00000000000000000
    
34. ... ...
    
35. 99 0.00000000000000000

复制代码

此时相邻两项的比值${p_{10}(t)}/{p_9(t)}\approx 3.043828$，

这就能引出一个猜想：

石子数每增加至原来的$3.043828$倍，合并后的体积变成原来的$2$倍

为了验证这个猜想，继续模拟到$t=30+ln(3.043828)$，结果如下：

1. p0 0.00000000000003074
   
2. p1 0.00000000000009821
   
3. p2 0.00000000000030037
   
4. p3 0.00000000000091475
   
5. p4 0.00000000000278450
   
6. p5 0.00000000000847559
   
7. p6 0.00000000002579825
   
8. p7 0.00000000007852545
   
9. p8 0.00000000023901796
   
10. p9 0.00000000072752961
    
11. 10 0.00000000221447510
    
12. 11 0.00000000674048162
    
13. 12 0.00000002051686744
    
14. 13 0.00000006244981592
    
15. 14 0.00000019008648165
    
16. 15 0.00000057859034064
    
17. 16 0.00000176112730671
    
18. 17 0.00000536054791800
    
19. 18 0.00001631639245644
    
20. 19 0.00004966249532533
    
21. 20 0.00015114743146005
    
22. 21 0.00045991240189826
    
23. 22 0.00139846492028206
    
24. 23 0.00424347614759388
    
25. 24 0.01279492574262670
    
26. 25 0.03784516978919188
    
27. 26 0.10566979019345986
    
28. 27 0.24939134140238661
    
29. 28 0.37382343001426172
    
30. 29 0.19887738401769942
    
31. 30 0.01523738427381104
    
32. 31 0.00003361134752197
    
33. 32 0.00000000007286223
    
34. 33 0.00000000000000000
    
35. ... ...
    
36. 99 0.00000000000000000

复制代码

这个结果和上面的结果相比，可以发现体积分布并没有发生变化，只是整体平移了一个等级。

因此这个模拟可以证实我的猜想。

但这个模拟比较粗糙，取的是$dt=10^{-7}$

有没有可能把$dt→0$的极限情况精确地求解出来呢？

KeyTo9_Fans

根据

$p_0(0)=1$，

$p_0(t+dt)=(p_0(t)-2p_0(t)dt)/(1-dt)$

可以解出

$p_0(t)=e^{-t}$

把$t=30$代入，得

$p_0(30)=e^{-30}\approx 9.358*10^{-14}$

这个结果和上面模拟的结果：

1. p0 0.00000000000009358

复制代码

是一致的

至此，我们成功地解出了$p_0$的通式

但$p_1$、$p_2$、……、的通式好像越来越复杂的样子，不是很好解

比如，要解这个式子：

$p_1(t+dt)=(p_1(t)+(p_0^2(t)-2p_1(t)\sum_{i=1}^99 p_i(t))dt)/(1-dt)$

就得先把$p_0(t)=e^{-t}$代入，得到：

$p_1(t+dt)=(p_1(t)+(e^{-2t}-2p_1(t)(1-e^{-t}))dt)/(1-dt)$

然后化简一下，得到：

$p_1(t+dt)=p_1(t)+(e^{-2t}-p_1(t)(1-2e^{-t}))dt$

若记$x=p_1(t)$，则有

$x(0)=0$，

${dx}/{dt}=e^{-2t}-(1-2e^{-t})x$

这个微分方程我就不知道怎么解了……

用Matlab硬解，结果如下：

$x={e^{-t}(e^{2(1-e^{-t})}-1)}/2$

也就是

$p_1(t)={e^{-t}(e^{2(1-e^{-t})}-1)}/2$

把$t=30$代入，得

$p_1(30)={e^{-30}(e^{2(1-e^{-30})}-1)}/2\approx 2.9893*10^{-13}$

这个结果和上面模拟的结果：

1. p1 0.00000000000029893

复制代码

是一致的

至此，我们成功地解出了$p_1$的通式

但$p_2$、$p_3$、……、的通式好像越来越复杂的样子，不是很好解

比如，要解这个式子：

$p_2(t+dt)=(p_2(t)+(p_1^2(t)-2p_2(t)\sum_{i=2}^99 p_i(t))dt)/(1-dt)$

就得先把$p_1(t)={e^{-t}(e^{2(1-e^{-t})}-1)}/2$代入，得到：

$p_2(t+dt)=(p_2(t)+({e^{-2t}(e^{2(1-e^{-t})}-1)^2}/4-2p_2(t)(1-e^{-t}-{e^{-t}(e^{2(1-e^{-t})}-1)}/2))dt)/(1-dt)$

然后化简一下，得到：

$p_2(t+dt)=p_2(t)+({e^{-2t}(e^{2(1-e^{-t})}-1)^2}/4-p_2(t)(1-2e^{-t}-e^{-t}(e^{2(1-e^{-t})}-1)))dt$

若记$x=p_2(t)$，则有

$x(0)=0$，

${dx}/{dt}={e^{-2t}(e^{2(1-e^{-t})}-1)^2}/4-(1-2e^{-t}-e^{-t}(e^{2(1-e^{-t})}-1))x$

这个微分方程我更不知道怎么解了……

尝试用Matlab解，出来一个巨复杂的式子：

$p_2(t)=1/4 (e^{e^{2 - 2e^{-t}}/2 - e^{-t} - t}\int_0^t e^{-x-3e^{-x}-e^{2-2e^{-x}}/2}(e^{2e^{-x}}-e^2)^2 dx)$

经代入数据验证，发现这个式子是对的

但$p_3$、$p_4$、……、的通式好像越来越复杂的样子，更不好解了……