Brachistochrone curve(传说中的最速降线)

  Brachistochrone 来自于古希腊语：βράχιστος χρόνος (brakhistos, superlative of brakhus short + chronos time)，意思是时间最短. Brachistochrone curve 指的就是"时间最短"的曲线.
  关于最速降线问题的起源和发展参见：知乎贴，用一句话来概括就是：神仙打架！

  为了您更好的阅读体验，请使用电脑浏览.

一、问题重述

  一质量为 $m$ 的质点，在重力作用下从定点 $A$ 沿曲线下滑至定点 $B$，试确定一条曲线，使得质点由 $A$ 下滑到 $B$ 所需时间最短.
  假定 $B$ 比 $A$ 低，重力加速度为常数 $g$，不计摩擦力和其他阻力等因素.

二、数学描述

  建立坐标系，如图所示. 设曲线为 $y=y(x),(x\in [0,c])$，满足 $y(0)=0,y(c)=H$，且足够光滑，并设这样的函数组成的全体为集合 $C$，我们的目标是求函数 $y$，使得
y = arg min ⁡ C {   质点沿   y = y ( x )   由   A   到   B   的下滑时间   } y=\argmin _C \{\,\textbf{质点沿}\,y=y(x)\, \textbf{由}\,A\,\textbf{到}\,B\,\textbf{的下滑时间}\,\} y=Cargmin​{质点沿y=y(x)由A到B的下滑时间}

三、物理建模

一言一蔽之，曰：分割-近似-取极限.

分割：
  建立相同坐标系，设 $A(0,0),B(c,H)$，将带状区域 $0<y<H$ 用平行于 $x$ 轴的直线 $y=y_k=kH/n$ 分割成 $n$ 个小带状区域.

近似：
  在带状区域 $y_{k-1}<y<y_k$ 内，可近似认为 $v_k=\sqrt{2g{y}_k}$ (根据能量守恒定律可得)不变，近似认为曲线为直线.

分析：怎样的路线才是最快的？先来考虑下述问题.

  设质点从 $A_1$ 经直线 $l$ 到达 $A_2$，质点在 $l$ 上侧时速度为 $v_1$，下侧时速度为 $v_2$，请问质点应沿什么路线运动才最省时？

  因为质点在直线两侧时速度不变，所以运动轨迹应是折线，折点在直线 $l$ 上，问题转化为确定最佳折点的位置.

  设 $A_{1}O\perp l,A_{2}D\perp l,C$ 为折点，其余标注如图所示，只有 $x$ 为变量.

则质点由 $A_1$ 到 $A_2$ 所需时间为$$t=\frac{A_{1}C}{v_1}+\frac{C{A}_2}{v_2}=\frac{\sqrt{x^2+a^2}}{v_1}+\frac{\sqrt{(c-x{)}^2+b^2}}{v_2}$$对 $x$ 求导$$\frac{dt}{dx}=\frac{1}{v_1}\frac{2x}{\sqrt{x^2+a^2}}-\frac{1}{v_2}\frac{2(c-x)}{\sqrt{(c-x{)}^2+b^2}}$$唯一驻点满足$$\frac{dt}{dx}=0\Rightarrow \frac{1}{v_1}\frac{x}{\sqrt{x^2+a^2}}=\frac{1}{v_2}\frac{(c-x)}{\sqrt{(c-x{)}^2+b^2}}$$即$$\frac{\sin {\alpha }_1}{v_1}=\frac{\sin {\alpha }_2}{v_2}$$这就是光学中的 $\text{Snell}$ 折射定律，即光沿"最短路径"传播.

