由最速降线问题引出非线性问题的通用解法

The Generalized Solutions to Nonlinear Problems Induced by the Most Rapid Curved Descending Line Problem

摘要: 针对最速降线的问题，首先在不存在摩擦时的情况下推导出适用于求此类非线性微分方程极值的变分法，将其应用于在现实里用拉格朗日乘子法求解库仑摩擦最速降线的过程中，随后在给定的边界条件下应用打靶法结合牛顿法进行快速迭代逼近并利用mathematica建模求得数值解。数值计算在科学研究和工程技术中都起到很重要的作用，而非线性方程组的数值解法是计算数学的一个重要的研究内容。此套方法涉及到以时间和空间等物理量为参数的表达方式，还具有可快速编程代码化的优点。因此可适用于物理实验中求解绝大部分在一定精度要求下的非线性问题。

Abstract: In order to solve the problem of the most rapidly falling line, we first derive the variational method for finding the extrema of such nonlinear differential equations in the absence of friction, apply it to the process of solving the Coulomb frictional most rapidly falling line by Lagrange multiplier method in reality, and then apply the targeting method combined with Newton’s method for fast iterative approximation under the given boundary conditions and use mathematica modeling to find the numerical solution. Numerical computation plays an important role in both scientific research and engineering technology, and the numerical solution of nonlinear systems of equations is an important research element in computational mathematics. This method involves the expression of physical quantities such as time and space as parameters, and also has the advantage of being rapidly programmable and codable. Therefore, it can be applied to solve most of the nonlinear problems in physical experiments under certain accuracy requirements.

1. 引言

最速降线存在于意大利科学家伽利略在1630年提出一个分析学的基本问题：一个质点在重力作用下，从一个给定点到不在它垂直下方的另一点，如果不计摩擦力，问沿着什么曲线滑下所需时间最短。伽利略认为这曲线是圆，可是这是一个错误的答案。瑞士数学家约翰·伯努利在1696年再提出这个最速降线的问题并征求解答。次年已有多位数学家得到正确答案，其中包括牛顿、莱布尼兹、洛必达和伯努利家族的成员。这问题的正确答案是连接两个点上凹的唯一一段旋轮线。旋轮线与1673年荷兰科学家惠更斯讨论的摆线相同，因为钟表摆锤作一次完全摆动所用的时间相等，所以摆线(旋轮线)又称等时曲线，这也在文献 [1] 中有所证明。数学家十分关注最速降线问题，大数学家欧拉也在1726年开始发表有关的论著，在1744年最先给出了这类问题的普遍解法，并产生了变分法这一新数学分支。

争对更具现实意义的存在摩擦时的最速降线问题，文献 [2] 给出了相关解析解，但其中推导过程与结果相当复杂，公式代换繁琐，最终的解析表达式甚至失去了解析解的意义。并且其只争对本文所引入的一种非线性问题进行求解，其中的公式变换技巧只适用于特殊情况，不能应用于一般的非线性问题。文献 [3] 用遗传算法来对曲线进行精度优化，最终算法收敛得到无阻力无初速度、考虑粘滞阻力、考虑摩擦阻力三种不同情况下的最速降线。本文尝试区别于以上方法尝试给出一种新型的、更直观的适合于求解此类非线性问题的通用解法。

2. 不存在摩擦时的最速曲线

2.1. 泛函

想象一个滑块P存在于一点A，速率为零，在受重力作用下沿一条轨道去到一点不高于A的B，求所需时间最少的轨道方程。先从直观上来理解这个问题：因为两点之间直线最短，假设AB之间是一条直线，滑块能够走最短的路径，但是由于一开始的加速度并不大，其全程平均速度始终处于较低值；若一开始让斜坡变陡使滑块的速度变快，而由于小球的滑行距离相应变长，其到达时间也未必减小。因此，最速降线轨迹应介于两条轨迹之间。对此，进行具体数学分析。用图1表示物体运动状态。

Figure 1. The state of motion of the slider

图1. 滑块的运动状态

用s来表示从A到B的弧长，速率v为s对时间t的导数并且整个过程机械能守恒：

${\left(\text{d}s\right)}^2={\left(\text{d}x\right)}^2+{\left(\text{d}y\right)}^2=\sqrt{1+{\left(\frac{\text{d}y}{\text{d}x}\right)}^2\text{d}}x$ (1)

$v=\frac{\text{d}s}{\text{d}t}$ (2)

$\frac{1}{2}m{v}^2=mgy$ (3)

