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

2015年 10月 20日

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

寻找最速降线

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

曲线示意图。

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

创建 COMSOL Multiphysics 模型

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

(1)

\begin{align*}
x(t) &=& r (t-sin(t))\\
y(t) &=& r (-1+cos(t))
\end{align*}

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

在本文中,我们将从最速降线的解析解和一组已知 rt 值入手,根据后者,我们可以推导出 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)。至于 dxx 方向上的无穷小运动,大理石沿着曲线行进的路径长度为 ds = dx \sqrt{1+f'(x)^2}。通过这段长度所需的时间等于长度除以速度:d\tau = ds/v。因此,对于任何给定曲线 y = f(x),我们得出了总行程时间的表达式:

T = \int{d\tau}=\int{ds/v}\propto\int_0^{xB}{\sqrt{{1+f'(x)^2}\over{0-f(x)}}dx}

现在我们只需要让 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。软件会在初始值字段中自动填入全局参数列表中的对应值。接着,在“求解时输出”栏下,勾选绘图 复选框。这些设置均显示在下方截图中。

COMSOL Multiphysics 中的优化设置。

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

绘图显示了优化结果。

结语

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

博客分类


评论 (2)

正在加载...
航 李
2018-07-16

这个博客我按着教程来 为什么后面勾选了‘求解时输出’选项提示为‘扫描错误,未知属性:属性 soulation’ 这个原因是什么呀?

Tengyue Gao
2018-10-19

李航,您好!
感谢您的评论。
模型相关的问题,请您联系我们的技术支持团队:
在线支持中心:cn.comsol.com/support
Email: support@comsol.com
谢谢!

浏览 COMSOL 博客