回过头来继续看我们的问题

  考虑质点经过第 $k$ 层与第 $k+1$ 层，根据近似，质点在每层中速度不变，由上述结论知，要使时间最短，${\alpha }_k,{\alpha }_{k+1}$ 应满足$$\frac{\sin {\alpha }_k}{v_k}=\frac{\sin {\alpha }_{k+1}}{v_{k+1}}$$对任意的 $k$ 成立，则 $$\frac{\sin {\alpha }_k}{v_k}=C_1(\text{常数})$$取极限：
  令 $n\to \infty$，平行线间距趋于零，对于曲线上任意一点，都有$$\frac{\sin \alpha }{v}=C_1(\text{常数})$$其中 $\alpha$ 为该点切线与铅垂线夹角，如图所示.

$$\begin{array}{rl}& \because \alpha +\beta =\frac{\pi }{2},y^{\prime }=\tan \beta \\ & \therefore \sin \alpha =\cos \beta =\frac{1}{\sqrt{(\tan \beta {)}^2+1}}=\frac{1}{\sqrt{(y^{\prime }{)}^2+1}}\end{array}$$再由 $v=\sqrt{2gy}$，得$$\frac{\sin \alpha }{v}=\frac{1}{\sqrt{2gy}\sqrt{(y^{\prime }{)}^2+1}}=C_1\Rightarrow y(1+(y^{\prime }{)}^2)=C(\text{常数})$$

四、数学建模

  上面算是用物理方法得到的结论，下面让我们看看复杂而严谨的数学方法，首先需要一些数学基础，鉴于我也没学过，就抄一下.

1. 数学基础

(1) 泛函

  首先来介绍一种特殊的映射，泛函. 以前我们接触最多的映射是函数，即数集到数集的映射 $f:R\to R$. 所谓的泛函，指的是这样一种映射：

  设 $C$ 是一个函数集，对于 $C$ 中任何一个元素 $y(x)$，数集 $B$ 中都有唯一的一个元素 $J$ 与之对应，则称 $J$ 为 $y(x)$ 的泛函数，记作$$J=J[y(x)]$$  泛函的值是数，自变量是函数，泛函是函数的推广. 特别地，我们这里考虑的集合 $C$ 中的元素是这样的函数$$y(a)=y_1,y(b)=y_2$$其中 $y_1,y_2$ 为常数，$y(x)$ 充分光滑.

  几何上来看，$C$ 中的元素均是由点 $A(a,y_1)$ 到点 $B(b,y_2)$ 的光滑曲线.

  一般情况下，泛函可以用以下积分式来表示$$J[y(x)]={\int }_a^{b}F(x,y,y^{\prime })dx$$

(2) 泛函变分

  函数 $y=f(x)$ 的微分 $dy$ 指的是由于自变量 $x_0\to x_0+dx$ 所引起的因变量 $y$ 的变化量的线性近似 $$dy=f^{\prime }(x_0)dx=\frac{df}{dx}{|}_{x=x_0}dx$$  而泛函 $J[y(x)]$ 的变分 $\delta J$ 指的是由于自变量 $y\to y+\delta y$ 所引起的因变量 $J$ 的变化量的线性近似 $\delta J$.

  下面仿照函数微分来定义泛函变分.

  假定 $y(x)(\in C)$ 和 $F(x,y,y^{\prime })$ 充分光滑，考虑泛函$$J[y(x)]={\int }_a^{b}F(x,y,y^{\prime })dx$$当 $y$ 有微小变化 $\delta y$ ($\delta y$ 也是 $x$ 的函数) 时，函数 $F$ 的改变量约等于微分，即$$\Delta F=F(x,y+\delta y,y^{\prime }+\delta y^{\prime })-F(x,y,y^{\prime })\approx \frac{\partial F}{\partial y}\delta y+\frac{\partial F}{\partial y^{\prime }}\delta y^{\prime }=dF$$则其泛函的改变量约等于变分 $\delta J$ $$\begin{array}{rl}\Delta J& =J[y(x)+\delta y(x)]-J[y(x)]\\ & ={\int }_a^b[F(x,y+\delta y,y^{\prime }+\delta y^{\prime })-F(x,y,y^{\prime })]dx\\ & \approx {\int }_a^b(\frac{\partial F}{\partial y}\delta y+\frac{\partial F}{\partial y^{\prime }}\delta y^{\prime })dx\\ & =\delta J\end{array}$$