从以上三个约束，设A，B的x坐标分别为 $\alpha$ ， $\beta$ ，可以得到整个过程的时间T的函数为：

$\text{d}t=\frac{\text{d}s}{\sqrt{2gy}}=\sqrt{\frac{1}{2gy}\left[1+{\left(\frac{\text{d}y}{\text{d}x}\right)}^2\right]}\text{d}x\Rightarrow T\left(t\right)={\displaystyle {\int }_{\alpha }^{\beta }\sqrt{\frac{1}{2gy}\left[1+{\left(\frac{\text{d}y}{\text{d}x}\right)}^2\right]}}$ (4)

可知，T是x，y以及y'的函数，即：

$T={\displaystyle {\int }_{\alpha }^{\beta }F}\left(x,y,y^{\prime }\right)\text{d}x$ (5)

对于每一种y，y'的形式，都存在一个时间T，时间T只和y，y'的形式有关；T定义了一个函数到实数的映射，数学上把T称为y，y'的泛函。

2.2. 欧拉–拉格朗日方程的引入

假设存在一个y能使T取得最小值并其存在一阶导数，对这个曲线附加一个轻微的扰动，曲线的形式有了轻微的改变，质点下滑的总时间将会增加，同时，微扰越小，当微扰趋近于0的时候，总时间的改变也应该趋近于0，在两个端点处，微扰的值应该恒为0。这个思想和微积分当中记号d的思想很类似，但是有一个显著的区别，微扰对应的其实是一族函数，为了区别，把它记为 $\delta$ 。这就是用于求解泛函极值的变分，表达式为：

$\delta T=T\left(y+\delta y\right)-T\left(y\right)$ (6)

用(5)式将(6)式展开得：

$\delta T={\displaystyle {\int }_{\alpha }^{\beta }\left[F\left(x,y+\delta y,y^{\prime }+\delta y^{\prime }\right)-F\left(x,y,y^{\prime }\right)\right]\text{d}x}$ (7)

根据定义可证明微分运算和变分运算是可以互换顺序的 [4] 。即：

$\text{d}\left(\delta y\right)=\delta \left(\text{d}y\right)$ (8)

套用微积分公式：

$\begin{array}{c}F\left(x,y+\delta y,y^{\prime }+\delta y^{\prime }\right)=F\left(x,y,y^{\prime }\right)+\frac{\partial F\left(x,y,y^{\prime }\right)}{\partial y}\delta y+\frac{\partial F\left(x,y,y^{\prime }\right)}{\partial y^{\prime }}\delta y^{\prime }\\ =F\left(x,y,y^{\prime }\right)+\frac{\partial F\left(x,y,y^{\prime }\right)}{\partial y}\delta y+\frac{\partial F\left(x,y,y^{\prime }\right)}{\partial y^{\prime }}{\left(\delta y\right)}^{\prime }\end{array}$ (9)

得到：

$\delta T={\displaystyle {\int }_{\alpha }^{\beta }\left(\frac{\partial F}{\partial y}\delta y+\frac{\partial F}{\partial y^{\prime }}{\left(\delta y\right)}^{\prime }\right)\text{d}x}$ (10)

对(10)项右部分第二项套用分布积分法则，端点处得变分为0，得：

${\displaystyle {\int }_{\alpha }^{\beta }\frac{\partial F}{\partial y^{\prime }}{\left(\delta y\right)}^{\prime }\text{d}}x={\frac{\partial F}{\partial y^{\prime }}\delta y|}_{x=\alpha }^{\beta }-{\displaystyle {\int }_{\alpha }^{\beta }\left(\frac{\text{d}}{\text{d}x}\frac{\partial F}{\partial y^{\prime }}\delta y\right)\text{d}x}=-{\displaystyle {\int }_{\alpha }^{\beta }\left(\frac{\text{d}}{\text{d}x}\frac{\partial F}{\partial y^{\prime }}\delta y\right)\text{d}x}$ (11)

当T取最小值时(10)式为0：

$\delta T={\displaystyle {\int }_{\alpha }^{\beta }\left(\frac{\partial F}{\partial y}-\frac{\text{d}}{\text{d}x}\frac{\partial F}{\partial y^{\prime }}\right)\delta y\text{d}x}=0$ (12)

由于 $\delta y$ 是任意选取的，要保证(12)式为0则需括号内部为0，则：

$\frac{\partial F}{\partial y}-\frac{\text{d}}{\text{d}x}\frac{\partial F}{\partial y^{\prime }}=0$ (13)

