微分方程数值求解

mathe

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的数值解，我们可以利用那个解来估算这些参数，应该会比较容易找到一个较好的解。

mathe

此楼由于使用方程推导错误，结果不对，正确的看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)$的泰勒展开式

mathe

下一步可以挑选一个合适的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.

mathe

这个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-...$等

1. 
   
2. fun(a0,h,n)={
   
3.    local(a,b,b0);
   
4.    a=vector(n);b=vector(n);
   
5.    b0=a0*a0;
   
6.       a[1]=(h-a0)*(1+b0);
   
7.       b[1]=2.0*a0*a[1];
   
8.       a[2]=((h-a0)*b[1]+(1-a[1])*(1+b0))/2.0;
   
9.       b[2]=2*a0*a[2]+a[1]^2;
   
10.       for(s=3,n,
    
11.           a[s]=(h-a0)*b[s-1]+1*b[s-2];
    
12.           a[s]-=(1+b0)*a[s-1];
    
13.           for(t=1,s-2,a[s]-=b[t]*a[s-1-t]);
    
14.           a[s]/=s;
    
15.           b[s]=2*a0*a[s];
    
16.           for(t=1,s-1,b[s]+=a[t]*a[s-t])
    
17.       );
    
18.    a
    
19. }
    

复制代码

mathe

解决了$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$

mathe

1. fun(a0,h,n)={    local(a,b,b0);    a=vector(n);b=vector(n);    b0=a0*a0;       a[1]=(h-a0)*(1+b0);       b[1]=2.0*a0*a[1];       a[2]=((h-a0)*b[1]+(1-a[1])*(1+b0))/2.0;       b[2]=2*a0*a[2]+a[1]^2;       for(s=3,n,           a[s]=(h-a0)*b[s-1]+1*b[s-2];           a[s]-=(1+b0)*a[s-1];           for(t=1,s-2,a[s]-=b[t]*a[s-1-t]);           a[s]/=s;           b[s]=2*a0*a[s];           for(t=1,s-1,b[s]+=a[t]*a[s-t])       );    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[l+1-u];     );     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);    [a,a0,b,b0,c,c0,d,d0,e,e0,f,f0,g,g0,h,h0,m,m0,k] }  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[1]);    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[18],v_params[19]);    h0=use_fun(m0,5.0,k,4.2);    m=gen_u_from_v(h0,v_params[16],v_params[17]);    g0=use_fun(h0,4.2,m,3.5);    h=gen_u_from_v(g0,v_params[14],v_params[15]);    f0=use_fun(g0,3.5,h,3.0);    g=gen_u_from_v(f0,v_params[12],v_params[13]);    e0=use_fun(f0,3.0,g,2.6);    f=gen_u_from_v(e0,v_params[10],v_params[11]);    d0=use_fun(e0,2.6,f,2.2);    e=gen_u_from_v(d0,v_params[8],v_params[9]);    c0=use_fun(d0,2.2,e,1.7);    dd=gen_u_from_v(c0,v_params[6],v_params[7]);    b0=use_fun(c0,1.7,dd,1.0);    c=gen_u_from_v(b0,v_params[4],v_params[5]);    a0=use_fun(b0,1.0,c,0.5);    b=gen_u_from_v(a0,v_params[2],v_params[3]);    [a0,b,b0,c,c0,dd,d0,e,e0,f,f0,g,g0,h,h0,m,m0,k] } gen_u_from_v(u0,v0,v)={    local(u,n);    n=length(v);    u=vector(n);    u[1]=u0*v0-1;    for(h=2,n,       u[h]=u0*v[h-1]+v0*u[h-1];       for(s=1,h-2,           u[h]+=u[s]*v[h-1-s]       );       u[h]/=h    );    u }  fun_v(params, x)={    local(r);    if(x<0.25,        r=use_fun(0.0,0.0,params[1],x),    if(x<0.7,        r=use_fun(params[2],0.5,params[3],x),    if(x<1.4,        r=use_fun(params[4],1.0,params[5],x),    if(x<2.0,        r=use_fun(params[6],1.7,params[7],x),     if(x<2.4,        r=use_fun(params[8],2.2,params[9],x),    if(x<2.8,        r=use_fun(params[10],2.6,params[11],x),    if(x<3.2,        r=use_fun(params[12],3.0,params[13],x),    if(x<3.8,        r=use_fun(params[14],3.5,params[15],x),    if(x<4.5,        r=use_fun(params[16],4.2,params[17],x),        r=use_fun(params[18],5.0,params[19],x)    )))))))));    r }  fun_u(params, x)={    local(r);    if(x<0.7,        r=use_fun(params[1],0.5,params[2],x),    if(x<1.4,        r=use_fun(params[3],1.0,params[4],x),    if(x<2.0,        r=use_fun(params[5],1.7,params[6],x),     if(x<2.4,        r=use_fun(params[7],2.2,params[8],x),    if(x<2.8,        r=use_fun(params[9],2.6,params[10],x),    if(x<3.2,        r=use_fun(params[11],3.0,params[12],x),    if(x<3.8,        r=use_fun(params[13],3.5,params[14],x),    if(x<4.5,        r=use_fun(params[15],4.2,params[16],x),        r=use_fun(params[17],5.0,params[18],x)    ))))))));    r }

复制代码

mathe

上面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\),可以得到下面图形

mathe

自变量范围可限定在0到2pi-2d,d较小比如0.3

mathe

此楼数据由于使用v方程不对，所以结果不正确，正确结果查看53#.

现在可以计算出在d=0.590893319813969928977206512531098688032304589325 时取得最小期望a=3.55214056307013612938917610007896899689039684051
附件为对应Pari/Gp代码
在Pari/Gp下比如可以先输入
\p 100
设定计算精度为100位
然后用
\r a.gp
装入附件中gp文件
即可输出最小期望值
（还可以设定更高精度，不过更加花时间）。
代码中将v等函数分成很多段，是因为泰勒展开式的收敛半径不够大，只好这样表示