mathe 发表于 2008-4-8 18:14
我们同样假设曲线还是类似上面的形式，我们现在先只对$0
检验了一下，好像以前推导的结果中最后一项符号弄错了$-\theta u^3$应该改为$+\theta u^3$
推导方法如下:
$F(\theta,u,u')=\theta \sqrt{u^2+(1+u')^2}$
要求\(\frac{\partial F}{\partial u} =\frac{d}{d\theta}\frac{\partial F}{\partial u'} \)
注意最后一次对$\theta$求导时$u,u'$都要看成$\theta$的函数,而$u'$就是$\frac{du}{d\theta}$
而已知边界条件为$u(2\pi-2d)=\tan(d)$
而求出函数$u(\theta)$后，我们只要画出参数曲线\((\cos(\theta)-u(\theta)\sin(\theta),\sin(\theta)+u(\theta)\cos(\theta))\)
另外根据方程我们容易看出$u'(0)=-1$这个表明在曲线末端曲率半径就是$u(0)$,也就是这个端点到对应圆的切点的距离。这个说明曲线末端是垂直于对应的切线的。
然后我们做变量替换$v=\frac{1+u'}{u}$,得出微分方程$\frac{dv}{d\theta}=(1-v/{\theta})(v^2+1)$, 边界条件$v(0)=0$
如果再做替换$v=ctg(\phi)$,得到$\frac{d\phi}{d\theta}={\ctg(\phi)}/{\theta}-1$,其中几何意义中$\phi$是轨迹上一点的切线和这个点对应的到单位圆（我们的目标圆）的切线的夹角
上面$v$的微分方程和边界条件唯一确定函数$v$.
然后由于$u'-vu+1=0$,假设方程有两个不同的解$u_1,u_2$,得出$(u'_1-u'_2)=v(u_1-v_2)$所以$u_1-u_2=C\exp(\int_0^{\theta}v(x)dx)$
也就是任意算出一个满足条件的函数$u$以后，只要加上$V(\theta)=\exp(\int_0^{\theta}v(x)dx)$的常数倍，就得出$u$的通式，最后我们需要通过方程$u(2\pi-2d)=\tan(d)$计算出对应的d.
而不同的d。每个$u$的不同的通式将有一个不同的d,我们需要从中再挑选出使得总平均距离最小的一条，估计计算量会比较大。但是由于前面已经有个Fans的数值解，我们可以利用那个解来估算这些参数，应该会比较容易找到一个较好的解。
此楼由于使用方程推导错误，结果不对，正确的看57#wayne的结果
可以直接计算v在0的各阶导数，得到$v(x)={x^2}/{2!}-{x^3}/{3!}+{x^4}/{4!}-{x^5}/{5!}+{31x^6}/{6!}-{241x^7}/{7!}+{1221x^8}/{8!}-...$
各阶导数如下
0, 1, -1, 1, -1, 31, -241, 1221, -5057, 33827, -416297, 4802129, -45443873, 404400471, -4427537953,
62816403709, -931431911329, 12963426685387, -178373107004569, 2714248435152393, -47013139883211905,
858451701378457919, -15611830700377251281, 286962627703175278133, -5586716322493471583489,
117035901589119243966579, -2564920873554811193152009, 57222662384675344225171009, -1302531633633982449029492833,
30828952180889290227817222183, -764935521381792405708239398273, 19689970016096416665102700122861,
-519179364888791602076807261741537, 14009456425867369490474633448565019, -389685025651079724342614143622680697,
11217177617603204308620951372491589881, -332824456122265020964776577949071311041,
10119696641293671177424444311102709920783, -314880882090414144821297856390241509480241,
10052042792358929093756335576361193578661157, -329771398244249935956157195605301940824156225,
11100251507060822225427959786616878009196902851, -382309972766801777651620935242662551939590902249,
13458163355732262709457203393834120129909357503537, -484559587633785904893983116919052824766201694037921,
17854599930176026426585448186674826692081445877113591, -672780250004786081496684113315284393474026762599453921,
25889307583182394319969547849331499945604772204423273181
由此我们可以计算出$v(1)$的高精度表示，然后可以在$x=1$再次计算v的各阶导数得到展开式（这时收敛半径应该更大）然后可以得出更大的x的取值，
最后比如可以得出v在$x=\pi$处的展开式，这个展开式应该可以用来计算$[0,2\pi)$上$v$的高精度取值
然后下一步可以得出$V(x)=exp(\int_0^x v(t)dt)$的泰勒展开式
下一步可以挑选一个合适的d,计算对应的u(x),由于我们已经有$u(2\pi-2d)=\tan(d)$,直接根据微分方程$u'(x)-v(x)u(x)+1=0$可以直接根据$v(x)$在$2\pi-2d$的展开式计算出$u(x)$在$2\pi-2d$的展开式，得到一个特解$u_0(x)$。然后一般情况下都有$u(x)=u_0(x)+C*V(x)$,而根据$u(2\pi-2d)=\tan(d)$可以求出其对应的d.
这个pari/gp代码应该可以用来计算满足初始条件$v(h)=a0$时$v$在h的展开式的前n项（除常数项a0)。
比如计算fun(0,0,10)可以得出
%2 = [0, 0.5000000000000000000000000000, -0.1666666666666666666666666667, 0.0416
6666666666666666666666667, -0.008333333333333333333333333333, 0.0430555555555555
5555555555556, -0.04781746031746031746031746032, 0.03028273809523809523809523810
, -0.01393573633156966490299823633, 0.009321814373897707231040564374]
代表在v(0)=0时在x=0展开式为$0*x+0.5*x^2-0.166666..*x^3+0.0146666...*x^4-...$等
fun(a0,h,n)={
local(a,b,b0);
a=vector(n);b=vector(n);
b0=a0*a0;
a=(h-a0)*(1+b0);
b=2.0*a0*a;
a=((h-a0)*b+(1-a)*(1+b0))/2.0;
b=2*a0*a+a^2;
for(s=3,n,
a=(h-a0)*b+1*b;
a-=(1+b0)*a;
for(t=1,s-2,a-=b*a);
a/=s;
b=2*a0*a;
for(t=1,s-1,b+=a*a)
);
a
}
解决了$v$的展开式以后，我们可以用之先得出$u$任意一个特解$u_0$的展开式以及$V$的展开式等。
此后我们知道目标函数为
$A(d)+\int_0^{2\pi-2d} \theta u(\theta)\sqrt{v^2(\theta)+1}d\theta$其中$A(d)=\ln(\frac{1+\sin(d)}{1-\sin(d)})+\frac{2\pi-2d}{\cos(d)}$
$A(d)+\int_0^{2\pi-2d}\theta u_0(\theta)\sqrt{v^2(\theta)+1}d\theta+C\int_0^{2\pi-2d}\theta V(\theta)\sqrt{v^2(\theta)+1}d\theta$
然后记$B(x)=\int_0^x\theta u_0(\theta)\sqrt{v^2(\theta)+1}d\theta,D(x)=\int_0^x\theta V(\theta)\sqrt{v^2(\theta)+1}d\theta$
同样可以推导出B,D的泰勒展开式
而根据$u(\2pi-d)=\tan(d)$可以得出$C=\frac{\tan(d)-u_0(2\pi-2d)}{V(2\pi-2d)}$
代入以后目标函数变成只含一个变量d的函数$"average"(d)=A(d)+B(2\pi-2d)+\frac{\tan(d)-u_0(2\pi-2d)}{V(2\pi-2d)}D(2\pi-2d)$。然后求导为0即可解得d 。
也就是
$\frac{(2\pi-2d)\sin(d)}{\cos^2(d)}-2(2\pi-2d)u_0(2\pi-2d)\sqrt{v^2(2\pi-2d)+1}-2\frac{\tan(d)-u_0(2\pi-2d)}{V(2\pi-2d)}(2\pi-2d) V(2\pi-2d)\sqrt{v^2(2\pi-2d)+1}+\frac{V(2\pi-2d)(\sec^2(d)+2u'_0(2\pi-2d))+2(\tan(d)-u_0(2\pi-2d))V'(2\pi-2d)}{V^2(2\pi-2d)}D(2\pi-2d)=0$
fun(a0,h,n)={ local(a,b,b0); a=vector(n);b=vector(n); b0=a0*a0; a=(h-a0)*(1+b0); b=2.0*a0*a; a=((h-a0)*b+(1-a)*(1+b0))/2.0; b=2*a0*a+a^2; for(s=3,n, a=(h-a0)*b+1*b; a-=(1+b0)*a; for(t=1,s-2,a-=b*a); a/=s; b=2*a0*a; for(t=1,s-1,b+=a*a) ); a }use_fun(a0,h,a,x)={ local(y,s,l); y=x-h;s=0.0;l=length(a); for(u=1,l, s*=y;s+=a; ); s*=y;s+=a0; s }gen_params(n)={ local(a,a0,b,b0,c,c0,d,d0,e,e0,f,f0,g,g0,h,h0,m,m0,k); a=fun(0.0,0.0,n); a0=use_fun(0.0,0.0,a,0.5); b=fun(a0,0.5,n); b0=use_fun(a0,0.5,b,1.0); c=fun(b0,1.0,n); c0=use_fun(b0,1.0,c,1.7); d=fun(c0,1.7,n); d0=use_fun(c0,1.7,d,2.2); e=fun(d0,2.2,n); e0=use_fun(d0,2.2,e,2.6); f=fun(e0,2.6,n); f0=use_fun(e0,2.6,f,3.0); g=fun(f0,3.0,n); g0=use_fun(f0,3.0,g,3.5); h=fun(g0,3.5,n); h0=use_fun(g0,3.5,h,4.2); m=fun(h0,4.2,n); m0=use_fun(h0,4.2,m,5.0); k=fun(m0,5.0,n); }gen_u_params(v_params,d)={ local(n,a0,b,b0,c,c0,dd,d0,e,e0,f,f0,g,g0,h,h0,m,m0,k); n=length(v_params); e=2*Pi-2*d; a0=fun_v(v_params,e);a=fun(a0,e,n); b=gen_u_from_v(tan(d),a0,a);m0=use_fun(tan(d),e,b,5.0); k=gen_u_from_v(m0,v_params,v_params); h0=use_fun(m0,5.0,k,4.2); m=gen_u_from_v(h0,v_params,v_params); g0=use_fun(h0,4.2,m,3.5); h=gen_u_from_v(g0,v_params,v_params); f0=use_fun(g0,3.5,h,3.0); g=gen_u_from_v(f0,v_params,v_params); e0=use_fun(f0,3.0,g,2.6); f=gen_u_from_v(e0,v_params,v_params); d0=use_fun(e0,2.6,f,2.2); e=gen_u_from_v(d0,v_params,v_params); c0=use_fun(d0,2.2,e,1.7); dd=gen_u_from_v(c0,v_params,v_params); b0=use_fun(c0,1.7,dd,1.0); c=gen_u_from_v(b0,v_params,v_params); a0=use_fun(b0,1.0,c,0.5); b=gen_u_from_v(a0,v_params,v_params); } gen_u_from_v(u0,v0,v)={ local(u,n); n=length(v); u=vector(n); u=u0*v0-1; for(h=2,n, u=u0*v+v0*u; for(s=1,h-2, u+=u*v ); u/=h ); u }fun_v(params, x)={ local(r); if(x<0.25, r=use_fun(0.0,0.0,params,x), if(x<0.7, r=use_fun(params,0.5,params,x), if(x<1.4, r=use_fun(params,1.0,params,x), if(x<2.0, r=use_fun(params,1.7,params,x), if(x<2.4, r=use_fun(params,2.2,params,x), if(x<2.8, r=use_fun(params,2.6,params,x), if(x<3.2, r=use_fun(params,3.0,params,x), if(x<3.8, r=use_fun(params,3.5,params,x), if(x<4.5, r=use_fun(params,4.2,params,x), r=use_fun(params,5.0,params,x) ))))))))); r }fun_u(params, x)={ local(r); if(x<0.7, r=use_fun(params,0.5,params,x), if(x<1.4, r=use_fun(params,1.0,params,x), if(x<2.0, r=use_fun(params,1.7,params,x), if(x<2.4, r=use_fun(params,2.2,params,x), if(x<2.8, r=use_fun(params,2.6,params,x), if(x<3.2, r=use_fun(params,3.0,params,x), if(x<3.8, r=use_fun(params,3.5,params,x), if(x<4.5, r=use_fun(params,4.2,params,x), r=use_fun(params,5.0,params,x) )))))))); r }
上面pari/gp代码可分段计算v,u0函数
函数gen_params返回计算v需要的中间数据，然后将其传给fun_v可得出v
然后比如选d=acos(1/1.07),使用函数gen_u_params可以产生计算u0需要的参数，然后传给fun_u可计算对应的u0
此后可以直接使用u0作为u画图得出一条近似曲线，只是我画的纵横比不对，看上去有点怪。
而计算最优d和对应最优值，还得计算V,B,D等函数，工作量将很大
同样使用wayne 39#的代码,将\(-\theta u(t)^3\)改为\(\theta u(t)^3\),可以得到下面图形
自变量范围可限定在0到2pi-2d,d较小比如0.3
此楼数据由于使用v方程不对，所以结果不正确，正确结果查看53#.
现在可以计算出在d=0.590893319813969928977206512531098688032304589325 时取得最小期望a=3.55214056307013612938917610007896899689039684051
附件为对应Pari/Gp代码
在Pari/Gp下比如可以先输入
\p 100
设定计算精度为100位
然后用
\r a.gp
装入附件中gp文件
即可输出最小期望值
（还可以设定更高精度，不过更加花时间）。
代码中将v等函数分成很多段，是因为泰勒展开式的收敛半径不够大，只好这样表示