这就是适用于求解泛函极值问题的适用条件——欧拉–拉格朗日方程。在套用公式解这类问题之前还应证明其解的必要性、存在性、唯一性和充分性 [5] 。

3. 存在摩擦时的最速降线

3.1. 模型一般化

以上的求解过程是基于无摩擦阻力假设之上的，在机械工程、设施建设、物理实验等实际情况中，往往应考虑摩擦阻尼对理想状况产生的影响，假设滑块所受到的摩擦力属于库仑摩擦力 [6] 。其由摩擦因素 $\mu$ 和正压力共同作用，而正压力不仅受到滑块所处轨道角度引起重力分量变化的影响，还受物体做圆弧运动时所受到的向心力的反作用力的影响，两者相加后即为物体给予轨道的正压力F。因此在图1中引入角度 $\theta$ ，将其定义为滑块运动方向与水平方向的夹角。同时任意化A点和B点的位置做出图2。

Figure 2. General most rapid descent line

图2. 一般性最速降线

3.2. 拉格朗日乘子法的引入

由上文可知，每一种y，y'对应一个函数值t，于是可将函数视为三维函数。以图3三维函数视图和图4等高线视图为例：

Figure 3. Three bit function view

图3. 三维函数视图

Figure 4. Contour view

图4. 等高线视图

设目标函数F为图2函数的高度h并受到函数 $g\left(y,y^{\prime }\right)=0$ 的约束，g函数始终处于三维函数表面上，那么函数的极值则存在于众多与g函数相切的等高线的高度h中。设极值点为M，将上述用数学语言表示成：

