数值计算表明，除了n=4以外，其他都是取切线？？？s0()=
{
sqrt(2+sqrt(5))
}
f0(x)=
{
(sqrt(x^4+x^2-1)+x)/(x^2-1)
}
df0(x)=
{
local(t);
t=sqrt(x^4+x^2-1);
-((x^2+1)*t+3*x^3-x)/((x^4-2*x^2+1)*t)
}
pl(n)=
{
local(s,b,e);
b=s0();
e=2*b;
s=-1.0/(n-1.0);
while(df0(e)<=s,e=e*2);
solve(X=b,e,df0(X)-s)
}
getsr(n)=
{
local(b,ex,ey,e);
b=s0();
ex=pl(n);
ey=f0(ex);
e=(ex+(n-1)*ey)/n;
}
getub(n)=
{
local(b,e);
b=pl(n);
e=2*b;
while(e+(n-1)*f0(e)<n*s0(),e=2*e);
}
fx(n,b,e,s)=
{
solve(X=b,e,X+(n-1)*f0(X)-n*s)
}
hx(n,b,e,s)=
{
local(x,y);
x=fx(n,b,e,s);
y=(n*s-x)/(n-1);
log(x+1/x)+(n-1)*log(y+1/y)
}
pf(n)=
{
local(ub,sr);
ub=getub(n);
sr=getsr(n);
plot(s=sr+0.000001,sr,n*log(s/n+n/s)-hx(n,ub,ub,s));
}
呵呵，上面代码错了一个小地方。（主要由于hujunhua和我使用的记号不同迷惑了我）
改了以后就可以看出解总是存在而且唯一
s0()=
{
sqrt(2+sqrt(5))
}
f0(x)=
{
(sqrt(x^4+x^2-1)+x)/(x^2-1)
}
df0(x)=
{
local(t);
t=sqrt(x^4+x^2-1);
-((x^2+1)*t+3*x^3-x)/((x^4-2*x^2+1)*t)
}
pl(n)=
{
local(s,b,e);
b=s0();
e=2*b;
s=-1.0/(n-1.0);
while(df0(e)<=s,e=e*2);
solve(X=b,e,df0(X)-s)
}
getsr(n)=
{
local(b,ex,ey,e);
b=s0();
ex=pl(n);
ey=f0(ex);
e=(ex+(n-1)*ey)/n;
}
getub(n)=
{
local(b,e);
b=pl(n);
e=2*b;
while(e+(n-1)*f0(e)<n*s0(),e=2*e);
}
fx(n,b,e,s)=
{
solve(X=b,e,X+(n-1)*f0(X)-n*s)
}
hx(n,b,e,s)=
{
local(x,y);
x=fx(n,b,e,s);
y=(n*s-x)/(n-1);
log(x+1/x)+(n-1)*log(y+1/y)
}
pf(n)=
{
local(ub,sr);
ub=getub(n);
sr=getsr(n);
plot(s=sr+0.000001,sr,n*log(s+1/s)-hx(n,ub,ub,s));
}
添加一个函数计算边界s的数值解，采用hujunhua的记号，即s是所有$x_i$的平均值:sf(n)=
{
local(ub,sr);
ub=getub(n);
sr=getsr(n);
solve(s=sr+0.0000000001,sr,n*log(s+1/s)-hx(n,ub,ub,s))
}(06:56) gp > sf(3)
%6 = 2.002867408063500179411985161
(06:56) gp > sf(4)
%7 = 1.928358916474980428360084936
(06:56) gp > sf(5)
%8 = 1.864007650000524850777355005
(06:56) gp > sf(6)
%9 = 1.810514381228163166929447075
(06:56) gp > sf(7)
%10 = 1.765782350018739829894444919
(06:56) gp > sf(8)
%11 = 1.727864312049133923676339702
(06:56) gp > sf(9)
%12 = 1.695269424440409338819703692
(06:56) gp > sf(10)
%13 = 1.666891983135502143928659033
(06:56) gp > sf(20)
%14 = 1.501387153741824255507396749
(06:56) gp > sf(100)
%15 = 1.250349255508031551479807498
(06:56) gp > sf(1000)
%16 = 1.090580657781168832151185931
(06:57) gp > sf(2)
%17 = 2.058171027371492250321981048
将plot函数改成ploth可以画出高精度的图（但是速度很慢）
如图，4个图分别是n=3,4,100,1000的情况
现在我们记$f(x)=ln(x+1/x),x>0$,那么我们现在的简化后的问题是已知
$x+(n-1)y=n*s$,求$f(x)+(n-1)f(y)$的最小值
我们可以改目标函数为$g(y)=f(n*s-(n-1)*y)+(n-1)f(y)$
于是$g'(y)=-(n-1)f'(x)+(n-1)f'(y)$,所以在$f'(y)=f'(x)$时才有$g'(y)=0$
同样$g''(y)=(n-1)^2f''(x)+(n-1)f''(y)$
其中$f'(y)={y^2-1}/{y^3+y},f''(y)={1+4y^2-y^4}/{y^2(y^2+1)^2}$
23# mathe
能保证精度吗？
我是说，有误差分析吗
我们应该改为在约束条件${(x+(n-1)y=ns),(f'(x)=f'(y)):}$之下，求目标函数$f(x)+(n-1)f(y)-nf(s)$的情况，而我们需要证明这个目标函数关于s是单调的，也就是说，上面目标函数只能在边界取到最值。
对这个函数采用拉格朗日极值法，得到
${(f'(x)+lambda_1+lambda_2*f''(x)=0),((n-1)f'(y)+(n-1)lambda_1-lambda_2*f''(y)=0),(-nf'(s)-nlambda_1=0):}$
于是$lambda_1=-f'(s)$,再由$f'(x)=f'(y)$这个已知条件得出$f'(s)=f'(x)$或$(n-1)f''(x)+f''(y)=0$
由于我们已经知道在我们讨论的范围内必然有$f'(s)!=f'(x)$,所以只能取第二种可能。
也就是我们需要证明方程
${(f'(x)=f'(y)),((n-1)f''(x)+f''(y)=0):}$没有符合要求的解
23# mathe
能保证精度吗？
我是说，有误差分析吗
wayne 发表于 2010-7-8 12:48 http://bbs.emath.ac.cn/images/common/back.gif
计算精度应该足够高了，而且计算结果也很合理。现在余下问题主要就在于如果能够证明差值函数是关于s的单调增函数就好办了
23#的结果与14#的结果，小数点前5位是相同的。我没带本本（太重，2.7公斤），无法精算，只能用execel或者几何画板计算，可以保证前5位。
现在终于敢说：wayne在8#算的6.0043029605431965656237025541486不正确。因为6.00430.../3=2.00143...<2.00286...
呵呵，现在找到解决方法了。
我们查看hujunhua定义的曲线$x^2y^2=(x+y)^2+1$,它等价于$f'(x)=f'(y)$
把后面的表达式看成y关于x的隐函数，我们可以得到$dy/dx={f''(x)}/{f''(y)}$
于是27#最后的要求相当于在曲线$x^2y^2=(x+y)^2+1$上取一点，这个点处切线的斜率为$-1/{n-1}$,于是这个唯一的极值点对应与我们一直关注的那条切线的切点，也就是说正好在我们讨论的范围的边界上。而对于我们感兴趣的那个区域，整个区域里面不在存在任何极值点，所以如果将目标函数$f(x)+(n-1)f(y)-nf(s)$看成s的隐函数，在我们讨论的范围里面它是单调的，这个正好解释了24#的图像。也就是说，我上面用二分法来计算总是不会出问题的，而我们讨论的s的分界点是唯一存在的。