极地出逃问题

简介

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

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

这个问题应该是在2007年4月之前就被提出了, KeyTo9のFans最早给出了一个精度相当不错得数值解 : 3.549260。 后来mathe尝试使用变分法,给出了一个最佳路线在极坐标下的微分方程θu2u=(1+u)3+2θu(1+u)2(u+θ)u(1+u)+θu3\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+3+7π61+\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-d的半径先行走到圆外某处,直到到达一条参数为2π2d2\pi-2d的切线(也就是对应切点的角度为2π2d2\pi-2d)。 对于用户行走的线路上每个点,我们定义一个函数u为这个点到对应切点的距离。通过这种参数设置,用户行走线路的坐标可以表示成 x(θ)=cos(θ)u(θ)sin(θ)x(\theta)=\cos(\theta)-u(\theta)\sin(\theta) y(θ)=sin(θ)+u(θ)cos(θ)y(\theta)=\sin(\theta)+u(\theta)\cos(\theta) 我们首先可以想到的解题方法是采用动态规划: 将θ\theta均匀划分成N份,并且同样对于u可能出现的值,在一定范围内均匀划分成M份, 于是曲线上一段线的长度可以用折线长度 L(θi)=(x(θi+1)x(θi))2+(y(θi+1)y(θi))2L(\theta_i)=\sqrt{(x(\theta_{i+1})-x(\theta_{i}))^2+(y(\theta_{i+1})-y(\theta_{i}))^2} 来近似表示,而动态规划过程要将表达式 θiL(θi)\sum\theta_i L(\theta_i) 最小化(这个部分不包含起始直线那一段)。 百度贴吧中KeyTo9のFans最先给出一个近似解3.549260, 估计就是采用这种方法或类似方法进行计算的。 mathe说他重复采用这种方法在将角度和u都分成1000份左右时,可以得出类似精度的结果。但是他发现此后即使再细分,精度也无法提升了。

变分法分析

我们同样假设曲线还是类似上面的形式,我们现在先只对0θ2π2d0\le\theta\le2\pi-2d部分进行分析。 由于dsdθ=(dxdθ)2+(dydθ)2\frac{ds}{d\theta}=\sqrt{(\frac{dx}{d\theta})^2+(\frac{dy}{d\theta})^2}, 而目标函数θiL(θi)\sum\theta_i L(\theta_i)的连续表示方式为 02π2dθdsdθdθ=02π2dθu2+(1+u2)dθ\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(θ)u(\theta)。 mathe利用变分法 得到了下面的微分方程 θu2u=(1+u)3+2θu(1+u)2(u+θ)u(1+u)θu3\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.

尝试微分方程求解

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

