摆线知识

圆形是我们学习数学以来最先接触的特殊图形——它没法被分解成有限个三角形。我们如果在一个滚轮的边上一定点放置一支笔,然后让滚轮沿着墙壁滚动,就能够在墙上画出摆线(Cycloid)的图像:

不妨再用几何的方法来得到其参数方程:

经过一定的分析,我们最终可以得到一组方程来描述摆线(其中a为圆的半径):

\left\{ \begin{aligned} x&=a(\theta-\sin\theta) \\ y&=a(1-\cos\theta) \end{aligned} \right.

一、摆线的基本性质

根据三角函数的性质以及摆线的绘制过程,我们知道0\le\theta\le2\pi内(即滚轮旋转一周)的曲线构成摆线的一个弧。所以我们可以根据这个性质来求出摆线一个周期的弧长。

第一步:计算弧微分

\begin{aligned} & {\mathrm{d}x\over\mathrm{d}\theta}=a(1-\cos\theta) \\ & {\mathrm{d}y\over\mathrm{d}\theta}=a\sin\theta \end{aligned}

根据勾股定理(\mathrm{d}s)^2=(\mathrm{d}x)^2+(\mathrm{d}y)^2,我们得到:

\begin{aligned} {\mathrm{d}s\over\mathrm{d}\theta} &=\sqrt{a^2(1-\cos\theta)^2+a^2\sin^2\theta} \\ &=a\sqrt{(1-\cos\theta)^2+\sin^2\theta} \\ &=a\sqrt{1-2\cos\theta+\cos^2\theta+\sin^2\theta} \\ &=a\sqrt{1-2\cos\theta+1} \\ &=a\sqrt{2-2(1-2\sin^2(\theta/2))} \\ &=a\sqrt{4\sin^2(\theta/2)} \\ &=2a\sin\left(\theta\over2\right) \end{aligned}

由于我们\theta的取值范围内正弦函数不会取负值,因此我们可以放心的去掉绝对值符号。

第二步:计算曲线积分

\int_C\mathrm{d}s=\int_0^{2\pi}2a\sin\left(\theta\over2\right)\mathrm{d}\theta=\left.-4a\cos\left(\theta\over2\right)\right|^{2\pi}_0=8a

因此,一个滚轮转动一周形成摆线的弧长刚好是滚轮半径的八倍。

二、最速降线问题

最速降线问题,即找到两个定点间连接的曲线,使得物体能够在重力作用下花最短的时间到达终点。早在1696年,著名数学家约翰·伯努利(Johann Bernoulli)就提出了这个问题。当时欧洲各地的数学家都使用了不同的方法来求解这个问题。然而,伯努利本人提出这个挑战的初衷是为了凸显自己方法的高明从而体现他自己很聪明。他本人花了近两个星期找到这个方法。然而他向欧洲发表这个问题后仅仅一天,已经退休牛顿就给出了他的解法。牛顿在后来的书信里曾说道:"我不喜欢让外国人嘲弄我的数学能力"。后来,伯努利的学生欧拉(Leonard Euler)不仅给出了最速降线问题的另一种解法,还设计了欧拉—拉格朗日方程(Euler—Lagrange equation)来解决类似的问题。本节中我将着重讲述伯努利和欧拉的解法。

在直接进入最速降线问题的推导前,我们需要提及费马原理(Fermat's principle),即在两点间,光总是选择耗时最短的路径。利用这个原理,我们可以研究一下光的折射性质:

现在我们想研究一下入射角与折射角的关系。不同于求证光的反射定律,我们需要用上述的费马原理来求解\alpha\beta的关系。假设光在白色区域内的传播速度为v_1,而它在蓝色区域的传播速度为v_2,现在我们要找到合适的x,使得光沿着红色路线的传播时间最少(即找到最快传播路径)。

本篇文章将使用微积分里优化函数的方法来求解。根据勾股定理,我们不难得出:

l_1=\sqrt{x^2+a^2} \\ l_2=\sqrt{(d-x)^2+b^2}

现在,我们再列出要优化的时间函数:

T={l_1\over v_1}+{l_2\over v_2}={\sqrt{x^2+a^2}\over v_1}+{\sqrt{(x-d)^2+b^2}\over v_2}

现在对x求导,得到:

{\mathrm{d}T\over\mathrm{d}x}={x\over v_1\sqrt{x^2+a^2}}-{d-x\over v_2\sqrt{(d-x)^2+b^2}}

我们知道,当导数为0或未定义时,函数取极值。然而在在T函数的定义域内,函数只在导数为0时达到最小值。因此,我们能够列出这样的方程来求解x:

& {x\over v_1\sqrt{x^2+a^2}}-{d-x\over v_2\sqrt{(d-x)^2+b^2}}=0 \\ & {x\over v_1\sqrt{x^2+a^2}}={d-x\over v_2\sqrt{(d-x)^2+b^2}}

根据原图,我们不难得出\sin\alpha={x\over\sqrt{x^2+a^2}},\cos\alpha={d-x\over\sqrt{(d-x)^2+b^2}}。因此,代入到原式,我们就得到了入射角\alpha与折射角\beta的关系:

{\sin\alpha\over v_1}={\sin\beta\over v_2} \Rightarrow{\sin\alpha\over\sin\beta}={v_1\over v_2}

因此,入射角与折射角的正弦比刚好等于两种介质下的光速比,这个结论又叫做斯涅尔定律(Snell's law)。虽然斯涅尔定律是一个光学定律,但是在力学中它依然适用。

三、伯努利解法——用斯涅尔定律求解最速降线

我们设物体最初距离地面的高度为h,而物体的当前速度和高度分别为v=v(t)y=y(t)。根据(理想状态下)机械能守恒,我们可以对物体的速度和当前高度建立如下关系:

{1\over2}mv^2+mgy=mgh

由于物体的质量恒不为0,所以我们可以直接让等式两侧与之相除:

{1\over2}v^2+gy=gh

根据基本代数操作,我们不难得出物体速度关于其纵坐标的表达式:

v=\sqrt{2g(h-y)}

根据幂函数的性质,物体的纵坐标其速度一一对应。如果我们把每个纵坐标都当作一个不同的介质而每个不同的介质下的光速都为\sqrt{2g(h-y)},我们可以把最速降线问题理解成光在多个不同介质下的传播问题。根据费马原理,光总是沿着耗时最短的路径传播,所以我们如果用光学的原则求解出耗时最短的传播路径,就可以得到最速降线问题的答案。

如果我们把物体的运动空间用平行于x轴的直线进行分割,就可以发现当"光线"从第j个区间穿越到第j+1个区间进行折射时,折射角满足斯涅尔定律。因此:

{\sin\alpha_j\over v_j}={\sin\alpha'_j\over v_{j+1}}

由于我们所做的发现均垂直于x轴(即平行于y轴),所以这些法线都相互平行。于是,根据两直线平行,内错角相等的原则,我们得到\alpha'_j=\alpha_{j+1},因此,我们可以把上述等式进行这样的拓展:

{\sin\alpha_1\over v_1}={\sin\alpha_2\over v_2}=\dots={\sin\alpha_j\over v_j}={\sin\alpha_{j+1}\over v_{j+1}}=\dots={\sin\alpha_N\over v_N}=C

因此每一个区间内入射角的正弦于区间内的传播速度成正比,而比例常数设为C。因此,我们可以再把速度关于y的表达式代入,得到:

{\sin\alpha\over\sqrt{2g(h-y)}}=C

经过一定代数操作,我们可以得到y关于\sin\alpha的表达式:

\begin{aligned} {\sin^2\alpha\over2g(h-y)}&=C^2 \\ {\sin^2\alpha\over2gC^2}&=h-y \\ y&=-{\sin^2\alpha\over2gC^2}+h \end{aligned}

现在我们再来看一看入射角于物体坐标微小变化的关系:

我们不难发现{\mathrm{d}x\over\mathrm{d}y}=\tan\alpha,因此我们可以通过\alpha构造一个最速降线的参数方程。现在我们先让y\alpha求导,得到下式:

{\mathrm{d}y\over\mathrm{d}\alpha}=-{\sin\alpha\cos\alpha\over gC^2}

根据链式法则,我们可以得到x关于\alpha的微分方程:

{\mathrm{d}x\over\mathrm{d}\alpha}={\mathrm{d}x\over\mathrm{d}y}{\mathrm{d}y\over\mathrm{d}\alpha}=\tan\alpha\cdot\left(-{\sin\alpha\cos\alpha\over gC^2}\right)=-{\sin^2\alpha\over gC^2}

现在两侧积分,我们可以得到:

\begin{aligned} x&=-{1\over gC^2}\int\sin^2\alpha\mathrm{d}\alpha \\ &=-{1\over4gC^2}\int(2-2\cos2\alpha)\mathrm{d}\alpha \\ &=-{1\over4gC^2}(2\alpha-\sin2\alpha)+x_0 \end{aligned}

现在我们可以将y的表达式改写,得到:

y=-{1\over4gC^2}(1-\cos2\alpha)+h

为了美观,我们可以更换自变量\theta=2\alpha,得到最速降线最后的参数方程:

\left\{ \begin{aligned} x&=-{1\over4gC^2}(\theta-\sin\theta)+x_0 \\ y&=-{1\over4gC^2}(1-\cos\theta)+h \end{aligned} \right.

我们不难发现最速降线的参数方程画出的图像刚刚好就是翻转后的摆线:

因此,最速降线的答案就是摆线。以上便是伯努利的解法,现在我们再来看看他的学生——欧拉——的解法。

四、欧拉的解法——变分法

我们先来看一下什么是泛函(functional)。泛函一般指从一个空间到数域的映射。

对于一个实变函数f,f\mapsto f(x_0)就是一个基本的泛函,它把函数空间映射到了实数域。而另外一种常见泛函则由积分构成:

J[y]=\int_{x_0}^{x_1}F(x,y,y')\mathrm{d}x

或许读者看到这里会有些摸不着头脑,但是我们这些内容会为后面的求解提供铺垫。现在我们来看一看泛函极值问题:找到一个满足y(x_0)=y_0,y(x_1)=y_1y,使得泛函J[y]=\int_{x_0}^{x_1}F(x,y,y')\mathrm{d}x取得极值。

在传统的单变量微积分里的优化问题中,我们一般都会通过寻找导数的零点来寻找函数的极值。在泛函里我们会寻找变分(variation)的零点来求得泛函的极值。如果我们假设y就是最终满足泛函达到极值的函数,而u=y+\eta表示的是y领域内的函数(其中\eta(x_0)=\eta(x_1)=0以便于所有的u满足y的边界条件)。则:

\delta J\equiv J[y+\eta]-J[y]=\int_{x_0}^{x_1}\left[F(x,y+\eta,y'+\eta')-F(x,y,y')\right]\mathrm{d}x

对于足够"微小"的y,我们可以用泛函核FF(x,y,y')附近的一阶泰勒展开来"估计"F(x,y+\eta,y'+\eta):

\begin{aligned} \delta J&=\int_{x_0}^{x_1}\left[F(x,y,y')+\nabla F\cdot(0,\eta,\eta')-F(x,y,y')\right]\mathrm{d}x \\ &=\int_{x_0}^{x_1}\nabla F\cdot(0,\eta,\eta')\mathrm{d}x \\ &=\int_{x_0}^{x_1}\left[{\partial F\over\partial y}\eta+{\partial F\over\partial y'}\eta'\right]\mathrm{d}x \end{aligned}

现在我们对被积函数里的第二项进行分部积分,得:

\delta J=\int_{x_0}^{x_1}{\partial F\over\partial y}\eta\mathrm{d}x+\left[{\partial F\over\partial y'}\eta\right]^{x_1}_{x_0}-\int_{x_0}^{x_1}\eta{\mathrm{d}\over\mathrm{d}x}{\partial F\over\partial y}\mathrm{d}x

根据\eta满足\eta(x_0)=\eta(x_1)=0的性质,我们可以得到泛函J变分的最简形式:

\delta J=\int_{x_0}^{x_1}\eta\left[{\partial F\over\partial y}-{\mathrm{d}\over\mathrm{d}x}{\partial F\over\partial y'}\right]\mathrm{d}x

由于在J[y]处取得极值,所以根据"直觉",我们可以得到\delta J=0。我们也不难能看出(变分法基本引理),\delta J=0当且仅当被积函数里\eta的系数在区间内恒为0。因此,我们不严谨地推导出了:

{\partial F\over\partial y}-{\mathrm{d}\over\mathrm{d}x}{\partial F\over\partial y'}=0

即大名鼎鼎的欧拉—拉格朗日方程(Euler—Lagrange equation),简称E-L方程。知乎有很多其它的推导过程,但是它们大多数都设u=y+\varepsilon\eta然后求泛函对\varepsilon的全微分。虽然这样更加的严谨,但是作者个人最容易理解的推导还是以上u=y+\eta的方法。

虽然E-L方程可以解决大多数泛函极值问题,但是对于满足{\partial F\over\partial x}=0的泛函,直接代入E-L方程会比较产生比较复杂的微分方程。因此,我们决定再简化一下难度,对E-L方程进行移项:

{\partial F\over\partial y}={\mathrm{d}\over\mathrm{d}x}{\partial F\over\partial y'}

对两侧乘上y'并积分,得到:

\begin{aligned} \int{\partial F\over\partial y}y'\mathrm{d}x&=\int y'{\mathrm{d}\over\mathrm{d}x}{\partial F\over\partial y'}\mathrm{d}x \\ \int{\partial F\over\partial y}y'\mathrm{d}x&=y'{\partial F\over\partial y'}-\int{\partial F\over\partial y'}y''\mathrm{d}x \\ \int\left[{\partial F\over\partial y}y'+{\partial F\over\partial y'}y''\right]\mathrm{d}x&=y'{\partial F\over\partial y'}+C \\ \int{\mathrm{d}F\over\mathrm{d}x}\mathrm{d}x&=y'{\partial F\over\partial y'}+C \\ F-y'{\partial F\over\partial y}&=C \\ \end{aligned}

因此,我们可以在优化满足F=F(y,y')的泛函使用我们得到的结论,即贝尔特拉米等式(Beltrami's identity)。我们已经将所有需要的工具准备完毕,现在可以开始进入最速降线问题的求解了。

用变分法求解最速降线问题

由于我们要找到一个曲线,使得物体沿曲线路径滑落的时间最短,我们需要设计一个描述时间的泛函:

T=\int_C{\mathrm{d}s\over v}

利用弧微分公式和前面提到的能量守恒,我们不难得出:

{\mathrm{d}s}=\sqrt{1+y'^2}\mathrm{d}x \\ v=\sqrt{2g(h-y)}

因此,我们的得到了这样的时间泛函:

T=\int_{x_0}^{x_1}\underbrace{\sqrt{1+y'^2}\over\sqrt{2g(h-y)}}_F\mathrm{d}x

因此,我们的最速降线问题变成变分问题,即找到一个满足边界条件的y,使得泛函T取得极值。由于泛函T的核F满足{\partial F\over\partial x}=0,所以我们可以使用更加方便的贝尔特拉米等式来求解。首先,我们求F关于y'的偏导数:

{\partial F\over\partial y'}={y'\over\sqrt{2g(h-y)}\sqrt{1+y'^2}}

代入到贝尔特拉米等式,得到:

\begin{aligned} F-y'{\partial F\over\partial y'}&=C \\ {\sqrt{1+y'^2}\over\sqrt{2g(h-y)}}-{y'^2\over\sqrt{2g(h-y}\sqrt{1+y'^2}}&=C \\ {1+y'^2-y'^2\over\sqrt{2g(h-y)}\sqrt{1+y'^2}}&=C \\ {1\over\sqrt{2g(h-y)}\sqrt{1+y'^2}}&=C \end{aligned}

此时,我们如果令y'\equiv{\mathrm{d}y\over\mathrm{d}x}=\cot\alpha,则根据三角恒等式\csc^2\alpha=1+\cot^2\alpha,原式变成:

{\sin\alpha\over\sqrt{2g(h-y)}}=C

我们可以发现这个方程和我们用斯涅尔定律得出的方程完全一样,所以其实伯努利和欧拉的方法都能够说明最速降线问题的答案是摆线。然而,摆线的故事才刚刚开始。摆线不仅是最速降线问题的答案,还是等时降落问题(Tautochrone problem)的解。

摆线与等时降落问题——弹簧与硬核的阿贝尔

等时降落问题,即找到一个合适的曲线,使得物体从曲线的任何位置受重力作用(理想状态下无摩擦)滑落到终点的时间不变。

当然,在我们有头绪之前。我们先来研究研究弹簧振子的简单调谐运动(Simple harmonic motion)。其中弹簧振子的位置x满足如下微分方程:

{\mathrm{d}^2x\over\mathrm{d}t^2}+\omega^2x=0

我们设初始条件:x(0)=x_0,\left.\mathrm{d}x\over\mathrm{d}t\right|_{t=0}=0,即物体的初速度为0而初始位置为x_0。因此,在我们对微分方程进行拉氏变换时能够得到:

\begin{aligned} s^2\mathcal{L}\{x\}-sx_0+\omega^2\mathcal{L}\{x\}&=0 \\ \mathcal{L}\{x\}&={sx_0\over s^2+\omega^2} \\ \mathcal{L}\{x\}&=x_0\mathcal{L}\{\cos(\omega t)\} \\ x&=x_0\cos(\omega t) \end{aligned}

根据余弦函数的性质,我们不难发现x的第一个零点在t={\pi\over2\omega},而且这个零点与x_0的值无关。这意味着(理想状态下)无论x_0多大,{\pi\over2\omega}秒后x总能回到0。意思就是说无论你把弹簧拉到多长,弹簧总能够在指定时间后回到原点。

利用简谐运动求解等时降落问题

我们引入广义坐标s来描述示意图终点与物体当前位置间曲线的弧长,则根据等时降落曲线满足:s变成0的耗时与s_0无关,这意味着s可能会满足简谐运动方程,即:

{\mathrm{d}^2s\over\mathrm{d}t^2}=-\omega^2s

现在我们在来看看作用在物体上的加速度。

由于蓝色部分的分量刚好被约束抵消,所以物体在曲线上的加速度刚好等于-g\sin\theta,因此我们不难发现{\mathrm{d}^2s\over\mathrm{d}t^2}=-g\sin\theta。再代入到s满足的简谐运动方程,我们得到了等时降落曲线的微分方程:

-g\sin\theta=-\omega^2s

现在我们对等式两侧两侧微分,则

g\cos\theta\mathrm{d}\theta=\omega^2\mathrm{d}s \Rightarrow\mathrm{d}s={g\cos\theta\over\omega^2}\mathrm{d}\theta

读者可以根据上图\theta以及三角函数本身的定义推出一下面两个式子:

\begin{aligned} \mathrm{d}x&=\cos\theta\mathrm{d}s \\ \mathrm{d}y&=\sin\theta\mathrm{d}s \end{aligned}

分别将以上两个等式代入到原来的微分方程中,我们可以得到两个方程:

\left\{ \begin{aligned} \mathrm{d}x&={g\over\omega^2}\cos^2\theta\mathrm{d}\theta \\ \mathrm{d}y&={g\over\omega^2}\sin\theta\cos\theta\mathrm{d}\theta \end{aligned} \right.

现在对两侧分别进行积分,得到:

\begin{aligned} x&={g\over\omega^2}\int\cos^2\theta\mathrm{d}\theta \\ &={g\over4\omega^2}\int(2+2\cos2\theta)\mathrm{d}\theta \\ &={g\over4\omega^2}(2\theta+\sin2\theta)+x_0 \end{aligned}

以及

\begin{aligned} y&={g\over\omega^2}\int\sin\theta\cos\theta\mathrm{d}\theta \\ &={g\over2\omega^2}\sin^2\theta+y_0 \\ &={g\over4\omega^2}(1-\cos2\theta)+y_0 \end{aligned}

依照上面最速降线的传统,我们可以进行换元:\phi=2\theta,则我们得到的等时降落曲线方程为:

\left\{ \begin{aligned} x&={1\over4g\omega^2}(\phi+\sin\phi)+x_0 \\ y&={1\over4g\omega^2}(1-\cos\phi)+y_0 \end{aligned} \right.

因此,我们利用简谐运动方程发现(翻转后的)摆线是等时降落问题的答案。这就是等时降落问题的一个传统解法,但是阿贝尔用却能量守恒给出了一个更加硬核的解答。

阿贝尔的硬核解法——积分方程法

与简谐运动方程法一样,阿贝尔也采用了同样的广义坐标s。现在把s代入到机械能守恒方程里(h为物体初始位置的纵坐标),得到(符号是因为物体下滑时s单调递减):

{1\over2}m\left(\mathrm{d}s\over\mathrm{d}t\right)^2+mgy=mgh \Rightarrow{\mathrm{d}s\over\mathrm{d}t}=-\sqrt{2g(h-y)}

根据问题的定义,我们知道等时降落问题的曲线满足物体降落时间与出发位置无关,即曲线满足\int_{y=h}^{y=0}\mathrm{d}t=T(T为常数)。利用链式法则,我们发现\mathrm{d}t={\mathrm{d}t\over\mathrm{d}s}{\mathrm{d}s\over\mathrm{d}y}\mathrm{d}y=-{1\over\sqrt{2g}\sqrt{h-y}}{\mathrm{d}s\over\mathrm{d}y}\mathrm{d}y,所以最后代入到原来的积分式,得到:

\sqrt{2g}T=\int_0^h(h-y)^{-1/2}{\mathrm{d}s\over\mathrm{d}y}\mathrm{d}y

于是我们得到了著名的阿贝尔积分方程(Abel integral equation)。在上一篇文章里,我使用分数阶微积分的方法求解了这个方程,但这次我打算使用拉普拉斯变换来求解:

\mathcal{L}\{\sqrt{2g}T\}=\mathcal{L}\left\{\int_0^h(h-y)^{-1/2}{\mathrm{d}s\over\mathrm{d}y}\mathrm{d}y\right\}

我们不难发现等式右侧刚好是一个卷积,因此根据卷积定理(Convolution theorem),右侧可以改写成两个拉氏变换的乘积:

\mathcal{L}\{\sqrt{2g}T\}=\mathcal{L}\{y^{-1/2}\}\mathcal{L}\left\{\mathrm{d}s\over\mathrm{d}y\right\} [公式]

再利用前期文章中总结出来的拉氏变换性质\mathcal{L}\{t^a\}={\Gamma(a+1)\over s^{a+1}},我们得到:

\begin{aligned} {\sqrt{2g}T\over s}&={\Gamma(1/2)\over s^{1/2}}\mathcal{L}\left\{\mathrm{d}s\over\mathrm{d}y\right\} \\ \mathcal{L}\left\{\mathrm{d}s\over\mathrm{d}y\right\}&={\sqrt{2g}T\over s^{1/2}}{1\over\Gamma(1/2)} \\ \mathcal{L}\left\{\mathrm{d}s\over\mathrm{d}y\right\}&={\sqrt{2g}T\over\Gamma^2(1/2)}\mathcal{L}\{y^{-1/2}\} \\ {\mathrm{d}s\over\mathrm{d}y}&={\sqrt{2g}T\over\pi\sqrt{y}} \end{aligned}

对等式两侧平方,得到\left(\mathrm{d}s\over\mathrm{d}y\right)^2={2gT^2\over\pi^2y}。根据弧微分公式(\mathrm{d}s)^2=(\mathrm{d}x)^2+(\mathrm{d}y)^2,倘若我们设{\mathrm{d}x\over\mathrm{d}y}\cot\theta,则\left(\mathrm{d}s\over\mathrm{d}y\right)^2=\csc^2\theta。这意味着以下方程是成立的:

\csc^2\theta={2gT^2\over\pi^2y} \Rightarrow y={2gT^2\over\pi^2}\sin^2\theta={gT^2\over\pi^2}(1-\cos2\theta)

现在我们以及得到等时降落曲线方程的一半了,另一半可以通过对{\mathrm{d}x\over\mathrm{d}\theta}的积分求得:

{\mathrm{d}x\over\mathrm{d}\theta}={\mathrm{d}x\over\mathrm{d}y}{\mathrm{d}y\over\mathrm{d}\theta}=\cot\theta{4gT^2\over\pi^2}\sin\theta\cos\theta={4gT^2\over\pi^2}\cos^2\theta
x={4gT^2\over\pi^2}\int\cos^2\theta\mathrm{d}\theta={gT^2\over\pi^2}(2\theta+\sin2\theta)+x_0

依照传统,我们通过\phi=2\theta换元能得到等时降落曲线简洁的参数方程:

\left\{ \begin{aligned} x&={gT^2\over\pi^2}(\phi+\sin\phi)+x_0 \\ y&={gT^2\over\pi^2}(1-\cos\phi) \end{aligned} \right.

我们不难发现,两种方法都可以说明等时降落问题的答案,如同最速降线问题,都是摆线。

参考地址:https://zhuanlan.zhihu.com/p/126421949