最速降线问题是一个古老的物理和数学相结合的问题,最初由伽利略在1630年提出.伽利略认为最速降线的形状为圆弧,但这是错误的.后来约翰·伯努利、牛顿、莱布尼兹、雅克比·伯努利等人解决了这个问题,结果是旋轮线(摆线).对最速降线的研究很有意义,甚至促进出现了一个新的数学分支——变分法.最速降线符合费马原理,即光传播的路径光程取极值,这也是约翰·伯努利的解决方法.在最速降线问题中,一个纯粹的力学问题和光学中的折射定律产生了密切的联系,体现了物理学的内在和谐之美.最速降线问题的极值思想、违反人直觉的特性(如小球运动轨迹可能包括下降段和上升段),使其广泛应用于各类物理科学演示实验中.在教学中,最速降线能够使学生对物体运动的力学规律和其背后的数学物理原理产生深刻认识.虽然利用变分法可以得到最速降线的解析解,但是如果能对最速降线进行数值求解和模拟,画出各种条件下小球运动的路径轨迹,直观比较小球沿不同曲线下滑的时间,对提高教学效果会有非常大的帮助.

一些文献中对最速降线的数值求解研究,或者比较单一且没有进行详细的讨论,或者物理意义不够突出(如遗传算法、粒子群算法等),算法意义大于物理意义[1,2,3,4,5,6,7,8].基于这些情况,本文使用IgorPro,从4个方面系统介绍了最速降线数值求解的方法:1)解析法,利用摆线方程数值求解并进行求解讨论;2)微分方程法,利用最速降线微分方程进行数值求解并进行求解讨论;3)定义求解法,利用最速降线定义进行数值求解并进行求解讨论;4)通过数值方法利用解析公式求解不同轨迹所用时间,并与最速降线比较.

IgorPro是一个科学软件,主要用于数据处理和可视化.通过对最速降线问题的求解,可使学生掌握利用IgorPro进行图表绘制、数值方程求解、微分方程求解和通过程序设计处理实验数据的方法.

1、最速降线问题

1.1问题描述

在竖直平面上建立坐标系Oxy,在此平面上任意给定两点,一小球(可视为质点)在重力的作用下沿一曲线从一较高点向不在正下方的较低点滑动(忽略摩擦力),问曲线取什么形状时,小球滑行的时间最短.这就是最速降线问题.为了方便,取小球的初始位置为原点,终点为B(x1,y1),如图1所示.

图1最速降线

1.2解析求解分析

如图1,取OB上一点P(x,y),当小球位于P点时,根据能量守恒,小球的速度为

υ=2gy−−−√υ=2gy(1)

g为重力加速度.取P附近一小段ds,则小球通过该小段需要时间dt为

t是y(x)的函数.于是最速降线问题就是要寻找这样一个函数f(x)=y(x),使得式(3)中的积分的值最小.利用欧拉-朗格拉日第二方程

解微分方程(7),得到如下解:

这是摆线的参数方程,其几何意义如图1所示.一半径为r的圆(旋轮)与Ox轴相切,初始切点位置在原点处,当圆沿着Ox轴向右滚动(滚动角度为θ)时,圆上切点处点的轨迹就是最速降线.将y(0)=0,y(x1)=y1代入方程(8),可以解得r和θmax(θmax表示θ的最大值).利用r可以计算出摆线上任意一点的坐标,利用θmax可计算出摆线的最大范围.从图1可以看出,随着B点位置的不同,最速降线除了下降段外,还会出现上升段.是否会出现上升段,取决于y1/x1的比值.由摆线方程(8)可知,摆线最低点处θ=π,此时有y/x=2/π,因此当y1/x1≥2/π时,最速降线只有下降段,当y1/x1<2/π时,还会有上升段.

1.3利用定义数值求解分析

将OB在y方向均匀分为N段,每段高度为h,h=y1/π,如图2所示.

图2最速降线问题数值分析示意图

2、使用IgorPro进行数值求解

2.1使用摆线方程数值求解

使用摆线方程数值求解直接利用了最速降线是摆线的结论,不同于微分方程法和定义法,不具有“第一原理性”.数值求解主要解决两个问题:1)利用(x1,y1)求解摆线方程的参数r和θmax,计算数值解;2)证明解的存在性和唯一性.

1)摆线方程数值求解

将(x1,y1)代入摆线方程(8),有

这是一个二元超越方程组,只能用数值方法求解.由方程组(17)可得:

1−cosθθ−sinθ=x1y11-cosθθ-sinθ=x1y1(18)

方程转化为一元方程,解得θmax后代入式(17)任何一个方程即可得到r.

方程(18)可使用FindRoots命令求解,方法如下:

这里取初始坐标为(0,0),终止坐标为(x1,y1),wBx和wBy里分别存放解的x和y坐标.图3给出了y1/x1取不同值时的解,可以看出当y1/x1=0.3时,最速降线含有上升段,而当y1/x1=1时不包含上升段.这与前面结论是一致的.

