\(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\)
可以看出,石子数目越大，曲线和指数函数偏离好像越大 记$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$，结果如下：
p0 0.00000000000009358
p1 0.00000000000029893
p2 0.00000000000091428
p3 0.00000000000278434
p4 0.00000000000847554
p5 0.00000000002579824
p6 0.00000000007852545
p7 0.00000000023901800
p8 0.00000000072752974
p9 0.00000000221447550
10 0.00000000674048285
11 0.00000002051687117
12 0.00000006244982729
13 0.00000019008651625
14 0.00000057859044596
15 0.00000176112762728
16 0.00000536054889374
17 0.00001631639542636
18 0.00004966250436453
19 0.00015114745896772
20 0.00045991248557135
21 0.00139846517444711
22 0.00424347691642835
23 0.01279492803884426
24 0.03784517638466895
25 0.10566980697559133
26 0.24939136999902001
27 0.37382343066513757
28 0.19887733788537268
29 0.01523737438557692
30 0.00003361129914276
31 0.00000000007286201
32 0.00000000000000000
... ...
99 0.00000000000000000
此时相邻两项的比值${p_{10}(t)}/{p_9(t)}\approx 3.043828$，
这就能引出一个猜想：
石子数每增加至原来的$3.043828$倍，合并后的体积变成原来的$2$倍
为了验证这个猜想，继续模拟到$t=30+ln(3.043828)$，结果如下：
p0 0.00000000000003074
p1 0.00000000000009821
p2 0.00000000000030037
p3 0.00000000000091475
p4 0.00000000000278450
p5 0.00000000000847559
p6 0.00000000002579825
p7 0.00000000007852545
p8 0.00000000023901796
p9 0.00000000072752961
10 0.00000000221447510
11 0.00000000674048162
12 0.00000002051686744
13 0.00000006244981592
14 0.00000019008648165
15 0.00000057859034064
16 0.00000176112730671
17 0.00000536054791800
18 0.00001631639245644
19 0.00004966249532533
20 0.00015114743146005
21 0.00045991240189826
22 0.00139846492028206
23 0.00424347614759388
24 0.01279492574262670
25 0.03784516978919188
26 0.10566979019345986
27 0.24939134140238661
28 0.37382343001426172
29 0.19887738401769942
30 0.01523738427381104
31 0.00003361134752197
32 0.00000000007286223
33 0.00000000000000000
... ...
99 0.00000000000000000
这个结果和上面的结果相比，可以发现体积分布并没有发生变化，只是整体平移了一个等级。
因此这个模拟可以证实我的猜想。
但这个模拟比较粗糙，取的是$dt=10^{-7}$
有没有可能把$dt→0$的极限情况精确地求解出来呢？ 根据
$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}$
这个结果和上面模拟的结果：
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}$
这个结果和上面模拟的结果：
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$、……、的通式好像越来越复杂的样子，更不好解了……
页:
1
[2]