(3) 泛函极值

在一元函数 $y=y(x)$ 中，$y$ 在 $x=x_0$ 处取极值的必要条件是$$\frac{dy}{dx}{|}_{x=x_0}=0$$下面证明泛函 $J$ 取极值的必要条件是 $$\delta J=0$$或者$$\frac{\partial F}{\partial y}-\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)=0$$
  设泛函 $J$ 在 $y=y(x)$ 处取极值，取 $\delta y$ 的特殊形式：$$\delta y=ϵ\phi (x)$$其中 $ϵ$ 为任意小的实数，$\phi (x)$ 为充分光滑的任意函数，且满足$$\phi (a)=\phi (b)=0$$则函数 $y+\delta y\in C$. 令 $$\begin{array}{rl}G(ϵ)& =J[y(x)+ϵ\phi (x)]\\ & =J[y(x)]+\delta J\\ & =J[y(x)]+ϵ{\int }_a^b(\frac{\partial F}{\partial y}\phi (x)+\frac{\partial F}{\partial y^{\prime }}{\phi }^{\prime }(x))dx\end{array}$$(上式中第二个等号应该是约等号，但因为 $ϵ$ 是任意小的实数，同时为了证明能够顺利进行，暂且允许这么一点瑕疵)
  由 $J$ 在 $y=y(x)$ 处取极值可知，$G(ϵ)$ 在 $ϵ=0$ 处取极值，则有$$\frac{dG}{dϵ}{|}_{ϵ=0}={\int }_a^b(\frac{\partial F}{\partial y}\phi (x)+\frac{\partial F}{\partial y^{\prime }}{\phi }^{\prime }(x))dx=0$$两端同乘 $ϵ$，则有$$\delta J={\int }_a^b(\frac{\partial F}{\partial y}ϵ\phi (x)+\frac{\partial F}{\partial y^{\prime }}ϵ{\phi }^{\prime }(x))dx=0$$ 下证第二个必要条件：$$\frac{\partial F}{\partial y}-\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)=0$$
$$\begin{array}{rl}0& ={\int }_a^b(\frac{\partial F}{\partial y}\phi (x)+\frac{\partial F}{\partial y^{\prime }}{\phi }^{\prime }(x))dx\\ & ={\int }_a^b\frac{\partial F}{\partial y}\phi (x)dx+{\int }_a^b\frac{\partial F}{\partial y^{\prime }}d(\phi (x))\text{(应用分部积分)}\\ & ={\int }_a^b\frac{\partial F}{\partial y}\phi (x)dx+\frac{\partial F}{\partial y^{\prime }}\phi (x){|}_a^b-{\int }_a^b\phi (x)\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)dx(\phi (a)=\phi (b)=0)\\ & ={\int }_a^b[\frac{\partial F}{\partial y}-\frac{d}{dx}(\frac{\partial F}{\partial y^{\prime }}\left)\right]\phi (x)dx\end{array}$$
这里需要介绍一个引理：

  设 E = { g   ∣   g ( a ) = g ( b ) = 0   ,   g ∈ C 1 [ a , b ] } \small E=\{g\,|\,g(a)=g(b)=0\,,\,g\in C^1[a,b]\} E={g∣g(a)=g(b)=0,g∈C1[a,b]}，即由在 $[a,b]$ 上连续且端点函数值为零的函数组成的集合，若 $[a,b]$ 上的连续函数 $f$ 满足：

  $\forall g\in E$，都有$${\int }_a^{b}f(x)g(x)dx=0$$则 $f(x)\equiv 0$.

  套用该引理，由 $\phi (x)$ 的任意性，可得$$\frac{\partial F}{\partial y}-\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)=0$$这个式子被称为欧拉—拉格朗日方程，简记为 $\text{E}$-$\text{L}$ 方程.

  至此，泛函取极值的两个必要条件证毕！

