最速降线与悬链线的大致推导过程

作者：

Ch4ot1c_M@dn3ss , 2023-05-27 00:11:02 , 所有人可见 , 阅读 71

4

3

上高中以后一直好奇最速下降曲线是怎么算出来的。到了大学怎么还能留着问题呢？今天就~~用漏洞百出的证法~~把它扬了。

前提

假定无摩擦，起点为 $(0,0)$，终点是 $(x_0, y_0)$，整套系统只有重力场和最速降线，并建立 $y$ 轴竖直向下，$x$ 轴水平向右的左手坐标系。

分析过程

一

对质点而言，由经典力学动能定理有，$\dfrac{1}{2}mv^2=mgy$。消去质量，并变形得到 $v^2=2gy$。

同时又由于质点在最速降线上运动，此时也有 $\dfrac{\sqrt{(\mathrm{d}x)^2 + (\mathrm{d}y)^2}}{\mathrm{d} t} = v $。进一步变形得到 $\dfrac{ \sqrt{1 + (y’)^2}\mathrm{d}x}{\mathrm{d} t} = v $。

那么对于时间的微分而言，有 $\mathrm{d} t=\sqrt{\dfrac{ 1 + (y’)^2}{2gy}}\mathrm{d}x$。

积分回来就有

$$t=\int_0^{x_0}\sqrt{\dfrac{ 1 + (y’)^2}{2gy}}\mathrm{d}x$$

问题便是选取合适的 $y$ 使得 $t$ 在积分之后取到最小值。

不难发现此时的函数涉及到了三个变量：自变量 $x$，因变量 $y$ 以及因变量的导数 $y’$。

除此之外还能注意到 $y(0) = 0, y(x_0) = y_0$。然后就似乎没有额外的限制条件了。

二

我们从物理常识可以知道待求的曲线一定存在、存在极值、连续、可导且二阶可导。

而为了研究如何找到这样的函数，我们可以很自然地类比微分寻找整个函数的极值点一样 引入函数的变动量（或者称之为函数的变分） $\delta(x)$，在两点间 所有可能 的函数为 $\bar{y}$，最优解为 $y$。

那么，$\bar{y}=y+\delta$。从而有 $\bar{y}’=y’+\delta’$。相应地，要求 $\delta(0)=\delta(x_0) = 0$ 以固定起始点和终点。

而为了使用到极限和微分求极值的理论，我们还需要对能式子中任意性最大的 $\delta(x)$ 做变换 $\delta(x) = \xi\eta (x)$。其中 $\xi$ 是常数。

具体而言，我们考虑三个变量所构成的多元函数 $\bar{F}(x,y+\xi\eta ,y’ + \xi\eta’ )=\bar{F}(x,\bar{y},\bar{y}’)$ 以及形成最优解函数的多元函数 $F(x,y,y’)$，其中最优解函数存在于 $\lim\limits_{\xi \to 0}\bar{F}(x,\bar{y},\bar{y}’)$ 所确定的函数族之中。

同时，$\bar{F}(x,\bar{y},\bar{y}’)$ 也是一个 关于 $\xi$ 的函数。

$$I(\xi) = t = \int_0^{x_0}\bar{F}(x,\bar{y},\bar{y}’)\mathrm d x$$

最优解可通过 $\xi$ 逐渐趋近以至于成为零得到，且这个过程隐含一个条件：当 $\xi$ 趋近于零时，$I(\xi)$ 整体对 $\xi$ 的微分一定是零。这样才能保证 $I(\xi)$ 取到极值。

那么整个过程可以被表示为
$$\lim_{\xi \to 0}\dfrac{\mathrm d I(\xi)}{\mathrm d \xi} = \lim_{\xi \to 0}\Big(\dfrac{\mathrm d }{\mathrm d \xi}\int_0^{x_0}\bar{F}(x,\bar{y},\bar{y}’)\mathrm d x\Big) = \lim_{\xi \to 0}\int_0^{x_0}\dfrac{\mathrm d \bar{F}}{\mathrm d \xi}\mathrm d x = 0$$

根据函数性质，以及链式法则，可以得到

