Python求解最速降线问题

问题重述和分析

所谓最速降线问题是：设 A 和 B 是铅直平面上不在同一铅直线上的两点，在所有连结 A 和 B 的平面曲线中，求一曲线，使质点在重力作用下，初速度为零时，沿此曲线从A滑行至B的时间最短（忽略摩擦和阻力的影响）；

将A点取为坐标原点，B点取为$(x_1,y_1)$，如下图所示：

根据能量守恒定律，质点在曲线$y(x)$上任意取一点处的速度满足：$$\frac{1}{2}m[(\frac{dx}{dt}{)}^2+(\frac{dy}{dt}{)}^2]=mgy$$可以获得：$$dt=\sqrt{\frac{d{x}^2+d{y}^2}{2gy}}=\sqrt{\frac{1+(\frac{dy}{dx}{)}^2}{2gy}}dx$$问题变成求$y=y(x)$，使得泛函最小，同时满足边界条件：$y(0)=0,y(x_1)=y_1$，泛函为：$$J(y(x))={\int }_0^{x_1}\sqrt{\frac{1+(\frac{dy}{dx}{)}^2}{2gy}}dx$$为了解决上述问题，我们取出主部，使用欧拉方程来做，主部即：$$\sqrt{\frac{1+(\frac{dy}{dx}{)}^2}{2gy}}=\sqrt{\frac{1+{\dot{y}}^2}{2gy}}$$

欧拉方程化简

对于符号表达的化简，我们可以使用SymPy，安装命令（conda install sympy），下面导入相关包并做一些输出设置（在jupyter notebook中打印latex公式）：

from sympy import *
from sympy.physics import mechanics
import numpy as np
import matplotlib.pyplot as plt

mechanics.mechanics_printing()

定义符号，并引入主部：

# 符号定义
x = Symbol('x')
y = Function('y')(x) # y是x的函数
g = Symbol('g')
C = Symbol('C')

# 主部, diff表示求导
f = sqrt(1 + diff(y)**2) / sqrt(2 * g * y)

对于主部$f$，回忆欧拉方程为：$$\frac{\partial f}{\partial y}-\frac{d}{dx}(\frac{\partial f}{\partial \dot{y}})=0$$所以下面，我们获取$\frac{\partial f}{\partial y}$和$\frac{\partial f}{\partial \dot{y}}$：

# 对y求导
lhs = diff(f, y)
# 对y'求导
rhs = diff(f, diff(y))

可以打印出结果：

由于欧拉方程等号右边为零，我可以消去$\frac{\partial f}{\partial y}$和$\frac{\partial f}{\partial \dot{y}}$的公共因子，得到：

# 考虑到欧拉方程中的形式,提取因子化简lhs和rhs
lhs, rhs = lhs.subs(sqrt(g), 1) * sqrt(2), rhs.subs(sqrt(g), 1) * sqrt(2)

表示出欧拉方程，simplify()用于化简符号表达到分子分母形式：

# 欧拉方程
eq1 = (diff(rhs, x) - lhs).simplify()

进一步化简，只提取 “分子项=0”，很明显这是和上面方程等价的，得到新方程的表达式为：

# 获取 分子=0 这个方程
subs_cons2 = 2 * y**(3/2) * (diff(y, x)**2 + 1)**(3/2)
eq2 = (eq1 * subs_cons2).expand()

对方程降阶，代价是引入常数$C$：

eq3 = ((eq2) * (diff(y, x))).expand()
eq4 = y + y * diff(y)**2 - C

所以我们得到了化简后的方程：$$-C+y{\dot{y}}^2+y=0$$现在获取符号上的解：

solve(eq4, diff(y))

结果为$$\dot{y}=[-\sqrt{\frac{C}{y}-1},\sqrt{\frac{C}{y}-1}]$$由于质点是往下滑行，所以$\dot{y}$应当是一个负数，所以可以确定：$\dot{y}=-\sqrt{\frac{C}{y}-1}$；

四阶龙格库塔计算数值解

当前大部分计算工具（matlab，scipy，sympy，numpy）对于微分方程的api，多数只适用于常微分方程（Ordinary Differential Equation，ODE），面对微分方程$-C+y{\dot{y}}^2+y=0$，我们难以找到一个合适的api获得解析解。因此，我将用数值分析的方式去解决此问题，针对上面的表达：$$\dot{y}=-\sqrt{\frac{C}{y}-1}$$这个形式让我想到了我可以用四阶龙格库塔去计算，下面简要回顾四阶龙格库塔：