2. 数学分析

下面对最速降线问题重新进行数学上的分析.

  设质点经过曲线上点 $P(x,y)$ 时速度为 $v$，由能量守恒定律$$mgy=\frac{1}{2}m{v}^2$$解得 $v=\sqrt{2gy}$. 弧微分 $ds=\sqrt{1+(y^{\prime }{)}^2}dx$，经过这一小段弧时近似认为速度不变，则所需时间为 $$dt=\frac{ds}{v}=\frac{\sqrt{1+(y^{\prime }{)}^2}}{\sqrt{2gy}}dx$$从 $A$ 下滑到 $B$ 所需总时间 $$T=\int dt={\int }_0^c\frac{\sqrt{1+(y^{\prime }{)}^2}}{\sqrt{2gy}}dx=\frac{1}{\sqrt{2g}}{\int }_0^c\frac{\sqrt{1+(y^{\prime }{)}^2}}{\sqrt{y}}dx$$令$$J[y(x)]={\int }_0^c\frac{\sqrt{1+(y^{\prime }{)}^2}}{\sqrt{y}}dx$$则 $J$ 是 $y(x)$ 的泛函，被积函数为 $F(y,y^{\prime })$ 的形式，显然也是 $F(x,y,y^{\prime })$ 的形式.

  我们的目标是求 $$y=\underset{C}{\mathrm{argmin}}T=\underset{C}{\mathrm{argmin}}J[y(x)]$$则 $J[y(x)]$ 在 $y$ 处取极值，由第二必要条件($\text{E}$-$\text{L}$ 方程)：$$\frac{\partial F}{\partial y}-\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)=0$$于是$$\begin{array}{rl}\frac{d}{dx}(F(y,y^{\prime })-y^{\prime }\frac{\partial F}{\partial y^{\prime }})& =\frac{\partial F}{\partial y}{y}^{\prime }+\frac{\partial F}{\partial y^{\prime }}{y}^{\prime \prime }-(y^{\prime \prime }\frac{\partial F}{\partial y^{\prime }}+y^{\prime }\frac{d}{dx}(\frac{\partial F}{\partial y^{\prime }}\left)\right)\\ & =y^{\prime }(\frac{\partial F}{\partial y}-\frac{d}{dx}(\frac{\partial F}{\partial y^{\prime }}\left)\right)\\ & =0\end{array}$$则$$F(y,y^{\prime })-y^{\prime }\frac{\partial F}{\partial y^{\prime }}=C_1(\text{常数})$$将 $$F(y,y^{\prime })=\frac{\sqrt{1+(y^{\prime }{)}^2}}{\sqrt{y}},\frac{\partial F}{\partial y^{\prime }}=\frac{1}{\sqrt{y}}\frac{y^{\prime }}{\sqrt{1+(y^{\prime }{)}^2}}$$ 代入，得$$\frac{\sqrt{1+(y^{\prime }{)}^2}}{\sqrt{y}}-y^{\prime }\frac{y^{\prime }}{\sqrt{y}\sqrt{1+(y^{\prime }{)}^2}}=\frac{1}{\sqrt{y}\sqrt{1+(y^{\prime }{)}^2}}=C_1(\text{常数})$$两边平方，则$$y(1+(y^{\prime }{)}^2)=C(\text{常数})$$

五、曲线求解

通过两种方法得到的结果相同，即最速曲线应满足$$y(1+(y^{\prime }{)}^2)=C(\text{常数})$$下面进行求解.

  根据上图 $y^{\prime }=\tan \beta =\cot \alpha$，因为质点一直向下滑(如果向上，肯定不是最优曲线)，所以 $0\le \alpha \le \pi /2$.

  令 $2R=C$，令 $y^{\prime }=dy/dx=\cot (\theta /2)$，则 $\theta =2\alpha$，所以 $0\le \theta \le \pi$.

  设 $x$ 是 $\theta$ 的函数，即 $x=x(\theta )$，我们的目标是求 $x(\theta )$.