$\{\begin{array}{l}{\text{d}f|}_M=\nabla f=\frac{\partial f}{\partial y}\text{d}x+\frac{\partial f}{\partial y\prime }\text{d}y=0\\ {\text{d}g|}_M=\nabla g=\frac{\partial g}{\partial y}\text{d}x+\frac{\partial g}{\partial y\prime }\text{d}y=0\end{array}$ (14)

在图3中，等高线 $l_1,l_2,l_3$ 所代表的高度 $h_1,h_2,h_3$ 即为函数 $f=\left(y,y^{\prime }\right)$ 的极小值的可能取值， $M=\min \left(h_1,h_2,h_3\right)$ 。此时从几何意义上讲，函数f的等高线的法向量与约束函数g在切点上的法向量平行且成倍数关系，从数学意义上讲，要让方程组(16)同时成立，方程组的系数必须线性相关，则：

$\nabla f=\lambda \cdot \nabla g$ (15)

$\lambda$ 为拉格朗日乘子。拓广到高维问题：如果还有多于一条的约束条件限制，那么函数 $f\left(x,y_1,y_2,\cdots ,y_m\right)$ 取得极值M时，会与所有约束条 $g_j\left(x,y_1,y_2,\cdots ,y_m\right)=0$ 所描述的曲面相切，并且方程组中的每一个方程与其他方程的系数都线性相关，由此推导出辅助函数F：

$\nabla f={\displaystyle \underset{j}{\sum }{\lambda }_j\nabla g_j}\Rightarrow F\equiv f-{\displaystyle \underset{j}{\sum }{\lambda }_j{g}_j}$ (16)

定义辅助函数的目的是能将原有的约束条件自然的隐藏在辅助函数F中，则最初的条件极值问题转化为求辅助函数F的极值问题。

3.3. 建模

首先推导在摩擦力约束下的运动方程，滑块在切向和法向的运动方程为：

$m\frac{\text{d}v}{\text{d}t}+mg\sin \theta +\mu F_n=0$ (17)

$mv\frac{\text{d}\theta }{\text{d}t}+mg\cos \theta -F_n=0$ (18)

联立得到：

$\frac{\text{d}v}{\text{d}t}+\mu v\frac{\text{d}\theta }{\text{d}t}+g\left(\sin \theta +\mu \cos \theta \right)=0$ (19)

另外存在两个速度约束下的运动方程：

$\frac{\text{d}x}{\text{d}t}-v\cos \theta =0$ (20)

$\frac{\text{d}y}{\text{d}t}-v\sin \theta =0$ (21)

随即结合拉格朗日乘子法定义拉式量L [7] ：

$L=1+\alpha \left(\frac{\text{d}x}{\text{d}t}-v\cos \theta \right)+\beta \left(\frac{\text{d}y}{\text{d}t}-v\sin \theta \right)+\gamma \left[\frac{\text{d}v}{\text{d}t}+\mu v\frac{\text{d}\theta }{\text{d}t}+g\left(\sin \theta +\mu \cos \theta \right)\right]$ (22)

其中 $\alpha ,\beta ,\gamma$ 为拉格朗日乘子，定义时间泛函：

$T={\displaystyle {\int }_0^{T}L\text{d}t}$ (23)

记自变量 $q=\left\{x,y,\theta ,v,\alpha ,\beta ,\gamma \right\}$ ，再欧拉–拉格朗日方程(13)得到 $\frac{\partial L}{\partial q}-\frac{\text{d}}{\text{d}t}\frac{\partial L}{\partial \frac{\text{d}q}{\text{d}t}}=0$ 由此产生七个方程：

$\frac{\text{d}\alpha }{\text{d}t}=0$ (24)

$\frac{\text{d}\beta }{\text{d}t}=0$ (25)

$\mu \left(\gamma \frac{\text{d}v}{\text{d}t}+v\frac{\text{d}\gamma }{\text{d}t}\right)-v\left(\alpha \sin \theta -\beta \cos \theta \right)-g\gamma \left(\cos \theta -\mu \sin \theta \right)=0$ (26)

$\gamma \mu \frac{\text{d}\theta }{\text{d}t}-\frac{\text{d}\gamma }{\text{d}t}-\left(\alpha \cos \theta +\beta \sin \theta \right)=0$ (27)

$\frac{\text{d}x}{\text{d}t}-v\cos \theta =0$ (28)

$\frac{\text{d}y}{\text{d}t}-v\sin \theta =0$ (29)

$v\mu \frac{\text{d}\theta }{\text{d}t}+\frac{\text{d}v}{\text{d}t}+g\left(\sin \theta +\mu \sin \theta \right)=0$ (30)

现在考虑起始位置的边界条件：

$\{\begin{array}{l}x\left(0\right)=x_0\\ y\left(0\right)=y_0\\ v\left(0\right)=0\end{array}$ (31)

和末位置的边界条件：

$\{\begin{array}{l}x\left(T\right)=x_1\\ y\left(T\right)=y_1\end{array}\Rightarrow \{\begin{array}{l}x\left(T\right)-x_1=0\\ y\left(T\right)-y_1=0\end{array}$ (32)

同时在此类可动边界的变分问题上还应附加横截条件 [2] ：

$\{\begin{array}{l}\gamma \left(T\right)=0\\ 1+g\gamma \left(\sin \theta +\mu \cos \theta \right)-v\left(\alpha \cos \theta +\beta \sin \theta \right)=0\end{array}$ (33)

3.4. 打靶法结合牛顿法

如上所示的边值问题要想直接求解会非常困难，因此使用打靶法将边值问题转为初值问题进行求解，在猜测初始值时，方程(26)，(27)，(30)可以联立得到一个只含 $\theta$ 的代数方程，因此 $\theta$ 由其他参数所制约，于是只需要猜测初值 $T,{\alpha }_0,{\beta }_0,{\gamma }_0$ ，随后将原方程组(24)至(30)按初值问题进行求解。首先记多项式组 $F=\left\{x-x_1,y-y_1,\gamma ,1+g\gamma \left(\sin \theta +\mu \cos \theta \right)-v\left(\alpha \cos \theta +\beta \sin \theta \right)\right\}$ 。

注意此时 $F=F\left[T,{\alpha }_0,{\beta }_0,{\gamma }_0\right]$ ，若初值问题的解：

$\{\begin{array}{l}x\left(T\right)\hfill \\ y\left(T\right)\hfill \\ \gamma \left(T\right)\hfill \\ 1+g\gamma \left(\sin \theta +\mu \cos \theta \right)-v\left(\alpha \cos \theta +\beta \sin \theta \right)\hfill \end{array}$ (34)

在 $\{\begin{array}{l}T=T_0\hfill \\ \alpha =a_0\hfill \\ \beta ={\beta }_0\hfill \\ \gamma ={\gamma }_0\hfill \end{array}$ 的值满足 $F\left[T_0,{\alpha }_0,{\beta }_0,{\gamma }_0\right]=0$ 或 $\left|F\left[T_0,{\alpha }_0,{\beta }_0,{\gamma }_0\right]\right|<\epsilon$ ，其中 $\epsilon$ 为允许的误差界限。这样，方程(34)即为边值问题的近似解。但往往 $F\left[T_0,{\alpha }_0,{\beta }_0,{\gamma }_0\right]\ne 0$ ，对 $T_0,a_0,{\beta }_0,{\gamma }_0$ 做适当调整。为此，可以采用解非线性方程组时的牛顿迭代格式 [8] 来快速逼近。此时迭代格式为：

$\left\{T_{i+1},{\alpha }_{i+1},{\beta }_{i+1},{\gamma }_{i+1}\right\}=\left\{T_i,{\alpha }_i,{\beta }_i,{\gamma }_i\right\}-J_i^{-1}\cdot F_i$ (35)

其中 $J_i$ 为雅克比矩阵 [9] ， $J=\nabla F$ 。在实际计算雅克比矩阵时应利用差分 [10] 去替代导数。继续令 $\{\begin{array}{l}T=T_1\hfill \\ \alpha =a_1\hfill \\ \beta ={\beta }_1\hfill \\ \gamma ={\gamma }_1\hfill \end{array}$ ，再求解此初值问题，计算得到它的解为：

$\{\begin{array}{l}x\left(T_1\right)\hfill \\ y\left(T_1\right)\hfill \\ \gamma \left(T_1\right)\hfill \\ 1+g{\gamma }_1\left(\sin \theta +\mu \cos \theta \right)-v\left(a_1\cos \theta +{\beta }_1\sin \theta \right)\hfill \end{array}$ (36)

由于曲线终点为 $\left(x_1,y_1\right)$ ，为了不混淆，将第二次迭代值的下标改为2，此时 $\{\begin{array}{l}x\left(T_1\right)=x_2\hfill \\ y\left(T_1\right)=y_2\hfill \\ \gamma \left(T_1\right)={\text{γ}}_2\hfill \\ 1+g{\gamma }_1\left(\sin \theta +\mu \cos \theta \right)-v\left(a_1\cos \theta +{\beta }_1\sin \theta \right)={\delta }_2\hfill \end{array}$ ，若 $F\left[T_1,{\alpha }_1,{\beta }_1,{\gamma }_1\right]=0$ 或 $\left|F\left[T_1,{\alpha }_1,{\beta }_1,{\gamma }_1\right]\right|<\epsilon$ ，则方程组(34)作为边值问题的近似解；否则，再对 $T_1$ 做适当的调整。如此重复计算，直至 $F\left[T_k,{\alpha }_k,{\beta }_k,{\gamma }_k\right]=0$ 或 $\left|F\left[T_k,{\alpha }_k,{\beta }_k,{\gamma }_k\right]\right|<\epsilon$ 时，便以 $\{\begin{array}{l}x\left(T_k\right)\hfill \\ y\left(T_k\right)\hfill \\ \gamma \left(T_k\right)\hfill \\ 1+g{\gamma }_k\left(\sin \theta +\mu \cos \theta \right)-v\left(a_k\cos \theta +{\beta }_k\sin \theta \right)\hfill \end{array}$ 作为边值问题的近似解。

一般情况下，需要根据具体物理情况猜测迭代初值，若猜测得当，牛顿法可以在一定精度要求下于20步内收敛。上述过程可用图5表示。

Figure 5. Progressive approach process

图5. 逐步逼近过程

下面给出一组具体算例：以知滑块从点 $A\left(-\left(3\pi /2+1\right),1\right)$ 以速度v = 0下滑至点B (0, 0)求当摩擦系数为 $\mu =0,0.05,0.1,0.15,0.17$ 时的最速曲线数值解。利用mathematica，将上述过程写为程序从而得到图6。

Figure 6. Numerical and analytical solutions of the maximum velocity curve for friction coefficients $\mu =0,0.05,0.1,0.15,0.17$

图6. 摩擦系数为 $\mu =0,0.05,0.1,0.15,0.17$ 时的最速曲线数值解与解析解

在图6中，自下往上为摩擦系数逐渐增大的最速降线的情况，这与文献 [3] 在摩擦系数 $\mu$ 为0至0.2的情况相符合，其中曲线为解析解，散点为数值解，可以看到两者重合，由此证实了该方法的可行性。

4. 结论

本文通过最速降线的问题给出了一种非线性问题的通用解法，该方法利用打靶法将边值问题转化为初值问题进行求解，随后引入牛顿法进行快速迭代逼近，在数学软件Mathematica符号计算能力和数值求解能力的双重优秀功能加持下，可快速求解绝大部分非线性问题，其中牛顿法也是一种易于快速编程建模的方法。由于复杂非线性问题的解析解较为难求甚至不可求解，因此此法适合在物理实验等现实领域中在较短时间内提供可操作性强的数值解。