一维空间中的吃豆子问题

mathe

我觉得是因为模型是错误的，所以计算结果会和试验结果不同

mathe

23楼方法比如离散成100个点，g,F,W的规模都在5000个点左右。于是积分变成求和（边界可以加上1/2的权值），给定g求F,W变成求线性方程组而给定F,W求g是一个线性规划问题。但是同时作为变量是非线性规划，很难有很好的方法。于是我们可以试验迭代.如果迭代过程不收敛，可以试着使用松弛因子法

mathe

推导了一下，发现对于稳定状态下，最优的时候还是选择g(a,b)=(a+b)/2
我觉得计算机模拟发现(a+b)/2不是最优解的原因是开始的时候，密度分布函数还没有稳定，我们不妨试着把前面部分数据抛弃以后再看看

mathe

我们可以认为使用某个策略g(a,b)经过充分长时间后，分布密度函数为p(a,b,w),其中a,b是两个豆子的坐标0<=a<=b<=1,而w是wayne的坐标。
于是平均长度为
$\int_0^1\int_a^1 C(a,b)dbda$
其中
$C(a,b)=\int_0^a (a-w)p(a,b,w)dw+\int_a^{g(a,b)}(w-a)p(a,b,w)dw+\int_{g(a,b)}^b(b-w)p(a,b,w)dw+\int_b^1(w-b)p(a,b,w)dw$
我们注意到上面表达式中对于给定的a,b,第一项和最后一项是常数，而中间两项取决于g,而对于给定a,b有唯一的g(a,b)使上面表达式最小，而这个表达式关于g的导数是
$(2g-a-b)p(a,b,g)$
所以$g={a+b}/2$时取到最小值,这时二阶导数取值为$2p(a,b,(a+b)/2)>0$，所以是极小值点而不是极大值点

mathe

我现在用Pari/Gp大概计算了一下。但是离散计算时，网格分割的多一点，Pari/Gp就堆栈溢出了。
现在在Windows下只能计算到15×15的网格，所以精度估计有点低。
zgg的表示方法不错，也就是假设给定g(a,b)策略后，稳定下来后，每次吃掉一个豆子时，余下一个豆子c和wayne的位置w可以有一个分布函数p(c,w),而新出来的豆子位置u是满足均匀分布的，可以先不予理会。
而由对称性我们知道p(c,w)=p(1-c,1-w)，所以我们可以光考虑0<w<c<1时的密度函数就可以了。
而我们可以得出关系式
$p(c,w)=\int_0^{g(w,c)} p(c,x)dx +\int_0^{g(w,c)} p(w,x)dx$
即
$p(c,w)=\int_0^{g(w,c)}p(c,x)dx+\int_0^w p(w,x) dx + \int_{1-g(w,c)}^{1-w} p(1-w,x)dx$
离散化p(c,w)后对于给定的g上面就变成一个n(n-1)/2阶线性迭代式，我们需要求最大的特征值对应的特征向量（理论上特征值应该是1），但是我离散化边缘处理的不是很好，特征值略微小于1.
计算特征向量计算量不小，但是上面方程写成迭代方程$p=M*p$后，我们可以通过迭代的方法，任意选择$p_0$,计算$p_n=M^n*p_0$，其中n充分大即可。
由于实际离散化后特征值略微小于1，我们还需要对计算结果$p_n$再次归一化，也就是乘上一个正数使得分量之和为1.

mathe

而得出$p(c,w)$以后，在给定$g(a,b)=(a+b)/2$时计算平均移动位置比较简单。我们同样只考虑$0<w<c<1$的情况（另外一半完成对称，所以结果相同，就不考虑了），也就是我们已经知道$\int_0^1\int_0^c p(c,w)dwdc=1/2$
而对于固定的$c,w,0<w<c<1$,根据u的不同位置，我们可以得出平均移动距离为
在$2w>c$时
$d(c,w)=\int_0^{2w-c}(c-w)du+\int_{2w-c}^w(w-u)du+\int_w^c(u-w)du+\int_c^1(c-w)du=-c^2+2wc-w^2+c-w$
在$2w<c$时为
$d(c,w)=\int_0^w(w-u)du+\int_w^c(u-w)du+\int_c^1(c-w)du=-c^2/2+w^2+c-w$
而最终平均移动距离为
$2\int_0^1\int_0^cp(c,w)d(c,w)dwdc$

mathe

对应Pari/Gp代码，不过由于n最大只能使用15,计算结果精度不高
比如
(09:15) gp > average_weight(10,200)
%34 = 0.2571836267388383841198304306
(09:17) gp > average_weight(11,200)
%35 = 0.2541804007241987873686495312
(09:17) gp > average_weight(12,200)
%36 = 0.2513520585091040215503666987
(09:17) gp > average_weight(13,200)
%37 = 0.2492168669439747355725463563
(09:17) gp > average_weight(14,200)
%38 = 0.2471598401292312838397486659
(09:17) gp > average_weight(15,200)
%39 = 0.2455242992227459358290402694

1. index2to1(a,b)=
   
2. {
   
3.     a*(a-1)/2+b+1
   
4. }
   
5. 
   
6. index2wn(n,a,b)=
   
7. {
   
8.     if(a<b, index2to1(a,b),
   
9.        index2to1(n-a,n-b))
   
10. }
    
11. 
    
