找回密码
 欢迎注册
楼主: KeyTo9_Fans

[原创] 合并石子问题(2)

[复制链接]
发表于 2022-8-3 07:28:05 | 显示全部楼层
设$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\)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2022-8-12 21:29:52 | 显示全部楼层
a.png
可以看出,石子数目越大,曲线和指数函数偏离好像越大

点评

是的,曲线在理论值的附近来回震荡  发表于 2022-8-14 09:56
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2022-8-14 10:39:56 | 显示全部楼层
记$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$的极限情况精确地求解出来呢?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2022-8-15 20:34:35 | 显示全部楼层
根据

$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$、……、的通式好像越来越复杂的样子,更不好解了……

点评

石子总数在变少,所以分母要有个(1-dt),才是新的比例,分母里的dt不能忽略  发表于 2022-8-18 09:04
可是按你的公式,应该$p_0(t)=exp(-2t)$才对呀?  发表于 2022-8-16 20:30
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-6-13 05:57 , Processed in 0.043562 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表