由 $$1+{\cot }^2\frac{\theta }{2}=1+\frac{{\cos }^2\frac{\theta }{2}}{{\sin }^2\frac{\theta }{2}}=\frac{{\sin }^2\frac{\theta }{2}+{\cos }^2\frac{\theta }{2}}{{\sin }^2\frac{\theta }{2}}=\frac{1}{{\sin }^2\frac{\theta }{2}}=\frac{2}{1-\cos \theta }$$所以$$y(1+(y^{\prime }{)}^2)=y(1+{\cot }^2\frac{\theta }{2})=\frac{2y}{1-\cos \theta }=2R\Rightarrow y=R(1-\cos \theta )$$由$$\frac{dy}{dx}=\cot \frac{\theta }{2}=\frac{\cos \frac{\theta }{2}}{\sin \frac{\theta }{2}},\frac{dy}{d\theta }=R\sin \theta =2R\sin \frac{\theta }{2}\cos \frac{\theta }{2}$$则$$x^{\prime }(\theta )=\frac{dx}{d\theta }=\frac{dy}{d\theta }/\frac{dy}{dx}=2R\sin \frac{\theta }{2}\cos \frac{\theta }{2}\frac{\sin \frac{\theta }{2}}{\cos \frac{\theta }{2}}=2R{\sin }^2\frac{\theta }{2}=R(1-\cos \theta )$$所以$$x(\theta )=\int x^{\prime }(\theta )d\theta =\int R(1-\cos \theta )d\theta =R(\theta -\sin \theta )+C_1$$