12. index1to2(n)=
    
13. {
    
14.     local(a,b);
    
15.     a=floor((sqrt(8*n-7)+1)/2);
    
16.     b=n-a*(a-1)/2-1;
    
17.     [a,b]
    
18. }
    
19. 
    
20. tcount(n)=
    
21. {
    
22.     n*(n+1)/2
    
23. }
    
24. 
    
25. buildm(n,g)=
    
26. {
    
27.     local(t,M,i,b,s);
    
28.     t=tcount(n);
    
29.     M=matrix(t,t);
    
30.     for(u=1,t,
    
31.        i=index1to2(u);
    
32.        b=floor(g[u]);
    
33.        if(b<i[1],
    
34.          for(v=0,b-1,M[u,index2to1(i[1],v)]=1.0/n);
    
35.          M[u,index2to1(i[1],b)]=(g[u]-b)/n,
    
36.          for(v=0,i[1]-1,M[u,index2to1(i[1],v)]=1.0/n);
    
37.          for(v=i[1],b-1,M[u,index2to1(n-i[1],n-v-1)]=1.0/n);
    
38.          M[u,index2to1(n-i[1],n-b-1)]=(g[u]-b)/n
    
39.        );
    
40.        if(b<i[2],
    
41.          for(v=0,b-1,M[u,index2to1(i[2],v)]=1.0/n);
    
42.          M[u,index2to1(i[2],b)]=(g[u]-b)/n,
    
43.          for(v=0,i[2]-1,M[u,index2to1(i[2],v)]=1.0/n);
    
44.          for(v=i[2],b-1,M[u,index2to1(n-i[2],n-v-1)]=1.0/n);
    
45.          M[u,index2to1(n-i[2],n-b-1)]=(g[u]-b)/n
    
46.        );
    
47.     );
    
48.     M
    
49. }
    
50. 
    
51. create_g(n)=
    
52. {
    
53.     local(g,t,i);
    
54.     t=tcount(n);
    
55.     g=vector(t);
    
56.     for(u=1,t,
    
57.        i=index1to2(u);
    
58.        g[u]=(i[1]+i[2])/2.0
    
59.     );
    
60.     g
    
61. }
    
62. 
    
63. build_sm(n)=
    
64. {
    
65.     buildm(n,create_g(n))
    
66. }
    
67. 
    
68. normalize(v)=
    
69. {
    
70.     local(mv,l);
    
71.     l=length(v);mv=0.0;
    
72.     for(u=1,l, mv+=v[u]);
    
73.     for(u=1,l, v[u]/=mv);
    
74.     v
    
75. }
    
76. cost(a,b)=
    
77. {
    
78.    if(a>2*b, a-a^2/2+b^2-b, 2*a*b+a-a^2-b-b^2)
    
79. }
    
80. build_sv(n, i)=
    
81. {
    
82.     local(m,v);
    
83.     m=buildm(n,create_g(n))^i;
    
84.     v=vector(tcount(n),x,1);
    
85.     normalize(v*m~)
    
86. }
    
87. 
    
88. average_weight(n,i)=
    
89. {
    
90.     local(v,t,s,c,w);
    
91.     v=build_sv(n,i);
    
92.     t=tcount(n);s=0.0;
    
93.     for(u=1,t,
    
94.        h=index1to2(u);
    
95.        s+=v[u]*cost(h[1]/(n+0.0),h[2]/(n+0.0));
    
96.     );
    
97.     s
    
98. }

复制代码

mathe

发现可以通过增加内存来扩展点的数目
(09:17) gp > allocatemem(40000000)
(09:19) gp > average_weight(20,200)
%41 = 0.2393806053175264654561361443
(09:21) gp > allocatemem(400000000)
(09:23) gp > average_weight(30,100)
%44 = 0.2331910035260221266557374996
由此可以看出这种近似计算误差还是很大，只能得出一个近似的估计值

mathe

而另外一种更好的方法是通过符号运算得出p的近似表达式，我们可以使用迭代式
$p_{n+1}(c,w)=\int_0^{g(w,c)}p_n(c,x)dx+\int_0^w p_n(w,x) dx + \int_{1-g(w,c)}^{1-w} p_n(1-w,x)dx$
而我们比如可以选择$p_0(c,w)=1$.
而我上面的数值计算发现迭代次数对精度影响好像不是很大，所以比如我们可以算出$p_100(c,w)$,对于$g(a,b)=(a+b)/2$这应该是一个关于c,w的100次多项式，系数在$100^2$项数量级，所以问题应该不大，只是最好直接使用支持符号运算的软件，Pari/Gp就不大适合了

KeyTo9_Fans

我等小辈功力不足，看来还是需要mathe大师亲自出马才有进展。

"由此可以看出这种近似计算误差还是很大，只能得出一个近似的估计值"

虽然如此，但已经可以得到比试验结果的精确程度更高的结果了。

记

$a_10=0.2571836267388383841198304306$
$a_11=0.2541804007241987873686495312$
$a_12=0.2513520585091040215503666987$
$a_13=0.2492168669439747355725463563$
$a_14=0.2471598401292312838397486659$
$a_15=0.2455242992227459358290402694$
$a_16=...$
......
$a_30=...$

通过分析该数列，可以得到精确程度为$10^{-6}$以上的极限值，比试验结果的精确程度高。