简介

在百度数学吧的东方角落和KeyTo9のFans提出过很多非常漂亮的数学题,比如随机游走问题 就是他提出的。
这里我们要讨论另外一个有人匿名提出而KeyTo9のFans最早给出不错解答的很有意思的难题:极地出逃问题原题在百度数学吧 ,在数学研发论坛 详细讨论了这个问题引出的一个微分方程。

极地出逃问题说:
一个人在冰天雪地里迷路了,已知离他一英里的地方有一条笔直的高速公路,但他完全失去了方向,不知道这条公路在哪个方向。如果天气恶劣,能见度只有一米(更精确的描述是能见度为0,也就是说只有到了公路才知道)。
问:采取怎样的方式才能尽快找到高速公路?尽快的定义是期望值最短,也就是说如果往不同方向出发,平均值最小。

这个问题应该是在2007年4月之前就被提出了,
KeyTo9のFans最早给出了一个精度相当不错得数值解 : 3.549260。
后来mathe尝试使用变分法,给出了一个最佳路线在极坐标下的微分方程\theta u^2u^{\prime\prime}=-(1+u^{\prime})^3+2\theta u(1+u^{\prime})^2-(u+\theta)u(1+u^{\prime})+\theta u^3。只是最初mathe给出的微分方程弄错了一个符号,大家一直未能通过这个方法求得合理得解答。
直到2015年,mathe重拾这个问题,最后经过复杂的分析,给这个微分方程找到一种比较合适的计算方法,得出极地出逃问题平均期望值的高精度结果为3.54925966919232488335307997505516088581017092471241443348333130503956135888576020826237。

动态规划求解

这个问题如果改为计算最差情况路径最短 ,那么就会相对比较容易,百度数学吧里面东方角落给出了这种情况的结果为1+\sqrt{3}+\frac{7\pi}6
如果考虑最坏情况下的最小值 MIN(MAX)
向任意方向直行一英里,然后开始沿以出发点为中心,一英里为
半径的圆走。走到3/4圈后开始沿切线方向走。这样,在最坏情况

下,需走2+3PI/2

我觉得不对,应该像24楼那样直行一英里多一点,再贴圆,最后切出,得到
 _
MIN(MAX)= 1 +√3 + 7π/6 < 2+3PI/2
minmax