在区间$[x_n,x_{n+1}]$内，取四个不同的点可以构造出四阶Runge-Kutta格式：

其中，$h=x_{n+1}-x_n$代表步长，$x_{n+\frac{1}{2}}=x_n+\frac{1}{2}h$，以及$f(x_n,y_n)$代表导数：$$f(x_n,y_n)=f(x_n,y(x_n))=\frac{y(x_{n+1})-y(x_n)}{h}$$由于我前面获得了$\dot{y}=-\sqrt{\frac{C}{y}-1}$，所以我其实获得了$f(x_n,y_n)$，然后边界条件之一有$y(0)=0$，因此我可以从起点开始逐步计算，从而得到一个带有符号常数$C$的轨迹，这里需要注意：我应该确保有$y(0)=-0.001$这样一个很小的数，不然分母为零无法算出$\dot{y}$。

机器携带符号计算是我不希望看到的，这很容易导致程序执行出错，所以我应该生成一个列表，其中保存了大量的常数$C$的不同数值，我们将每个常数代入龙格库塔中可以获得一个序列（即各条最速线的轨迹），剩下问题是选择哪个常数？注意到我还有一个边界条件：终点$(x_1,y_1)$，所以我可以这样做：

• 取出每个序列的$x_1$处对应的$y$值，检验并选择最接近$y_1$的那个$y$，用它对应的线作为最速线，而它对应的常数就是对于我们这个问题的最优常数；

另外需要注意一个问题，我需简单分析一下常数$C$的范围，这可以提高搜索的效率，观察$\dot{y}=-\sqrt{\frac{C}{y}-1}$，根号内的值应该大于0（因为质点是下滑的，导数应该小于0），所以常数$C$的值应该小于质点下滑的最小纵坐标处，也就是终点处。

我现在把上述过程实现，并且拟定终点的边界条件为$(1.6,-1.0)$：

xdst,ydst=1.6,-1.0

def dyexpr(y,const):
    # y'的表达,const为常数
    return -((const/y)-1.0)**0.5

def rungekutta(C):
	"""
    四阶龙格库塔
    C为常数,且C<y_min
    """
    x = 0.0
    y = -0.001
    h = 0.001
    xs = [x]
    ys = [y]
    for i in range(2000):
        k1 = dyexpr(y, C)
        k2 = dyexpr(y + 0.5 * h * k1, C)
        k3 = dyexpr(y + 0.5 * h * k2, C)
        k4 = dyexpr(y + h * k3, C)
        x += h
        xs.append(round(x,3))
        y += h * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0
        ys.append(y)
    return xs,ys

# 用于尝试数值解的常数集合
C_list=[round(-1*cons,3) for cons in np.arange(1.0,2.0,0.01)]

# 每个常数都获得一组数值解
res_list=[]
for const in C_list:
    xs,ys=rungekutta(const)
    index=xs.index(xdst)
    temp=ys[index]
    if np.isnan(temp):
        temp=100.0
    res_list.append((temp-ydst)**2)

# 选出最优的常数
print(res_list)
final_const=C_list[res_list.index(min(res_list))]
print(final_const) # -1.02

用已知的解析解验证数值解

对于终点为$(1.6,-1.0)$的情况，已知解析解的参数方程如下：$$x=\frac{1}{2}(t-sint),y=\frac{1}{2}(1-cost)$$
我可以将它们可视化，从而对比两个解的情况：

# 绘制数值解的结果
xs,ys=rungekutta(final_const) # ys就是我们要的最速降线
print(len(ys))
plt.plot(xs,ys,label='Numerical Solution')

# 根据解析解公式,验证数值解
tlist=np.arange(0,1.2*np.pi,0.001)
x_analytic=[0.5*(t-np.sin(t)) for t in tlist]
y_analytic=[-0.5*(1-np.cos(t)) for t in tlist]
plt.plot(x_analytic,y_analytic,label='Analytic Solution')

# 标记终点
plt.axvline(1.6,c='r')

plt.legend()
plt.show()

可以看出两个解在最开始比较接近，到后面随着误差累加，使数值解与解析解出现偏离，红色的竖线用于标记终点的横坐标$x=1.6$