当时mathe猜测上面函数的问题是只含有一个参变量,而通常二阶微分应该包含两个变参(这样我们分别给定两个边界条件就可以得到一个唯一解)。 产生上面的原因是解方程过程中假设了函数u(θ)u(\theta)θ=0\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重新分析计算这个问题 ,这时他发现前面计算中出了一点错误,其中最后一项符号弄错了,θu3-\theta u^3应该改为+θu3+\theta u^3, 推导方法如下: F(θ,u,u)=θu2+(1+u)2F(\theta,u,u^{\prime})=\theta \sqrt{u^2+(1+u^{\prime})^2}, 要求Fu=ddθFu\frac{\partial F}{\partial u} =\frac{d}{d\theta}\frac{\partial F}{\partial u^{\prime}}, 注意最后一次对θ\theta求导时u,uu,u^{\prime}都要看成θ\theta的函数,而uu^{\prime}就是dudθ\frac{du}{d\theta} 而已知边界条件为u(2π2d)=tan(d)u(2\pi-2d)=\tan(d) 而求出函数u(θ)u(\theta)后,我们只要画出参数曲线(cos(θ)u(θ)sin(θ),sin(θ)+u(θ)cos(θ))(\cos(\theta)-u(\theta)\sin(\theta),\sin(\theta)+u(\theta)\cos(\theta)) 另外根据方程我们容易看出u(0)=1u^{\prime}(0)=-1这个表明在曲线末端曲率半径就是u(0)u(0),也就是这个端点到对应圆的切点的距离。这个说明曲线末端是垂直于对应的切线的。 然后我们做变量替换v=1+uuv=\frac{1+u^{\prime}}{u},得出微分方程dvdθ=(1vθ)(v2+1)\frac{dv}{d\theta}=(1-\frac v{\theta})(v^2+1)(*), 边界条件v(0)=0v(0)=0. 如果再做替换v=cot(ϕ)v=\cot(\phi),得到dϕdθ=cot(ϕ)θ1\frac{d\phi}{d\theta}=\frac{\cot(\phi)}{\theta}-1,其中几何意义中ϕ\phi是轨迹上一点的切线和这个点对应的到单位圆(我们的目标圆)的切线的夹角. 上面vv的微分方程和边界条件唯一确定函数vv. 然后由于uvu+1=0u^{\prime}-vu+1=0,假设方程有两个不同的解u1,u2u_1,u_2,得出(u1u2)=v(u1v2)(u^{\prime}_1-u^{\prime}_2)=v(u_1-v_2)所以u1u2=Cexp(0θv(x)dx)u_1-u_2=C\exp(\int_0^{\theta}v(x)dx) 也就是任意算出一个满足条件的函数uu以后,只要加上V(θ)=exp(0θv(x)dx)V(\theta)=\exp(\int_0^{\theta}v(x)dx)的常数倍,就得出uu的通式,最后我们需要通过方程u(2π2d)=tan(d)u(2\pi-2d)=\tan(d)计算出对应的d。 每个uu的不同的通式将对应一个不同的d,我们需要从中再挑选出使得总平均距离最小的一条,估计计算量会比较大。但是由于前面已经有个东方角落的数值解,我们可以利用那个解来估算这些参数,应该会比较容易找到一个较好的解。

(*)这步推导过程mathe再次犯了一个错误,多亏了计算高手数学星空的指正: θu2u=(1+u)3+2θu(1+u)2(u+θ)u(1+u)+θu3\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=1+uuv=\frac{1+u^{\prime}}{u},得出微分方程dvdθ=(θv)(v2+1)\frac{dv}{d\theta}=(\theta-v)(v^2+1) 应改成 dvdθ=(1vθ)(v2+1)\frac{dv}{d\theta}=(1-\frac{v}{\theta})(v^2+1).

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

解决了vv的展开式以后,我们可以用uu的一个特解u0u_0的展开式以及VV的展开式计算其通解。 此后我们知道目标函数为 A(d)+02π2dθu(θ)v2(θ)+1dθA(d)+\int_0^{2\pi-2d} \theta u(\theta)\sqrt{v^2(\theta)+1}d\theta其中A(d)=ln(1+sin(d)1sin(d))+2π2dcos(d)A(d)=\ln(\frac{1+\sin(d)}{1-\sin(d)})+\frac{2\pi-2d}{\cos(d)} A(d)+02π2dθu0(θ)v2(θ)+1dθ+C02π2dθV(θ)v2(θ)+1dθ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)=0xθu0(θ)v2(θ)+1dθ,D(x)=0xθV(θ)v2(θ)+1dθ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π2d)=tan(d)u(2\pi-2d)=\tan(d)可以得出C=tan(d)u0(2π2d)V(2π2d)C=\frac{\tan(d)-u_0(2\pi-2d)}{V(2\pi-2d)} 代入以后目标函数变成只含一个变量d的函数average(d)=A(d)+B(2π2d)+tan(d)u0(2π2d)V(2π2d)D(2π2d)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。 也就是 (2π2d)sin(d)cos2(d)2(2π2d)u0(2π2d)v2(2π2d)+12tan(d)u0(2π2d)V(2π2d)(2π2d)V(2π2d)v2(2π2d)+1+V(2π2d)(sec2(d)+2u0(2π2d))+2(tan(d)u0(2π2d))V(2π2d)V2(2π2d)D(2π2d)=0\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(θ)v(\theta)的图像

u 这是$u(\theta)$的图像 dT-d 这是这是目标函数关于d的导数的图像 T-d 这是目标函数关于d的图像,也就是距离的期望值 result 这个是最优解对应的图像

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

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

© 2021 Blog Boost Starter