对于这个期望路径最短问题,东方角落最先建议使用一种先直行再拐弯绕行 方案,并给出如下的示意图
ep1
假设最优结果所用曲线总是这种类型(只是后面绕圆的曲线可以任意选择),
我们将终点切线对应角度定义为角度0(也就是上面图中向上的方向为角度0),逆时针方向为角度增加的方向(所以顺时针行走过程中角度一直在变小)。
开始的时候,用户是沿着一个角度为-d的半径先行走到圆外某处,直到到达一条参数为2\pi-2d的切线(也就是对应切点的角度为2\pi-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_i L(\theta_i)
最小化(这个部分不包含起始直线那一段)。
百度贴吧中KeyTo9のFans最先给出一个近似解3.549260, 估计就是采用这种方法或类似方法进行计算的。
mathe说他重复采用这种方法在将角度和u都分成1000份左右时,可以得出类似精度的结果。但是他发现此后即使再细分,精度也无法提升了。

变分法分析

我们同样假设曲线还是类似上面的形式,我们现在先只对0\le\theta\le2\pi-2d部分进行分析。
由于\frac{ds}{d\theta}=\sqrt{(\frac{dx}{d\theta})^2+(\frac{dy}{d\theta})^2},
而目标函数\sum\theta_i L(\theta_i)的连续表示方式为
\int_0^{2\pi-2d} \theta \frac{ds}{d\theta} d\theta = \int_0^{2\pi-2d} \theta \sqrt{u^2+(1+{u^{\prime}}^2) }d\theta,
即我们要求使得上面积分的取到最小值时的函数u(\theta)
mathe利用变分法 得到了下面的微分方程
\theta u^2u^{\prime\prime}=-(1+u^{\prime})^3+2\theta u(1+u^{\prime})^2-(u+\theta)u(1+u^{\prime})-\theta u^3.

尝试微分方程求解

上面的微分方程是非线性的,而且在\theta=0附近不是很稳定,很难数值求解。
不过可以注意到在上面微分方程中,让\theta=0,可以得到u^{\prime}(0)=-1.
等式两边同时经过一次微分后让后让\theta=0代入,我们可以继续得到u^{\prime\prime}(0)=-\frac{u(0)}2.
同样,两边同时n次微分后,我们可以计算出u^{(n+1)}(0)的值(用u(0)来表示)。
mathe采用了C语言进行符号运算,算出了上面式子的各阶微分到前50项,并且解出在\theta=0时各阶微分的值。
惊奇的发现所有奇数次导数在\theta=0的取值是常数,同u(0)没有关系,而偶数次导数在\theta=0的取值同u(0)成正比。
由此我们知道,函数u可以写成 u(\theta)=w(\theta)+u(0)v(\theta),其中w是奇函数,v是偶函数。
于是我们只要分别计算出u(0)=1u(0)=2时候的函数u的表达式,对于所有其他情况下的函数u都可以用这两个函数的线性组合表示出来。
但是mathe采用上面的思路进行计算以后,发现得出的结果明显不合理
在得出这些参数以后,我们就可以计算对于不同的u(0),对应的曲线u(\theta)的图像。计算结果表明,通常很快,在\theta稍微大一点时,对应的u(\theta)马上就小于0了,这个同几何意义是矛盾的。实际上这个表示在\theta以后的角度中,用户都应该贴着圆周走。这个结果让mathe非常吃惊,因为它得到的模型同动态规划的结果已经不同了。
其中需要计算积分\int_0^{2\pi-2d} \theta \sqrt{u^2+(1+(u^{\prime})^2) } d\theta,如果采用数字积分的方法,不仅编程比较麻烦,而且精度比较难控制。为此,他再次将函数 \sqrt{u^2+(1+(u^{\prime})^2}的各阶导数通过计算机计算出来,然后将它表示成多项式形式(同样分段表示),然后这个定积分也就变成多项式计算了.

当时mathe猜测上面函数的问题是只含有一个参变量,而通常二阶微分应该包含两个变参(这样我们分别给定两个边界条件就可以得到一个唯一解)。
产生上面的原因是解方程过程中假设了函数u(\theta)\theta=0任意阶导数都存在。可能另外存在一些解的二阶导数趋向无穷了,所以这个方法就计算不出来了。
根据这个结果,得出的结论是最优的线路是先一直在圆周边界上走到最后只余下83度的角度,然后切换到对应的u(0)=1.714的曲线上,总共期望距离为4.0749987,远远大于前面动态规划得出的结果。

qxqcn曾经请一位精通Mathematica/Maple的朋友帮忙解决这个微分方程,没有成功。后来wayne试着使用Mathematica9.0.1同样软件也报错,他分析原因发现是由于初始值不当导致的。
尝试用无穷小代替0,得出了如下数值解
ep2
对应的Mathematica代码为

U=Block[{\[Epsilon]=$MachineEpsilon},NDSolveValue[{t u[t]^2u''[t]==-(1+u'[t])^3+2t u[t](1+u'[t])^2-(u[t]+t)u[t](1+u'[t])-t u[t]^3,u[\[Epsilon]]==1,u'[\[Epsilon]]==-1},u,{t,30}]]

mathe卷土重来

2015年1月,mathe重新分析计算这个问题 ,这时他发现前面计算中出了一点错误,其中最后一项符号弄错了,-\theta u^3应该改为+\theta u^3, 推导方法如下:
F(\theta,u,u^{\prime})=\theta \sqrt{u^2+(1+u^{\prime})^2},
要求\frac{\partial F}{\partial u} =\frac{d}{d\theta}\frac{\partial F}{\partial u^{\prime}},
注意最后一次对\theta求导时u,u^{\prime}都要看成\theta的函数,而u^{\prime}就是\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^{\prime}(0)=-1这个表明在曲线末端曲率半径就是u(0),也就是这个端点到对应圆的切点的距离。这个说明曲线末端是垂直于对应的切线的。
然后我们做变量替换v=\frac{1+u^{\prime}}{u},得出微分方程\frac{dv}{d\theta}=(1-\frac v{\theta})(v^2+1)(*), 边界条件v(0)=0.
如果再做替换v=\cot(\phi),得到\frac{d\phi}{d\theta}=\frac{\cot(\phi)}{\theta}-1,其中几何意义中\phi是轨迹上一点的切线和这个点对应的到单位圆(我们的目标圆)的切线的夹角.
上面v的微分方程和边界条件唯一确定函数v.
然后由于u^{\prime}-vu+1=0,假设方程有两个不同的解u_1,u_2,得出(u^{\prime}_1-u^{\prime}_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。
每个u的不同的通式将对应一个不同的d,我们需要从中再挑选出使得总平均距离最小的一条,估计计算量会比较大。但是由于前面已经有个东方角落的数值解,我们可以利用那个解来估算这些参数,应该会比较容易找到一个较好的解。

(*)这步推导过程mathe再次犯了一个错误,多亏了计算高手数学星空的指正:
\theta u^2u^{\prime\prime}=-(1+u^{\prime})^3+2\theta u(1+u^{\prime})^2-(u+\theta)u(1+u^{\prime})+\theta u^3
做变量替换v=\frac{1+u^{\prime}}{u},得出微分方程\frac{dv}{d\theta}=(\theta-v)(v^2+1) 应改成
\frac{dv}{d\theta}=(1-\frac{v}{\theta})(v^2+1).

v的不同解在v离1比较远时的图像是类似的,但是在0附近区别比较大。
而满足本问题的解应该有v(0)=0,v^{\prime}(0)=\frac12,而且可以看出这时v是一个奇函数。
设其在0的展开式为v(\theta)=\sum_{n=0}^{+\infty} a_{2n+1}\theta^{2n+1},我们得出:
\sum_{n=0}^{+\infty}(2n+1)a_{2n+1}\theta^{2n+1}=(\theta-\sum_{n=0}^{+\infty} a_{2n+1}\theta^{2n+1})(1+(\sum_{n=0}^{+\infty} a_{2n+1}\theta^{2n+1})^2).
展开比较系数就可以得出a_1=\frac12,(2n+2)a_{2n+1}=\sum_{k=0}^{n-1}a_{2k+1}a_{2n-2k-1}-\sum_{k=0}^{n-1}a_{2k+1}\sum_{m=0}^{n-k-1}a_{2m+1}a_{2(n-k-m)-1},n\ge 1.
而对应的方程u^{\prime}-vu+1=0初始条件u_0(0)=0的解也是奇函数。而V(\theta)=\exp(\int_0^{\theta} v(t)dt)是偶函数,u的通解为u(\theta)=u_0(\theta)+C\times V(\theta)

解决了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(2\pi-2d)=\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^{\prime}_0(2\pi-2d))+2(\tan(d)-u_0(2\pi-2d))V^{\prime}(2\pi-2d)}{V^2(2\pi-2d)}D(2\pi-2d)=0

最终的结果

更新代码以后,mathe得出结果如下
d = 0.59090489185668802875786709583348871021050937452181485289441819749025679734298579253579
a = 3.54925966919232488335307997505516088581017092471241443348333130503956135888576020826237
这个结果同fans的数据一致了
其中对应的图片
v
这是v(\theta)的图像

u
这是u(\theta)的图像

dT-d
这是这是目标函数关于d的导数的图像

T-d
这是目标函数关于d的图像,也就是距离的期望值

result
这个是最优解对应的图像

点击可以下载对应的Pari/gp代码

此外,v,u函数在\theta=0的收敛半径好像大于\pi,所以实际上,分段数目可以少很多,但是由于直接从以前代码修改过来,就保留了更多的分段了。
比较奇怪的是计算出来在\theta=0.5的收敛半径好像小于1,这个有点不符合道理(圆的边界上必然有函数的奇点,但是落在0为中心收敛圆内部)。
一个可能的解释就是计算误差引起的。可以看到满足v的微分方程的一族曲线中,除了这条v(0)=0的曲线,其余的在靠近\theta=0时,值都会发生剧烈变化,远远离开原点。
这个导致咋\theta=0.5时,如果有一点计算误差,对应到\theta=0处,会有很大的误差,所以就远离我们的目标曲线了。所幸,\theta越大,这种误差
影响越小。
另外,如果我们能够修改一下代码,更大范围使用在\theta=0处的展开式(可以有精确表达形式),会有更好的精度。
另外比较有意思的是取u_0(0)=0对应的u的解在\theta=0处展开式的前面15项竟然全部是正数(对应到\theta^{31}系数),但是可惜后面的就出现负数了。