微分方程数值求解

mathe

我们首先能够想到的是公路可能的分布正好是单位圆的所有切线，每条切线出现的概率相同。我们可以采用切点对应的角度$theta$作为参数了表示这些切线。
第一步肯定是按照直线从圆心走到圆的边界，然后可以继续按照直线向前走若干距离（这段距离可以是0），然后拐弯按某个方向绕圆周行走，可以参考百度中帖子第49楼KeyTo9给出的图:

我们现在不严格分析可能的情况，先假设所用曲线总是这种类型（只是后面绕圆的曲线可以任意选择）。
我们将终点切线对应角度定义为角度0(也就是上面图中向上的方向为角度0),逆时针方向为角度增加的方向（所以行走过程中角度一直在降低）。
所以开始的时候，用户是沿着一个角度为-d的半径行走到圆外某处，知道到达一条参数为$2pi-2d$的切线（也就是对应切点的角度为$2pi-2d$）
而对于用户行走的线路上每个点,我们定义一个函数u为这个点到对应切点的距离。通过这种方法，用户行走线路的坐标可以表示成
  $x(\theta)=cos(\theta)-u(\theta)sin(\theta)$
  $y(\theta)=sin(\theta)+u(\theta)cos(\theta)$

第一种可以想到的解题方法是采用动态规划：
我们将$\theta$均匀划分成N份，同样对于u可能出现的值，在一定范围内划分成M份，
于是曲线上一段线的长度可以用折线长度
  $L(\theta_i)=sqrt((x(\theta_{i+1})-x(\theta_{i}))^2+(y(\theta_{i+1})-y(\theta_{i}))^2)$
来近似表示
而动态规划过程将要将表达式
  $\sum\theta_iL(\theta_i)$
最小化（这个部分不包含起始直线那一段）
上面过程显然可以通过动态规划来做

mathe

在将角度和切线都分成1000份左右的时候，通过计算可以得出百度链接中49楼相类似的结果。
不过我发现在将参数进一步提高时，发现最优结果的长度要变大，分别从3.5...提高到3.6...和3.7...
这个说明上面部分计算过程的计算误差还是很显著的（其中一个明显的原因是计算曲线的长度采用折线来代替，导致结果偏小）。这个结果同理论最优结果到底相差多少还是很难知道。
为此，我决定看看采用理论分析的方法是否可以得到精确的结果和采用的方法。

mathe

