我觉得是因为模型是错误的，所以计算结果会和试验结果不同
23楼方法比如离散成100个点，g,F,W的规模都在5000个点左右。于是积分变成求和（边界可以加上1/2的权值），给定g求F,W变成求线性方程组而给定F,W求g是一个线性规划问题。但是同时作为变量是非线性规划，很难有很好的方法。于是我们可以试验迭代.如果迭代过程不收敛，可以试着使用松弛因子法
推导了一下，发现对于稳定状态下，最优的时候还是选择g(a,b)=(a+b)/2
我觉得计算机模拟发现(a+b)/2不是最优解的原因是开始的时候，密度分布函数还没有稳定，我们不妨试着把前面部分数据抛弃以后再看看
我们可以认为使用某个策略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$，所以是极小值点而不是极大值点
我现在用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.
而得出$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$
对应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.2455242992227459358290402694index2to1(a,b)=
{
a*(a-1)/2+b+1
}
index2wn(n,a,b)=
{
if(a<b, index2to1(a,b),
index2to1(n-a,n-b))
}
index1to2(n)=
{
local(a,b);
a=floor((sqrt(8*n-7)+1)/2);
b=n-a*(a-1)/2-1;
}
tcount(n)=
{
n*(n+1)/2
}
buildm(n,g)=
{
local(t,M,i,b,s);
t=tcount(n);
M=matrix(t,t);
for(u=1,t,
i=index1to2(u);
b=floor(g);
if(b<i,
for(v=0,b-1,M,v)]=1.0/n);
M,b)]=(g-b)/n,
for(v=0,i-1,M,v)]=1.0/n);
for(v=i,b-1,M,n-v-1)]=1.0/n);
M,n-b-1)]=(g-b)/n
);
if(b<i,
for(v=0,b-1,M,v)]=1.0/n);
M,b)]=(g-b)/n,
for(v=0,i-1,M,v)]=1.0/n);
for(v=i,b-1,M,n-v-1)]=1.0/n);
M,n-b-1)]=(g-b)/n
);
);
M
}
create_g(n)=
{
local(g,t,i);
t=tcount(n);
g=vector(t);
for(u=1,t,
i=index1to2(u);
g=(i+i)/2.0
);
g
}
build_sm(n)=
{
buildm(n,create_g(n))
}
normalize(v)=
{
local(mv,l);
l=length(v);mv=0.0;
for(u=1,l, mv+=v);
for(u=1,l, v/=mv);
v
}
cost(a,b)=
{
if(a>2*b, a-a^2/2+b^2-b, 2*a*b+a-a^2-b-b^2)
}
build_sv(n, i)=
{
local(m,v);
m=buildm(n,create_g(n))^i;
v=vector(tcount(n),x,1);
normalize(v*m~)
}
average_weight(n,i)=
{
local(v,t,s,c,w);
v=build_sv(n,i);
t=tcount(n);s=0.0;
for(u=1,t,
h=index1to2(u);
s+=v*cost(h/(n+0.0),h/(n+0.0));
);
s
}
发现可以通过增加内存来扩展点的数目
(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
由此可以看出这种近似计算误差还是很大，只能得出一个近似的估计值
而另外一种更好的方法是通过符号运算得出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就不大适合了
我等小辈功力不足，看来还是需要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}$以上的极限值，比试验结果的精确程度高。