简介
在百度数学吧的东方角落和KeyTo9のFans提出过很多非常漂亮的数学题，比如随机游走问题 就是他提出的。
这里我们要讨论另外一个有人匿名提出而KeyTo9のFans最早给出不错解答的很有意思的难题：极地出逃问题 。原题在百度数学吧 ，在数学研发论坛 详细讨论了这个问题引出的一个微分方程。
极地出逃问题说：
这个问题应该是在2007年4月之前就被提出了，
KeyTo9のFans最早给出了一个精度相当不错得数值解 : 3.549260。θ u 2 u ′ ′ = − ( 1 + u ′ ) 3 + 2 θ u ( 1 + u ′ ) 2 − ( u + θ ) u ( 1 + u ′ ) + θ u 3 \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 θ u 2 u ′ ′ = − ( 1 + u ′ ) 3 + 2 θ u ( 1 + u ′ ) 2 − ( u + θ ) u ( 1 + u ′ ) + θ u 3
动态规划求解
这个问题如果改为计算最差情况路径最短 ，那么就会相对比较容易，百度数学吧里面东方角落给出了这种情况的结果为1 + 3 + 7 π 6 1+\sqrt{3}+\frac{7\pi}6 1 + 3 + 6 7 π
我觉得不对，应该像24楼那样直行一英里多一点，再贴圆，最后切出，得到
_
ＭＩＮ（ＭＡＸ）= 1 +√3 + 7π/6 ＜ ２＋３ＰＩ／2
对于这个期望路径最短问题，东方角落最先建议使用一种先直行再拐弯绕行 方案，并给出如下的示意图− d -d − d 2 π − 2 d 2\pi-2d 2 π − 2 d 2 π − 2 d 2\pi-2d 2 π − 2 d x ( θ ) = cos ( θ ) − u ( θ ) sin ( θ ) x(\theta)=\cos(\theta)-u(\theta)\sin(\theta) x ( θ ) = cos ( θ ) − u ( θ ) sin ( θ ) y ( θ ) = sin ( θ ) + u ( θ ) cos ( θ ) y(\theta)=\sin(\theta)+u(\theta)\cos(\theta) y ( θ ) = sin ( θ ) + u ( θ ) cos ( θ ) θ \theta θ L ( θ i ) = ( x ( θ i + 1 ) − x ( θ i ) ) 2 + ( y ( θ i + 1 ) − y ( θ i ) ) 2 L(\theta_i)=\sqrt{(x(\theta_{i+1})-x(\theta_{i}))^2+(y(\theta_{i+1})-y(\theta_{i}))^2} L ( θ i ) = ( x ( θ i + 1 ) − x ( θ i ) ) 2 + ( y ( θ i + 1 ) − y ( θ i ) ) 2 ∑ θ i L ( θ i ) \sum\theta_i L(\theta_i) ∑ θ i L ( θ i )
变分法分析
我们同样假设曲线还是类似上面的形式，我们现在先只对0 ≤ θ ≤ 2 π − 2 d 0\le\theta\le2\pi-2d 0 ≤ θ ≤ 2 π − 2 d d s d θ = ( d x d θ ) 2 + ( d y d θ ) 2 \frac{ds}{d\theta}=\sqrt{(\frac{dx}{d\theta})^2+(\frac{dy}{d\theta})^2} d θ d s = ( d θ d x ) 2 + ( d θ d y ) 2 ∑ θ i L ( θ i ) \sum\theta_i L(\theta_i) ∑ θ i L ( θ i ) ∫ 0 2 π − 2 d θ d s d θ d θ = ∫ 0 2 π − 2 d θ u 2 + ( 1 + u ′ 2 ) 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 ∫ 0 2 π − 2 d θ d θ d s d θ = ∫ 0 2 π − 2 d θ u 2 + ( 1 + u ′ 2 ) d θ u ( θ ) u(\theta) u ( θ ) 变分法 得到了下面的微分方程θ u 2 u ′ ′ = − ( 1 + u ′ ) 3 + 2 θ u ( 1 + u ′ ) 2 − ( u + θ ) u ( 1 + u ′ ) − θ u 3 \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 θ u 2 u ′ ′ = − ( 1 + u ′ ) 3 + 2 θ u ( 1 + u ′ ) 2 − ( u + θ ) u ( 1 + u ′ ) − θ u 3
尝试微分方程求解
上面的微分方程是非线性的，而且在θ = 0 \theta=0 θ = 0 θ = 0 \theta=0 θ = 0 u ′ ( 0 ) = − 1 u^{\prime}(0)=-1 u ′ ( 0 ) = − 1 θ = 0 \theta=0 θ = 0 u ′ ′ ( 0 ) = − u ( 0 ) 2 u^{\prime\prime}(0)=-\frac{u(0)}2 u ′ ′ ( 0 ) = − 2 u ( 0 ) u ( n + 1 ) ( 0 ) u^{(n+1)}(0) u ( n + 1 ) ( 0 ) u ( 0 ) u(0) u ( 0 ) θ = 0 \theta=0 θ = 0 θ = 0 \theta=0 θ = 0 u ( 0 ) u(0) u ( 0 ) θ = 0 \theta=0 θ = 0 u ( 0 ) u(0) u ( 0 ) u ( θ ) = w ( θ ) + u ( 0 ) v ( θ ) u(\theta)=w(\theta)+u(0)v(\theta) u ( θ ) = w ( θ ) + u ( 0 ) v ( θ ) u ( 0 ) = 1 u(0)=1 u ( 0 ) = 1 u ( 0 ) = 2 u(0)=2 u ( 0 ) = 2 结果明显不合理 。u ( 0 ) u(0) u ( 0 ) u ( θ ) u(\theta) u ( θ ) θ \theta θ u ( θ ) u(\theta) u ( θ ) θ \theta θ ∫ 0 2 π − 2 d θ u 2 + ( 1 + ( u ′ ) 2 ) d θ \int_0^{2\pi-2d} \theta \sqrt{u^2+(1+(u^{\prime})^2) } d\theta ∫ 0 2 π − 2 d θ u 2 + ( 1 + ( u ′ ) 2 ) d θ u 2 + ( 1 + ( u ′ ) 2 \sqrt{u^2+(1+(u^{\prime})^2} u 2 + ( 1 + ( u ′ ) 2
当时mathe猜测上面函数的问题是只含有一个参变量，而通常二阶微分应该包含两个变参（这样我们分别给定两个边界条件就可以得到一个唯一解）。
产生上面的原因是解方程过程中假设了函数u ( θ ) u(\theta) u ( θ ) θ = 0 \theta=0 θ = 0
qxqcn曾经请一位精通Mathematica/Maple的朋友帮忙解决这个微分方程，没有成功。后来wayne试着使用Mathematica9.0.1同样软件也报错，他分析原因发现是由于初始值不当导致的。
他尝试用无穷小代替0，得出了如下数值解
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重新分析计算这个问题 ，这时他发现前面计算中出了一点错误，其中最后一项符号弄错了，− θ u 3 -\theta u^3 − θ u 3 + θ u 3 +\theta u^3 + θ u 3 F ( θ , u , u ′ ) = θ u 2 + ( 1 + u ′ ) 2 F(\theta,u,u^{\prime})=\theta \sqrt{u^2+(1+u^{\prime})^2} F ( θ , u , u ′ ) = θ u 2 + ( 1 + u ′ ) 2 ∂ F ∂ u = d d θ ∂ F ∂ u ′ \frac{\partial F}{\partial u} =\frac{d}{d\theta}\frac{\partial F}{\partial u^{\prime}} ∂ u ∂ F = d θ d ∂ u ′ ∂ F θ \theta θ u , u ′ u,u^{\prime} u , u ′ θ \theta θ u ′ u^{\prime} u ′ d u d θ \frac{du}{d\theta} d θ d u u ( 2 π − 2 d ) = tan ( d ) u(2\pi-2d)=\tan(d) u ( 2 π − 2 d ) = tan ( d ) u ( θ ) u(\theta) u ( θ ) ( cos ( θ ) − u ( θ ) sin ( θ ) , sin ( θ ) + u ( θ ) cos ( θ ) ) (\cos(\theta)-u(\theta)\sin(\theta),\sin(\theta)+u(\theta)\cos(\theta)) ( cos ( θ ) − u ( θ ) sin ( θ ) , sin ( θ ) + u ( θ ) cos ( θ ) ) u ′ ( 0 ) = − 1 u^{\prime}(0)=-1 u ′ ( 0 ) = − 1 u ( 0 ) u(0) u ( 0 ) v = 1 + u ′ u v=\frac{1+u^{\prime}}{u} v = u 1 + u ′ d v d θ = ( 1 − v θ ) ( v 2 + 1 ) \frac{dv}{d\theta}=(1-\frac v{\theta})(v^2+1) d θ d v = ( 1 − θ v ) ( v 2 + 1 ) v ( 0 ) = 0 v(0)=0 v ( 0 ) = 0 v = cot ( ϕ ) v=\cot(\phi) v = cot ( ϕ ) d ϕ d θ = cot ( ϕ ) θ − 1 \frac{d\phi}{d\theta}=\frac{\cot(\phi)}{\theta}-1 d θ d ϕ = θ c o t ( ϕ ) − 1 ϕ \phi ϕ v v v v v v u ′ − v u + 1 = 0 u^{\prime}-vu+1=0 u ′ − v u + 1 = 0 u 1 , u 2 u_1,u_2 u 1 , u 2 ( u 1 ′ − u 2 ′ ) = v ( u 1 − v 2 ) (u^{\prime}_1-u^{\prime}_2)=v(u_1-v_2) ( u 1 ′ − u 2 ′ ) = v ( u 1 − v 2 ) u 1 − u 2 = C exp ( ∫ 0 θ v ( x ) d x ) u_1-u_2=C\exp(\int_0^{\theta}v(x)dx) u 1 − u 2 = C exp ( ∫ 0 θ v ( x ) d x ) u u u V ( θ ) = exp ( ∫ 0 θ v ( x ) d x ) V(\theta)=\exp(\int_0^{\theta}v(x)dx) V ( θ ) = exp ( ∫ 0 θ v ( x ) d x ) u u u u ( 2 π − 2 d ) = tan ( d ) u(2\pi-2d)=\tan(d) u ( 2 π − 2 d ) = tan ( d ) u u u
（*）这步推导过程mathe再次犯了一个错误，多亏了计算高手数学星空的指正:θ u 2 u ′ ′ = − ( 1 + u ′ ) 3 + 2 θ u ( 1 + u ′ ) 2 − ( u + θ ) u ( 1 + u ′ ) + θ u 3 \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 θ u 2 u ′ ′ = − ( 1 + u ′ ) 3 + 2 θ u ( 1 + u ′ ) 2 − ( u + θ ) u ( 1 + u ′ ) + θ u 3 v = 1 + u ′ u v=\frac{1+u^{\prime}}{u} v = u 1 + u ′ d v d θ = ( θ − v ) ( v 2 + 1 ) \frac{dv}{d\theta}=(\theta-v)(v^2+1) d θ d v = ( θ − v ) ( v 2 + 1 ) d v d θ = ( 1 − v θ ) ( v 2 + 1 ) \frac{dv}{d\theta}=(1-\frac{v}{\theta})(v^2+1) d θ d v = ( 1 − θ v ) ( v 2 + 1 )
v的不同解在v离1比较远时的图像是类似的，但是在0附近区别比较大。v ( 0 ) = 0 , v ′ ( 0 ) = 1 2 v(0)=0,v^{\prime}(0)=\frac12 v ( 0 ) = 0 , v ′ ( 0 ) = 2 1 v v v v ( θ ) = ∑ n = 0 + ∞ a 2 n + 1 θ 2 n + 1 v(\theta)=\sum_{n=0}^{+\infty} a_{2n+1}\theta^{2n+1} v ( θ ) = ∑ n = 0 + ∞ a 2 n + 1 θ 2 n + 1 ∑ n = 0 + ∞ ( 2 n + 1 ) a 2 n + 1 θ 2 n + 1 = ( θ − ∑ n = 0 + ∞ a 2 n + 1 θ 2 n + 1 ) ( 1 + ( ∑ n = 0 + ∞ a 2 n + 1 θ 2 n + 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) ∑ n = 0 + ∞ ( 2 n + 1 ) a 2 n + 1 θ 2 n + 1 = ( θ − ∑ n = 0 + ∞ a 2 n + 1 θ 2 n + 1 ) ( 1 + ( ∑ n = 0 + ∞ a 2 n + 1 θ 2 n + 1 ) 2 ) a 1 = 1 2 , ( 2 n + 2 ) a 2 n + 1 = ∑ k = 0 n − 1 a 2 k + 1 a 2 n − 2 k − 1 − ∑ k = 0 n − 1 a 2 k + 1 ∑ m = 0 n − k − 1 a 2 m + 1 a 2 ( n − k − m ) − 1 , n ≥ 1 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 a 1 = 2 1 , ( 2 n + 2 ) a 2 n + 1 = ∑ k = 0 n − 1 a 2 k + 1 a 2 n − 2 k − 1 − ∑ k = 0 n − 1 a 2 k + 1 ∑ m = 0 n − k − 1 a 2 m + 1 a 2 ( n − k − m ) − 1 , n ≥ 1 u ′ − v u + 1 = 0 u^{\prime}-vu+1=0 u ′ − v u + 1 = 0 u 0 ( 0 ) = 0 u_0(0)=0 u 0 ( 0 ) = 0 V ( θ ) = exp ( ∫ 0 θ v ( t ) d t ) V(\theta)=\exp(\int_0^{\theta} v(t)dt) V ( θ ) = exp ( ∫ 0 θ v ( t ) d t ) u ( θ ) = u 0 ( θ ) + C × V ( θ ) u(\theta)=u_0(\theta)+C\times V(\theta) u ( θ ) = u 0 ( θ ) + C × V ( θ )
解决了v v v u u u u 0 u_0 u 0 V V V A ( d ) + ∫ 0 2 π − 2 d θ u ( θ ) v 2 ( θ ) + 1 d θ A(d)+\int_0^{2\pi-2d} \theta u(\theta)\sqrt{v^2(\theta)+1}d\theta A ( d ) + ∫ 0 2 π − 2 d θ u ( θ ) v 2 ( θ ) + 1 d θ A ( d ) = ln ( 1 + sin ( d ) 1 − sin ( d ) ) + 2 π − 2 d cos ( d ) A(d)=\ln(\frac{1+\sin(d)}{1-\sin(d)})+\frac{2\pi-2d}{\cos(d)} A ( d ) = ln ( 1 − s i n ( d ) 1 + s i n ( d ) ) + c o s ( d ) 2 π − 2 d A ( d ) + ∫ 0 2 π − 2 d θ u 0 ( θ ) v 2 ( θ ) + 1 d θ + C ∫ 0 2 π − 2 d θ V ( θ ) v 2 ( θ ) + 1 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 A ( d ) + ∫ 0 2 π − 2 d θ u 0 ( θ ) v 2 ( θ ) + 1 d θ + C ∫ 0 2 π − 2 d θ V ( θ ) v 2 ( θ ) + 1 d θ B ( x ) = ∫ 0 x θ u 0 ( θ ) v 2 ( θ ) + 1 d θ , D ( x ) = ∫ 0 x θ V ( θ ) v 2 ( θ ) + 1 d θ 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 ( x ) = ∫ 0 x θ u 0 ( θ ) v 2 ( θ ) + 1 d θ , D ( x ) = ∫ 0 x θ V ( θ ) v 2 ( θ ) + 1 d θ u ( 2 π − 2 d ) = tan ( d ) u(2\pi-2d)=\tan(d) u ( 2 π − 2 d ) = tan ( d ) C = tan ( d ) − u 0 ( 2 π − 2 d ) V ( 2 π − 2 d ) C=\frac{\tan(d)-u_0(2\pi-2d)}{V(2\pi-2d)} C = V ( 2 π − 2 d ) t a n ( d ) − u 0 ( 2 π − 2 d ) a v e r a g e ( d ) = A ( d ) + B ( 2 π − 2 d ) + tan ( d ) − u 0 ( 2 π − 2 d ) V ( 2 π − 2 d ) D ( 2 π − 2 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) a v e r a g e ( d ) = A ( d ) + B ( 2 π − 2 d ) + V ( 2 π − 2 d ) t a n ( d ) − u 0 ( 2 π − 2 d ) D ( 2 π − 2 d ) ( 2 π − 2 d ) sin ( d ) cos 2 ( d ) − 2 ( 2 π − 2 d ) u 0 ( 2 π − 2 d ) v 2 ( 2 π − 2 d ) + 1 − 2 tan ( d ) − u 0 ( 2 π − 2 d ) V ( 2 π − 2 d ) ( 2 π − 2 d ) V ( 2 π − 2 d ) v 2 ( 2 π − 2 d ) + 1 + V ( 2 π − 2 d ) ( sec 2 ( d ) + 2 u 0 ′ ( 2 π − 2 d ) ) + 2 ( tan ( d ) − u 0 ( 2 π − 2 d ) ) V ′ ( 2 π − 2 d ) V 2 ( 2 π − 2 d ) D ( 2 π − 2 d ) = 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 c o s 2 ( d ) ( 2 π − 2 d ) s i n ( d ) − 2 ( 2 π − 2 d ) u 0 ( 2 π − 2 d ) v 2 ( 2 π − 2 d ) + 1 − 2 V ( 2 π − 2 d ) t a n ( d ) − u 0 ( 2 π − 2 d ) ( 2 π − 2 d ) V ( 2 π − 2 d ) v 2 ( 2 π − 2 d ) + 1 + V 2 ( 2 π − 2 d ) V ( 2 π − 2 d ) ( s e c 2 ( d ) + 2 u 0 ′ ( 2 π − 2 d ) ) + 2 ( t a n ( d ) − u 0 ( 2 π − 2 d ) ) V ′ ( 2 π − 2 d ) D ( 2 π − 2 d ) = 0
最终的结果
更新代码以后，mathe得出结果如下v ( θ ) v(\theta) v ( θ )
u ( θ ) u(\theta) u ( θ )
点击可以下载对应的Pari/gp代码
此外，v,u函数在θ = 0 \theta=0 θ = 0 π \pi π θ = 0.5 \theta=0.5 θ = 0 . 5 θ = 0 \theta=0 θ = 0 θ = 0.5 \theta=0.5 θ = 0 . 5 θ = 0 \theta=0 θ = 0 θ \theta θ θ = 0 \theta=0 θ = 0 u 0 ( 0 ) = 0 u_0(0)=0 u 0 ( 0 ) = 0 θ = 0 \theta=0 θ = 0 θ 31 \theta^{31} θ 3 1
