问题重述和分析

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

将A点取为坐标原点，B点取为 $x_{1},y_{1})$ ，如下图所示：
fig1
根据能量守恒定律，质点在曲线 $y (x)$ 上任意取一点处的速度满足： $\frac{1}{2}m[(\frac{dx}{dt})^{2}+(\frac{dy}{dt})^{2}]=mgy$ 可以获得： $dt=\sqrt{\frac{dx^{2}+dy^{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))

可以打印出结果：
fig2

由于欧拉方程等号右边为零，我可以消去 $\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)

fig3
表示出欧拉方程，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格式：
fig4
其中， $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()