我们同样假设曲线还是类似上面的形式，我们现在先只对$0<=\theta<=2pi-2d$部分进行分析。
由
${ds}/{d\theta}=sqrt({{dx}/{d\theta}}^2+{{dy}/{d\theta}}^2)$
式子
$\sum\theta_iL(\theta_i)$
的连续表示方式就是
$\int_0^{2pi-2d} \theta {ds}/{d\theta} d\theta$ = $ \int_0^{2pi-2d} \theta sqrt(u^2+(1+u'^2) ) d\theta$
也就是我们要求上面积分的最小值，
然后利用变分法 （http://mathworld.wolfram.com/CalculusofVariations.html）就可以得到一楼中的微分方程
$\theta u^2u''=-(1+u')^3+2\theta u(1+u')^2-(u+\theta)u(1+u')-\theta u^3$

mathe

上面的微分方程是非线性的，很难去数值求解。而我也没有什么好的工具，所以这个微分方程让我相当头疼。
不过可以注意到在上面微分方程中，让$\theta=0$,我们马上可以得到$u'(0)=-1$
而两边同时经过一次微分后让后让$\theta=0$代入，我们可以得到$u''(0)=-{u(0)}/2$
同样，两边同时n次微分后，我们应该可以计算处$u^{(n+1)}(0)$的值（用$u(0)$来表示）
最终我采用了C语言进行符号运算，算出了上面式子的各阶微分到前50项，并且解出在$\theta=0$时各阶微分的值。
惊奇的发现所有奇数次导数在0的取值是常数，同$u(0)$没有关系，而偶数次导数在0的取值同$u(0)$成正比。
由此我们知道，函数u可以写成 $u(\theta)=w(\theta)+u(0)v(\theta)$
其中w是奇函数，v是偶函数。
而且只要我们分别计算出$u(0)=1$和$u(0)=2$时候的函数u的表达式，对于所有其他情况下的函数u都可以用这两个函数的线性组合表示出来

mathe

经计算，在$u(0)=1$和$u(0)=2$时，其在$\theta=0$的泰勒展开式的系数分别为
u(0)=1, u(0)=2
.1e1,.2e1
.-1e1,.-1e1
.-25e0,.-5e0
.166666666666666666667e0,.166666666666666666667e0
.546875e-1,.109375e0
.-354166666666666666667e-1,.-354166666666666666667e-1
.-130208333333333333333e-1,.-260416666666666666667e-1
.866815476190476190476e-2,.866815476190476190476e-2
.340016682942708333333e-2,.680033365885416666667e-2
.-232663432126322751323e-2,.-232663432126322751323e-2
.-948969523111979166667e-3,.-189793904622395833333e-2
.662781172729903198653e-3,.662781172729903198653e-3
.277329815758599175347e-3,.554659631517198350694e-3
.-196618716592471475284e-3,.-196618716592471475284e-3
.-837725509051022345787e-4,.-167545101810204469157e-3
.600682502332136375512e-4,.600682502332136375512e-4
.259416088808482077792e-4,.518832177616964155584e-4
.-187661782965473318626e-4,.-187661782965473318626e-4
.-819043055887853913995e-5,.-163808611177570782799e-4
.596710589150589134315e-5,.596710589150589134315e-5
.262646760239261690438e-5,.525293520478523380876e-5
.-19246660399963770934e-5,.-19246660399963770934e-5
上面所有数字中，小数点后面的负号要出现在小数点前面(gmp的输出功能中我没有找到好的方法，就这样输出了）
可以看出，这个系数减小的速度还是比较慢的，通过这个泰勒展开式如果计算$\theta$比较大的情况是很难的。
但是丢于$\theta$比较小的情况，使用这几项(20项)就应该可以计算出非常精确的值。
所以我们可以将区间$(0,2pi)$分成10份，上面的多项式仅仅用于第一份区间。然后我们将$\theta={pi}/5$代入多项式极其导数，算出在$u({pi}/5)$和$u('({pi}/5)$,然后我们再次使用14楼的类似方法，可以计算u在$\theta={pi}/5$的各阶导数，从而得到这里的泰勒展开式；继续下去，我们就可以通过这种方法得到函数w和v的分段表示（每一段都是一个高次多项式）

mathe

在得出这些参数以后，我们就可以计算对于不同的$u(0)$,对应的曲线$u(\theta)$的图像。计算结果表明，通常很快,在$\theta$稍微大一点时，对应的$u(\theta)$马上就小于0了，这个同几何意义是矛盾的。实际上这个表示在$\theta$以后的角度中，用户都应该贴着圆周走。这个结果让我非常吃惊，因为它得到的模型同动态规划的结果已经不同了。
最后我们还需要计算积分$ \int_0^{2pi-2d} \theta sqrt(u^2+(1+{u'}^2) ) d\theta$，如果采用数字积分的方法，不仅编程比较麻烦，而且精度比较难控制，为此，我再次将函数 $sqrt(u^2+(1+{u'}^2)$的各阶导数通过计算机计算出来，然后将它表示成多项式形式（同样分段表示），然后这个定积分也就变成多项式计算了

mathe

由于$u(theta)$总是在某个$theta$取0，我们解方程计算出这个$theta$,然后前面部分就沿着圆周走就可以了。
计算结果表明在$u(0)=6.434$时，期望值可以取到最小值3.83796896

无心人

你学学Standard ML语言吧
省的你用C复杂，用数学软件麻烦

mathe

我不认为一门新的语言能解决这个问题。

无心人

是
但能简化程序的复杂度