图3利用摆线方程数值求解结果。

当y1/x1取值小于2/π时，曲线包括下降段和上升段

2)解的存在性和唯一性

最速降线的解存在且唯一[6],这要求方程(18)总是有解,且解只有一个.为了验证此结论,可构造如下函数:

f(x)=x−sinx1−cosxx∈[0,2π]f(x)=x-sinx1-cosxx∈[0,2π](19)

利用下面的方法绘制f(x)曲线:

结果如图4所示,可以看出f(x)单调递增且取值范围为[0,∞),因此无论y=y1/x1如何取值,都能找到与f(x)唯一的一个交点,这个交点的横坐标就是方程(18)的解.这就证明了解的存在性和唯一性.

图4函数(19)对应函数曲线

2.2使用微分方程求解

由方程(7)可得

C=2r为常数,初值条件y(0)=0.微分方程(20)数值求解存在两个困难:1)y′的符号未定;2)待定常数C未知.由图1可知y′取值为先正后负,是否能取负值,取决于y1/x1的大小.在数值求解时,可实时判断y′取值,当y′=0时,之后的y′取值应都为负值.C作为待定参数,可通过拟合的方法求解,当C取到准确值时有y(x1)=y1.

定义如下函数:

函数F(t)表示将方程(20)中C看作变量(C→t,这里使用t是为了避免和x混淆)并求解微分方程(20)后y(x1)的值,显然Φ(t)=0的根就是待定系数C.这样拟合过程就是一个数值方程的求解问题.此微分方程完整的求解方法如下:

IgorPro函数f1计算式(21)函数Φ(t),BrachistochroneFunc计算方程(20)微分值.在使用求解方程命令FindRoots调用函数f1拟合求解Φ(t)=0时,需要指定根的范围.和2.1证明解的唯一性方法类似,利用数值方法计算t取不同值时Φ(t)的值,通过Φ(t)值域确定根的范围.可使用f1计算Φ(t),方法为ddd=f1(y1,x).结果如图5所示.

图5Φ(t)在t(C→t)取不同值时的分布

从图5可以看出随着y1/x1减小,Φ(t)=0的根逐渐减小.当y1/x1取图5所示值时,根的范围可取[0.2,2].

Φ(t)是一个十分特殊的函数,在Φ(t)里利用给定的初始值C(就是IgorPro函数f1的参数值t)对微分方程(20)求解,利用求得的解来计算Φ(t)函数值.因此Φ(t)=0是一个包含积分过程的复杂方程,这样的方程一般是无法通过解析方法求解的,但利用数值方法就可以突破这个限制.f1在计算Φ(t)时使用了微分方程求解命令IntegrateODE.IgorPro函数BrachistochroneFunc在计算微分值时,原点附近的斜率非常大,由于数值求解算法积分变量间隔有限,原点附近返回的中间解(非最终收敛解)会剧烈震荡,导致中间解计算出来的y′2是负值,如果直接使用此值计算,会出现对负数开平方的错误.此时可将y′人为设置为一个固定常数-1,这样既可以避免负数开方导致求解过程失败,又不会影响结果,IntegrateODE会继续自动调整积分变量间隔直至积分值收敛,此时返回的解是正常的.

当得到Φ(t)=0根的时候,该根对应的微分方程解就是待求最速降速曲线,因此这里的方法同时拟合出了参数C的值和方程(20)的解.结果如图6,图形标记是微分方程求解结果,实线是利用摆线方程进行数值求解结果(理论值).

图6利用微分方程数值求解最速降速曲线结果

2.3按照最速降线定义数值求解

由1.3,根据如下递推关系计算{xi}:

xn=x0+ch∑ni=1υi1−c2υ2i−−−−−−−√xn=x0+ch∑i=1nυi1-c2υi2(22)

递推关系中含有待定系数C,因此和微分方程求解一样,同样需要先得到C的值.当n=N时,有xN=x1,利用这个关系可以求解C.这里仍然通过求解数值方程方法来拟合求解C.定义如下方程

f(t)=t−xN−x0h∑Ni=1υi1−t2υ2i√=0f(t)=t-xΝ-x0h∑i=1Νυi1-t2υi2=0(23)

方程(23)的根就是待求C的值.求出c之后,代入递推关系式(22),依次计算出xi,就可得到最速降线的数值解.完整的方法如下.

wMathy和wMathx存放了解的数值结果.IgorPro函数f2计算方程(23)对应的函数f(t).使用FindRoots求解f(t)=0同样需要知道f(t)=0根的区间范围,和2.2类似,可对y1/x1取不同值时,f(t)取值进行数值模拟,结果如图7所示.

图7方程(24)中f(t)随t的取值

