如何求解两点间的最速降线

两点之间的最短路线不一定是直线。如果将用时最短作为从 A 点到达 B 点的最短路线评判标准，那么在重力作用下，高度不同的两点间的最短路线就是最速降线。在本篇博客文章中，我们将演示如何使用 COMSOL Multiphysics 的内置数学表达式和“优化模块”来求解最速降线。

寻找最速降线

想象一下：在弧形坡道（类似滑板公园）上滚下一块大理石，然后测定大理石从 A 点滚到B 点所需的时间。我们的目标是重新设计两点之间的坡道形状，使大理石的运动时间为最短。为了简单起见，我们假设一个理想情况：忽略摩擦，且将大理石视为无限小（一个质点）。

最速降线 是一条理想曲线，物体在该曲线上的下降速度最快。实际上，这种情况存在一个对应的解析解，或者,通过一系列推导工作，利用 COMSOL Multiphysics 的“偏微分方程”功能来求解此类问题。然而，毕竟我是“最小行动原则的虔诚信徒”——这是对“懒惰的物理学家”的一种“美称”——所以我们在本文中选择使用“优化模块”。

创建 COMSOL Multiphysics 模型

假设问题中的 A 点坐落于原点 (x, y) = (0, 0)，那么最速降线的解析解是一条具有下列形式的参数化曲线：

(1)

其中 r 参数是一个常数，参数 t 是参数化曲线的工作参数，在曲线上的 tA 与 tB 点之间呈线性变化。在求解该问题时，我们通常已知 B 点的位置，待求解 r 和 t。

在本文中，我们将从最速降线的解析解和一组已知 r 和 t 值入手，根据后者，我们可以推导出 B 点的位置。我们将展示如何使用 COMSOL Multiphysics 的优化功能和“优化模块”来计算解析解的近似值。

我们先建立一个空模型。事实上，我们不会在“模型开发器”中添加任何“组件”，也完全不需要几何、物理场或网格。这是一个有趣的“没有模型的模型”示例。

首先，我们定义一组全局参数。将常数参数 r 设为 1，将 B 点上的 tB 值设为 1.2 \pi。（因为 A 是原点，所以 tA=0）然后可以使用方程（1）计算 B 点的位置（xB，yB），如下方截图所示。

接着，我们使用插值函数来拟合最速降线。我们将几个插值点的 x 位置定义为 x1 ~ x4，并为 y 位置（ y1 ~ y4）给出相应的初始猜测值，使插值点开始时位于点 A 到点 B 之间的直线上。

这个插值函数可按照下方截图进行设置。

单击创建绘图 按钮后，就会生成插值函数绘图。为了方便读者看清，我们删除了两个外推插值图。随后，将解析解的线图添加到同一个绘图组中，以便稍后与数值解进行比较。

要创建参数化曲线图，首先在“全局定义”节点下创建一个虚拟解析函数，将变元的上限设置为 tB（“绘图参数”栏中）。该解析函数的唯一作用是提供 t 在 0 到 tB 之间的数值列表，据此绘制参数化曲线。单击创建绘图，然后将生成的线图 1 从一维绘图组 2 拖拽到一维绘图组 1 中（其名称会自动更改为“线图 2”）。在 y 轴数据中输入方程（1）的对应公式：r*(-1+cos(root.x))。同样地，在 x 轴数据中输入 r*(root.x-sin(root.x))。

请注意：在此例中，COMSOL Multiphysics 变量 root.x 对应方程（1）的参数 t。

单击绘图 按钮，软件立刻生成了解析解（绿色曲线）和初始猜测值（蓝色曲线和黑点）图表，如下所示。

计算行程时间

为了计算大理石从 A 点滚到 B 点所花费的时间，我们假设小球运动时无摩擦，从而使全部势能损失转化为动能，而动能与速度的平方成正比。所以如果使用公式 y = f(x) 来描述曲线，则瞬时速度与高度方向上损失的平方根成正比：v \propto \sqrt{0-f(x)}（上文假设了 A 点的原始高度为 0）。至于 dx 在 x 方向上的无穷小运动，大理石沿着曲线行进的路径长度为 ds = dx \sqrt{1+f'(x)^2}。通过这段长度所需的时间等于长度除以速度：d\tau = ds/v。因此，对于任何给定曲线 y = f(x)，我们得出了总行程时间的表达式：

现在我们只需要让 COMSOL 软件找到用时最短的那条曲线。

您或许注意到了，我们在表达式中使用了“正比”符号（\propto）来表示行程时间，这是因为我们忽略了公式中的大理石质量和重力加速度常数。这些数字组合在一起，只会基于常数因子缩短或增长运动时间，并不会影响最小化问题中的曲线形状。换句话说，最速降线完全不受大理石的重量的影响。

因为我们正在使用插值函数 int1 来拟合曲线 f(x)，因此可以使用上方公式对表示总行程时间的全局变量 T 进行定义：integrate(sqrt((1+(d(int1(x),x))^2)/max(0-int1(x),eps)),x,0,xB)。如下所示，分母中的其余表达式 max(... ,eps) 排除了被零除的错误情况。

优化求解器

现在，我们就要“命令”软件去最大程度地减少行程时间。首先创建一个空研究，然后单击右键，在研究下添加一个“优化”节点。我们可以使用“坐标搜索”、Nelder-Mead、BOBYQA 或 COBYLA 优化求解器来对付这个“无模型”优化问题。事实表明，COBYLA 的求解速度最快。

打开设置窗口后，我们在“目标函数”栏的表达式框中输入 T，它被默认设为最小化。然后，在“控制变量和参数”栏下，单击添加 按钮（“加号”图标），并添加参数 y1 - y4。软件会在初始值字段中自动填入全局参数列表中的对应值。接着，在“求解时输出”栏下，勾选绘图 复选框。这些设置均显示在下方截图中。

单击计算 后，观察优化求解器不断上下移动插值曲线，直至获得最小行程时间。即使我们使用的是只有四个插值点的简单线性插值曲线，下图中的结果也与解析解非常吻合。使用高阶插值和添加更多插值点能够进一步改进求解结果。

结语

我的导师 William Vetterling 就职于 ZINK 成像公司（ZINK Imaging），数年前他曾与我们共进午餐，一起谈论了怎样使用 COMSOL Multiphysics 去求解任何问题，因为它提供的数学工具可以处理科技研究中的各类方程。他后来基于这个想法在 COMSOL 用户年会 2012 波士顿站发表了题为“巴别图书馆”的主题演讲。我们举出了一个“没有模型的模型”案例，通过使用软件内置的数学表达式和优化功能来求解最速降线问题。我希望这项研究能够激发 COMSOL Multiphysics 的创造性应用，帮助人们应对更多的技术难题。