$$\dfrac{\mathrm d \bar{F}}{\mathrm d \xi} = \dfrac{\partial \bar{F}}{\partial x}\cdot\dfrac{\partial x}{\partial \xi} + \dfrac{\partial \bar{F}}{\partial \bar{y}}\cdot\dfrac{\partial \bar{y}}{\partial \xi} + \dfrac{\partial \bar{F}}{\partial \bar{y}’}\cdot\dfrac{\partial \bar{y}’}{\partial \xi} = 0 + \eta\dfrac{\partial \bar{F}}{\partial \bar{y}} + \eta’\dfrac{\partial \bar{F}}{\partial \bar{y}’}$$

进一步有

$$\begin{aligned}\lim_{\xi \to 0}\dfrac{\mathrm d I(\xi)}{\mathrm d \xi} &= \lim_{\xi \to 0}\int_0^{x_0} \eta\dfrac{\partial \bar{F}}{\partial \bar{y}} + \eta’\dfrac{\partial \bar{F}}{\partial \bar{y}’}\mathrm d x \\
&= \int_0^{x_0} \lim_{\xi \to 0}\Big(\eta\dfrac{\partial \bar{F}}{\partial (y + \xi\eta)} + \eta’\dfrac{\partial \bar{F}}{\partial (y’ + \xi\eta’)}\Big)\mathrm d x \\ \\
&= \int_0^{x_0}\Big( \eta\dfrac{\partial F}{\partial y} + \eta’\dfrac{\partial F}{\partial y’}\Big)\mathrm d x = \int_0^{x_0}\eta\dfrac{\partial F}{\partial y} \mathrm d x + \left. \eta\dfrac{\partial F}{\partial y’}\right|_0^{x_0} - \int_0^{x_0}\eta \mathrm d \dfrac{\partial F}{\partial y’} \\
&=\int_0^{x_0}\eta\Big(\dfrac{\partial F}{\partial y} - \dfrac{\mathrm d}{\mathrm d x }\dfrac{\partial F}{\partial y’}\Big) \mathrm d x + \color{red}{0} = 0\end{aligned}$$

标红是因为代入上下限始终有 $\eta(0) = \eta(x_0) = 0$。

然后我们考虑一个问题：当满足 $\eta(x_1) = \eta(x_2) = 0$ 时，积分 $\displaystyle\int_{x_1}^{x_2}\eta(x)f(x)\mathrm d x=0$ 会导出什么结论。

由于 $\eta(x)$ 随意，不妨构造为 $\eta(x) = -f(x)(x - x_1)(x - x_2)$ 以满足题目条件。此时可以有 $\displaystyle \int_{x_1}^{x_2} -[f(x)]^2(x - x_1)(x - x_2)\mathrm d x$。

因为 $[f(x)]^2 \geq 0;\forall x \in [x_1,x_2],-(x-x_1)(x- x_2) \geq 0$，如果存在值使得 $[f(x)]^2$ 大于零，则会直接导致算出的积分值大于零，进一步与已知的函数值产生矛盾，所以 唯一的可能 只有 $[f(x)]^2 = f(x) = 0$。

那么回归到刚刚的讨论，此时就可以有 $\dfrac{\partial F}{\partial y} - \dfrac{\mathrm d}{\mathrm d x }\dfrac{\partial F}{\partial y’} = 0$。

三

再来看一看 最优解所构成的多元函数 对于 $x$ 的全微分。

$$\dfrac{\mathrm d F}{\mathrm d x} = \dfrac{\partial F}{\partial x} + \dfrac{\partial F}{\partial y}\dfrac{\partial y}{\partial x} + \dfrac{\partial F}{\partial y’}\dfrac{\partial y’}{\partial x} = \dfrac{\partial F}{\partial x} + \dfrac{\partial F}{\partial y}y’ + \dfrac{\partial F}{\partial y’}y’'$$