从图7可以看出,随着y1/x1取值不同,f(t)值域发生变化,特别是当y1/x1=2/π时,f(t)最大值理论上是0(图7中对应曲线略微大于0是数值求解的误差),当y1/x1<2/π时,f(t)将恒小于0,这表示f(t)=0无解.故此方法的适用范围是y1/x1>2/π.在给定y1/x1后还需要知道f(t)=0根的区间,从图7可以看出该区间范围,表1列出了y1/x1不同时方程(23)根分布区间[xL,xH]的建议值.

表1方程f(t)=0根分布区间[xL,xH]随y1/x1建议值

求解结果如图8所示,图形标记是最速降线求解的结果,实线是利用摆线方程求解的结果(理论值).

图8利用定义法求解最速降线问题曲线结果

2.4数值求解最速降线和其他类型曲线时间

设曲线的形式为y=y(x),按照1.2,小球沿曲线y下滑需要的时间为

t=∫x101+y′2−−−−−−√dx2gy−−−√t=∫0x11+y′2dx2gy(24)

如果曲线的解析函数已知,式(24)的数值计算方法则非常简单,通过“微元法”积分便可转化为求和.“微元法”将曲线分为很多小段,精度和效率与“微元”大小及所用算法直接相关.IgorPro可直接对解析函数进行数学运算,由于采用了成熟算法和高度优化后的代码,精度和效率一般都远远优于简单使用“微元法”编写的程序.这里用Integrate1D函数计算式(24)的积分值.

取10条不同形状的二次曲线(y=ax2+bx+c,当a=0时为直线)与最速降线(取x1=1,y1=1,利用2.1介绍方法可解得旋轮圆半径r=0.572917)比较.先按照式(24)积分核函数编写二次曲线和最速降线解析表达式的IgorPro函数ParaBolicFun和BchFunc,然后利用Integrate1D计算每条曲线需要的时间t.

参数取值及计算结果如表2所示,曲线形状如图9所示.

图9小球沿不同曲线下滑轨迹,实线表示二次曲线,虚线表示最速降线

表2小球沿不同曲线轨迹下落需要的时间(x0,y0)=(0,0),(x1,y1)=(1,1)单位:m

表2中黑色粗体标出了曲线分别取直线和最速降线时所用的时间.通过观察图9可知,曲线在初始阶段应尽可能“陡峭”,以在短距离内获得最大速度,此后应尽可能“平缓”,从而能以较大速度通过后段.简单计算可知最速降线在原点的斜率为无穷大,与这里的分析是一致的.

3、结束语

利用数值算法求解微分方程(7)时,在靠近原点的地方,根据数值计算原理人为设定有限斜率值,避免了靠近原点的地方中间解剧烈震荡导致求解失败的问题.在求解待定参数时,将参量处理为变量,设计“参量”函数研究参数的分布范围,将微分方程和数值方程结合起来,通过拟合而不是“穷举”的方式去求解.计算小球沿不同曲线下滑的时间时,没有直接采用“微元法”,而是通过IgorPro的数学函数计算积分值,以提高计算精度和效率.本文的这些处理技巧,对于解决微分方程数值求解过程中的奇异性、数值计算过程中的最优化拟合、数值计算中的解析求解等问题,具有一定参考价值.

在教学中,数值计算和模拟方法能够直观、全面、多角度地对物理原理进行分析和展示,可以让学生获得具体而又感性的认识,而不是仅仅停留在抽象的概念上.本文从最速降线问题的数值求解和小球滑行时间的数值计算出发,详细讨论了求解过程中可能遇到的各种问题,例如考虑解是否存在及如何求解等,并尽可能利用图表将求解过程具体化.利用IgorPro的图表绘制和数据分析处理功能,能够使抽象的物理问题形象化,复杂的物理原理简单化,枯燥的物理公式图表化,这对于提高教学效果是很有帮助的,在教学中加强对数据处理工具的使用也是当前物理实验教学的一个基本要求.

参考文献：

[1]金乐,吴炳俊,刘祥树.用Mathematica软件设计最速降线模拟实验的研究[J].大学物理,2015,34(12):27-30.

[2]上海交通大学数学科学学院,寻找最速降线[EB/OL].

[3]王海军.巧解最速降线及其等时性[J].大学物理,2015,34(08):16-18.

[4]邵云.质点沿最速降线和首尾固定的两相连线段下落问题的研究[J].大学物理,2019,38(12):20-27.

[5]刘志勇,边军辉,刘培国.Mathematica在最速降线问题中的应用[J].西安文理学院学报(自然科学版),2013,16(01):91-94.

[6]谢建华.最速降线问题解充分性的证明[J].力学与实践,2009,31(03):82-84.

[7]师玉荣,马君,马丽珍.基于Matlab编程对最速降线问题的研究[J].物理与工程,2012,22(04):11-15.

[8]陈德锋,廖桂颖,王江涌.基于遗传算法的最速降线问题求解[J].力学研究,2015,04(04):76-88.

贾小文,贺秀良,范海英.利用IgorPro对最速降线问题进行数值求解[J].大学物理,2020,39(12):13-19+31.

基金:陆军军事交通学院博士启动经费(41741107)资助.