得到 $x,y$ 关于 $\theta$ 的参数方程$$\{\begin{array}{l}x(\theta )=R(\theta -\sin \theta )+C_1\\ y(\theta )=R(1-\cos \theta )\end{array},0\le \theta \le \pi$$还有两个待定系数 $R,C_1$，别忘了我们还有两个边界条件没用呢，即$$y(0)=0,y(c)=H$$先消去参数 $\theta$，作变量代换 $$\theta =\arccos (1-\frac{y}{R}),0\le \theta \le \pi$$则$$\begin{array}{rl}x& =R\arccos (1-\frac{y}{R})-R\sqrt{2\frac{y}{R}-\frac{y^2}{R^2}}+C_1\\ & =R\arccos (1-\frac{y}{R})-\sqrt{y(2R-y)}+C_1\end{array}$$将边界条件换种描述方式，即$$x(0)=0,x(H)=c$$将其代入方程，得$$\begin{array}{rl}x(0)& =R\arccos (1)+C_1=C_1=0\\ x(H)& =R\arccos (1-\frac{H}{R})-\sqrt{H(2R-H)}+0=c\end{array}$$则 $C_1=0,R$ 的值可以由 $H,c$ 确定. 参数方程变为$$\{\begin{array}{l}x(\theta )=R(\theta -\sin \theta )\\ y(\theta )=R(1-\cos \theta )\end{array},0\le \theta \le \pi$$其实 $\theta$ 的值只有在 $c/H=\pi /2$ 时才能取到 $\pi$ (看完下一部分你便知晓原因了)，其余情况下 $\theta <\pi$.

  可能有些读者已经看出来了，上述参数方程正是滚轮线(也叫摆线)，参数 $R$ 表示滚轮半径. 所以最速降线是滚轮线的一部分.

六、结果解释

滚轮线的参数方程推导过程如下：

  由于轨迹是无摩擦匀速滚动产生的，所以轮与地面接触点 $B$ 的速度始终为 $0$，$B$ 点为瞬时旋转中心.

  设转过的角度为 $\theta$，角速度为 $w$，则圆心 $A$ 的速度 $v=wR$，移动的距离 $EA=vt=Rwt=R\theta$. 设点 $C$ 坐标 $(x(\theta ),y(\theta ))$，则 $EC=x(\theta ),BD=y(\theta )$.

$$\begin{array}{rl}EA=EC+CA& \Rightarrow R\theta =x(\theta )+R\sin \theta \\ BA=BD+DA& \Rightarrow R=y(\theta )+R\cos \theta \end{array}$$所以滚轮线的参数方程为$$\{\begin{array}{l}x(\theta )=R(\theta -\sin \theta )\\ y(\theta )=R(1-\cos \theta )\end{array},0\le \theta \le 2\pi$$参数 $\theta$ 表示转过的角度.

  下面来验证滚轮线就是最速降线，根据物理建模的结论(即下面这个图，其中 $\alpha$ 为该点切线与铅垂线夹角)：

最速降线需满足 $$\frac{\sin \alpha }{v}=C(常数)$$ 由 $v=\sqrt{2gy}$，得 $$\frac{\sin \alpha }{\sqrt{y}}=C(常数)\Rightarrow \frac{{\sin }^2\alpha }{y}=C(常数)$$
对滚轮线的图作些变动

  因为 $BF$ 为直径，所以 $BC\perp CF$.
  由 $B$ 点瞬时速度为 $0$，所以 $C$ 点的速度方向与 $BC$ 垂直.
  即 $\overrightarrow{CF}$ 为速度方向，直线 $CF$ 为滚轮线的切线.

  则两张图中的 $\alpha$ 是一致的.
  因为 $CG//BF$，所以 $\angle AFC=\angle FCG=\alpha$.
  由弦切角定理，$\angle GBC=\angle AFC=\alpha$. 所以$$y=CG=BC\sin \alpha =BF\sin \alpha \sin \alpha =BF{\sin }^2\alpha =2R{\sin }^2\alpha$$所以 $$\frac{{\sin }^2\alpha }{y}=\frac{1}{2R}(\text{常数})$$验证完毕.

最速降线是滚轮线这一巧合已经足以让人称奇，下面来看一件更加巧合的事情.

  考虑轮子上点 $C$ 的速度，设此时 $C$ 的纵坐标也是 $y$，对应角度为 $\alpha$， 根据上述过程，有 $y=2R{\sin }^2\alpha$. 因为点 $B$ 是瞬时旋转中心，所以 $C$ 的速度 $v_c=w\cdot BC=w\cdot 2R\sin \alpha$.

  我们已经知道，如果质点沿最速降线下滑，则其在 $(x,y)$ 处的速度为 $v=\sqrt{2gy}$. 将 $y=2R{\sin }^2\alpha$ 代入，得 $v=\sqrt{2g\cdot 2R{\sin }^2\alpha }=2\sqrt{gR}\sin \alpha$.

  对比这两个式子 $$v_c=w\cdot 2R\sin \alpha$$$$v=2\sqrt{gR}\sin \alpha$$若 $w=\sqrt{\frac{g}{R}}$，则 $v=v_c$. 这说明了什么？我想大家应该已经清楚了. 即若轮子旋转的角速度与半径之间满足 $w=\sqrt{g/R}$，不仅得到的曲线是最速降线，且轮上点的运动速度与沿这条曲线运动(可能是下滑也可能是上滑)的质点的速度相同，$\text{How amazing it is!}$

(来源：下方B站视频)

七、补充资料

B站视频：最速降线问题——最快下滑路径为什么是旋轮线？（中英字幕）

维基百科：
1. Brachistochrone_curve(最速降线) (里面有"神仙"们的解法)
2. 旋轮线 (里面有其他有趣的性质)

参考文献：
1.于涛.数学物理方程与特殊函数[M].北京:科学出版社,2008.

Plus: 如有错误、可以改进的地方、或任何想说的，请在评论区留言！