注意到 $\dfrac{\mathrm d}{\mathrm d x}\Big( y’\dfrac{\partial F}{\partial y’}\Big) = \color{red}{y’'\dfrac{\partial F}{\partial y’}} + y’\dfrac{\mathrm d}{\mathrm d x}\dfrac{\partial F}{\partial y’}$，那么上式减左式得到

$$\dfrac{\mathrm d}{\mathrm d x}\Big(F - y’\dfrac{\partial F}{\partial y’}\Big) = \dfrac{\partial F}{\partial x} + y’\Big( \dfrac{\partial F}{\partial y} - \dfrac{\mathrm d}{\mathrm d x}\dfrac{\partial F}{\partial y’}\Big)$$

根据第二节的结论，此时有 $\dfrac{\mathrm d}{\mathrm d x}\Big(F - y’\dfrac{\partial F}{\partial y’}\Big) = \dfrac{\partial F}{\partial x}$。

当整个函数没有出现 $x$ 的时候，$\dfrac{\mathrm d}{\mathrm d x}\Big(F - y’\dfrac{\partial F}{\partial y’}\Big) = 0$。这个时候就有 $F - y’\dfrac{\partial F}{\partial y’} = C$。

四

将上一节算得的结论运用到 $t=\displaystyle\int_0^{x_0}\sqrt{\dfrac{ 1 + (y’)^2}{2gy}}\mathrm{d}x$ 中的 $\sqrt{\dfrac{ 1 + (y’)^2}{2gy}}$ 就有

$$
\sqrt{\dfrac{ 1 + (y’)^2}{2gy}} - y’\cdot\dfrac{2y’}{\sqrt{ 2gy}}\cdot\dfrac{1}{2\sqrt{ 1 + (y’)^2}}=\dfrac{1+(y’)^2 - (y’)^2}{\sqrt{ 2gy( 1 + (y’)^2)}} = C
$$

于是我们得到微分方程 $y( 1 + (y’)^2) = \dfrac{C}{2g}$。将 $\dfrac{C}{2g}$ 设为 $k$，移项开方整理得到

$$\dfrac{\mathrm d y}{\sqrt{\dfrac{k-y}{y}}} = \mathrm d x$$

做变量代换 $y = \dfrac{k}{2}(1-\cos\theta)$ 得到

$$\dfrac{k}{2}\sqrt{\dfrac{1-\cos\theta}{1+\cos\theta}}\sin\theta\mathrm d \theta = \dfrac{k}{2}(1-\cos\theta)\mathrm d \theta = \mathrm d x$$

进而有 $x = \dfrac{k}{2}(\theta-\sin\theta) + C$。

因此 $\begin{cases}x = \dfrac{k}{2}(\theta-\sin\theta) + C \\ y = \dfrac{k}{2}(1-\cos\theta)\end{cases}$，其中更具体的表达式需要由 $(0,0),(x_0,y_0)$ 确定。

对于 $(0,0)$，此时有 $\begin{cases}x = \dfrac{k}{2}(\theta-\sin\theta) \\ y = \dfrac{k}{2}(1-\cos\theta)\end{cases}$。而对于另一个坐标而言，既然已经找到了表达式，就不必多言了。

悬链线

在悬挂问题中，悬挂点是固定的两个点，绳子可看成经过这两个点且总长为 $L$ 的曲线。

根据能量最低原理，绳子对应曲线的形状将会使得绳子的质心最低，而曲线可用一个函数表示。

那么问题就变成求解函数的表示。

我们以地面水平方向为 $x$ 轴，两悬挂点水平方向中垂线为 $y$ 轴，则两悬挂坐标分别为 $-\dfrac{a}{2}$ 和 $\dfrac{a}{2}$ ，绳子对应的曲线函数为 $f(x)$，绳长 $L$，质量为 $m$，则其质心的 $y$ 坐标值 $y_\text{center} $ 不难根据质心方程列出¹

$$
y_\text{center} = \dfrac{\displaystyle\int_{0}^{L}f(x)\cdot\dfrac{m}{L} \mathrm{d}{s}}{m} = \dfrac{1}{L}\int_{-\frac{a}{2}}^{\frac{a}{2}} f(x)\sqrt{1+[f’(x)]^2} \mathrm{d}{x}
$$

沿用先前算得的结论，即可得到以下的微分方程。

$$f(x)\sqrt{1+[f’(x)]^2} - f(x)\dfrac{[f’(x)]^2}{\sqrt{1+[f’(x)]^2}}=\dfrac{f(x)}{\sqrt{1+[f’(x)]^2}}=C$$

进一步有

$$\dfrac{\mathrm d f(x)}{\mathrm d x}=\sqrt{\dfrac{[f(x)]^2 - C}{C}}$$

从而

$$\dfrac{\mathrm d f(x)}{\sqrt{[f(x)]^2 - C}}=\dfrac{\mathrm d x}{\sqrt{C}}$$

后面就不多解释了。

参考资料

悬链线 ↩

5 评论

99的手速加上1的运气 6天前

谢谢，最近想用代码模拟球滚动时每个0.01秒后的速度，但有点难写

Ch4ot1c_M@dn3ss 5天前

做一道题之前一定要学好预备知识并做好分析。不然自然会觉得难。

由前面的阐述现在知道 $\color{red}{x_0,y_0 > 0};k = \dfrac{C}{2g}$ 以及
$$\mathrm d t=\dfrac{\mathrm d s}{v} = \sqrt{\dfrac{\Big(\frac{k}{2}(1-\cos\theta)\Big)^2+\Big(\frac{k}{2}\sin\theta\Big)^2}{2g\cdot\frac{k}{2}(1-\cos\theta)}}\mathrm d \theta = \sqrt{\dfrac{k}{2g}}\mathrm d \theta = \dfrac{\sqrt{C}}{2g}\mathrm d \theta$$

从而有
$$\Delta t = \dfrac{\sqrt{C}}{2g} \Delta \theta = 0.01\Rightarrow \Delta \theta = \dfrac{0.02g}{\sqrt{C}}$$

所以现在的问题变成求 $C$。

$$\begin{cases}x=\dfrac{C}{4g}(\theta - \sin\theta) \\
y = \dfrac{C}{4g}(1-\cos\theta)
\end{cases} \Rightarrow
\begin{cases}x_0=\dfrac{C}{4g}(\theta^* - \sin\theta^*) \\
y_0 = \dfrac{C}{4g}(1-\cos\theta^*)
\end{cases}$$

由于我们希望 $\theta^*$ 也是已知量，不妨先算 $\theta^*$。此时可以有

$$ \theta^* -\sin\theta^* - \dfrac{x_0}{y_0}(1-\cos\theta^*) = 0 $$

直接算肯定算不出来，但这个式子可以用牛顿迭代或者二分计算数值解。

需要注意的是，因为 $x_0,y_0 > 0$，又当 $\dfrac{x_0}{y_0}$ 很小时（即 $\dfrac{x_0}{y_0} < 1$ 时），另一个零点会很靠近 $(0,0)$，而很大时零点个数就相应增多，但方便起见，只考虑 $0$ 以后的第一个解。具体来讲——摆线是周期函数。要 $0$ 以后的第一个根就相当于只下滑一次，不会出现下滑后减速又再次下滑的情况。

不妨从此处开始使得 $\lambda = \dfrac{x_0}{y_0}$。当 $x=0, x- \sin x - \lambda(1-\cos x) = 0$ 总是成立的。

因为在 $x = 2\pi $ 时，$2\pi - \lambda(1-\cos (2\pi)) = 2\pi > 0$。

所以只要找到函数值是负的情况就能说明函数在区间 $(0, 2\pi)$ 里还有一个零点。

我们的目的是找到 $x$ 使得 $0 < x-\sin x < \lambda(1-\cos x)$，也即 $\dfrac{x-\sin x}{1-\cos x} < \lambda$ 成立。

注意到（然而高考之类的考试如果真考了不要这样写，单调性之类的证明还是要的）这个函数单是在 $x\in (0,2\pi)$ 时，值域都已经是 $(0, +\infty)$，那么必然存在值使得 $\dfrac{x-\sin x}{1-\cos x} < \lambda$ 成立。

从而能保证 $\theta^*$ 是已知的。

对于牛顿迭代，有

$$\theta_{n}^* = \theta_{n-1}^* - \dfrac{\theta_{n-1}^* -\sin\theta_{n-1}^* - \lambda(1-\cos\theta_{n-1}^*)}{1-\cos\theta_{n-1}^* - \lambda\sin\theta_{n-1}^*}$$

但第一个极值点还需要通过用程序求导确定。此时有

$$\begin{aligned}&1 - \cos\theta - \lambda\sin\theta = 0 \\& \Rightarrow \csc\theta-\cot \theta = \tan\dfrac{\theta}{2} = \lambda
\\ &\therefore \theta_\text{init} = 2 \arctan\lambda
\end{aligned}$$

那么取初值为 $\pi + \arctan\lambda$ 是比较理想的。

而 $C=\dfrac{4gy_0}{1-\cos\theta^*},\Delta \theta = 0.01\sqrt{\dfrac{g(1-\cos\theta^*)}{y_0}}$。

还有 $\begin{cases}x = \dfrac{y_0}{1-\cos\theta^*}(\theta - \sin\theta)\\
y = \dfrac{y_0}{1-\cos\theta^*}(1-\cos\theta)
\end{cases}$。

每次要算速度，如果是大小，就用上面的 $\Delta \theta$ 带到 $v = \sqrt{2gy} = \sqrt{\dfrac{2gy_0}{1-\cos\theta^*}}\cdot\sqrt{1-\cos\theta}$ 中变化的 $\theta$ 就行。

方向就沿着切线方向，同样你可以通过链式法则算出来 $k_l=\dfrac{\frac{\mathrm d y}{\mathrm d \theta}}{\frac{\mathrm d x}{\mathrm d \theta}} = \dfrac{\sin\theta}{1-\cos\theta} \Rightarrow l: \dfrac{y - y(\theta)}{\sin\theta} = \dfrac{x-x(\theta)}{1-\cos\theta}$。或者写成方向向量的形式 $(1-\cos\theta,\sin\theta)$。

请注意，算出来的与角度有关的结果都是 弧度制。检查的时候注意单位以防止浪费时间。

Ch4ot1c_M@dn3ss 5天前

基于上一条评论描述，一种可行（？的用来计算每 $0.01s$ 速度变化情况的 python 3 代码如下。

"""最速降线计算每 0.01s 的速度增量"""

# 作者 Ch4ot1c_M@dn3ss
# 编写时间 2023.05.27

import math as m_
import matplotlib.pyplot as plt


def _ceil(_: float) -> int:
    return m_.ceil(_)


def _sin(_: float) -> float:
    """去除包头后的 sin 函数"""
    return m_.sin(_)


def _cos(_: float) -> float:
    """去除包头后的 cos 函数"""
    return m_.cos(_)


def _1_cos(_: float) -> float:
    """重复使用的 1 - cos(x)"""
    return 1.0 - _cos(_)


def _sqrt(_: float) -> float:
    """去除包头后的 sqrt 函数"""
    return m_.sqrt(_)


def div(__: float, _: float) -> float:
    """浮点数除法"""
    return __ / _


pi = 3.14159265358979323846264338327950288419716939937510
iterator_count: int = 64                                                # 64 次迭代限制
acceleration_g: float = 9.78046                                         # 重力加速度大小
des_x, des_y = map(float, input().split())                              # 终点坐标需要是 **正数**
try:
    assert (des_y > 0 and des_x > 0)
except AssertionError:
    print("Needed: x > 0, y > 0.")
    exit(0)
ratio = des_x / des_y


def newton_iteration(_value: float) -> float:
    global iterator_count, ratio
    __n_, cnt, __n = _value, 0, 0.0
    while cnt <= iterator_count:
        __n = __n_ - div(__n_ - _sin(__n_) - ratio * _1_cos(__n_), _1_cos(__n_) - ratio * _sin(__n_))
        __n_, cnt = __n, cnt + 1

    return __n


def init_theta(Lamb: float) -> float:
    def arc_tan(_alpha: float) -> float:
        return m_.atan(_alpha)
    return 2 * arc_tan(Lamb)


theta_ = newton_iteration(init_theta(ratio) / 2.0 + pi)                 # 以初值迭代，以保证不会找到 (0, 0)
print(f'The ideal Zero point, or the Theta-parameter of destination is {theta_}.')
Constant = div(4 * acceleration_g * des_y, _1_cos(theta_))              # 计算得到的常量
print(f"c / 4g = {div(Constant, 4 * acceleration_g)}\n")                # 查验用
d_theta = 0.01 * _sqrt(div(acceleration_g * _1_cos(theta_), des_y))     # 每次偏移角度

const, i = _sqrt(div(2 * acceleration_g * des_y, _1_cos(theta_))), 0    # 计算速度时需要的常量

velocity = [const * _sqrt(_1_cos(d_theta * i)) for i in range(0, _ceil(theta_ / d_theta))]
time_counter = [i * 0.01 for i in range(0, _ceil(theta_ / d_theta))]

for v in velocity:
    print('%.8lf, ' % v, end="")                                        # 保留八位小数输出
    i += 1
    if i % 10 == 0:
        print('')

fig1, = plt.plot(time_counter, velocity, "ob:", label='')               # 画图
plt.legend(handles=[fig1], labels=['v-t'], loc='best')
plt.show()

这个答复所用到的资料是这篇文章。一些细节或者疑问还需要你自己看了~ 我只能帮到这么多。

Ch4ot1c_M@dn3ss 5天前

等一下，我做别的事情突然想起来你的要求是滚动，那这的确太难了。但我还是希望这个滑动的解法能对你有一定帮助。

99的手速加上1的运气 4天前回复了 Ch4ot1c_M@dn3ss 的评论

谢谢大佬！