- 问题重述和分析
- 欧拉方程化简
- 四阶龙格库塔计算数值解
- 用已知的解析解验证数值解
所谓最速降线问题是：设 A 和 B 是铅直平面上不在同一铅直线上的两点，在所有连结 A 和 B 的平面曲线中，求一曲线，使质点在重力作用下，初速度为零时，沿此曲线从A滑行至B的时间最短（忽略摩擦和阻力的影响）；
将A点取为坐标原点，B点取为
(
x
1
,
y
1
)
(x_{1},y_{1})
(x1,y1)，如下图所示：
根据能量守恒定律，质点在曲线
y
(
x
)
y(x)
y(x)上任意取一点处的速度满足：
1
2
m
[
(
d
x
d
t
)
2
+
(
d
y
d
t
)
2
]
=
m
g
y
frac{1}{2}m[(frac{dx}{dt})^{2}+(frac{dy}{dt})^{2}]=mgy
21m[(dtdx)2+(dtdy)2]=mgy可以获得：
d
t
=
d
x
2
+
d
y
2
2
g
y
=
1
+
(
d
y
d
x
)
2
2
g
y
d
x
dt=sqrt{frac{dx^{2}+dy^{2}}{2gy}}=sqrt{frac{1+(frac{dy}{dx})^{2}}{2gy}}dx
dt=2gydx2+dy2
=2gy1+(dxdy)2
dx问题变成求
y
=
y
(
x
)
y=y(x)
y=y(x)，使得泛函最小，同时满足边界条件：
y
(
0
)
=
0
,
y
(
x
1
)
=
y
1
y(0)=0,y(x_{1})=y_{1}
y(0)=0,y(x1)=y1，泛函为：
J
(
y
(
x
)
)
=
∫
0
x
1
1
+
(
d
y
d
x
)
2
2
g
y
d
x
J(y(x))=int_{0}^{x_{1}}sqrt{frac{1+(frac{dy}{dx})^{2}}{2gy}}dx
J(y(x))=∫0x12gy1+(dxdy)2
dx为了解决上述问题，我们取出主部，使用欧拉方程来做，主部即：
1
+
(
d
y
d
x
)
2
2
g
y
=
1
+
y
˙
2
2
g
y
sqrt{frac{1+(frac{dy}{dx})^{2}}{2gy}}=sqrt{frac{1+dot{y}^{2}}{2gy}}
2gy1+(dxdy)2
=2gy1+y˙2
对于符号表达的化简，我们可以使用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 f f，回忆欧拉方程为： ∂ f ∂ y − d d x ( ∂ f ∂ y ˙ ) = 0 frac{partial f}{partial y}-frac{d}{dx}(frac{partial f}{partialdot{y}})=0 ∂y∂f−dxd(∂y˙∂f)=0所以下面，我们获取 ∂ f ∂ y frac{partial f}{partial y} ∂y∂f和 ∂ f ∂ y ˙ frac{partial f}{partialdot{y}} ∂y˙∂f：
# 对y求导 lhs = diff(f, y) # 对y'求导 rhs = diff(f, diff(y))
可以打印出结果：
由于欧拉方程等号右边为零，我可以消去 ∂ f ∂ y frac{partial f}{partial y} ∂y∂f和 ∂ f ∂ y ˙ frac{partial f}{partialdot{y}} ∂y˙∂f的公共因子，得到：
# 考虑到欧拉方程中的形式,提取因子化简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 C C：
eq3 = ((eq2) * (diff(y, x))).expand() eq4 = y + y * diff(y)**2 - C
所以我们得到了化简后的方程： − C + y y ˙ 2 + y = 0 -C+ydot{y}^{2}+y=0 −C+yy˙2+y=0现在获取符号上的解：
solve(eq4, diff(y))
结果为 y ˙ = [ − C y − 1 , C y − 1 ] dot{y}=[-sqrt{frac{C}{y}-1},sqrt{frac{C}{y}-1}] y˙=[−yC−1 ,yC−1 ]由于质点是往下滑行，所以 y ˙ dot{y} y˙应当是一个负数，所以可以确定： y ˙ = − C y − 1 dot{y}=-sqrt{frac{C}{y}-1} y˙=−yC−1 ；
四阶龙格库塔计算数值解当前大部分计算工具（matlab，scipy，sympy，numpy）对于微分方程的api，多数只适用于常微分方程（Ordinary Differential Equation，ODE），面对微分方程 − C + y y ˙ 2 + y = 0 -C+ydot{y}^{2}+y=0 −C+yy˙2+y=0，我们难以找到一个合适的api获得解析解。因此，我将用数值分析的方式去解决此问题，针对上面的表达： y ˙ = − C y − 1 dot{y}=-sqrt{frac{C}{y}-1} y˙=−yC−1 这个形式让我想到了我可以用四阶龙格库塔去计算，下面简要回顾四阶龙格库塔：
在区间
[
x
n
,
x
n
+
1
]
[x_{n},x_{n+1}]
[xn,xn+1]内，取四个不同的点可以构造出四阶Runge-Kutta格式：
其中，
h
=
x
n
+
1
−
x
n
h=x_{n+1}-x_{n}
h=xn+1−xn代表步长，
x
n
+
1
2
=
x
n
+
1
2
h
x_{n+frac{1}{2}}=x_{n}+frac{1}{2}h
xn+21=xn+21h，以及
f
(
x
n
,
y
n
)
f(x_{n},y_{n})
f(xn,yn)代表导数：
f
(
x
n
,
y
n
)
=
f
(
x
n
,
y
(
x
n
)
)
=
y
(
x
n
+
1
)
−
y
(
x
n
)
h
f(x_{n},y_{n})=f(x_{n},y(x_{n}))=frac{y(x_{n+1})-y(x_{n})}{h}
f(xn,yn)=f(xn,y(xn))=hy(xn+1)−y(xn)由于我前面获得了
y
˙
=
−
C
y
−
1
dot{y}=-sqrt{frac{C}{y}-1}
y˙=−yC−1
，所以我其实获得了
f
(
x
n
,
y
n
)
f(x_{n},y_{n})
f(xn,yn)，然后边界条件之一有
y
(
0
)
=
0
y(0)=0
y(0)=0，因此我可以从起点开始逐步计算，从而得到一个带有符号常数
C
C
C的轨迹，这里需要注意：我应该确保有
y
(
0
)
=
−
0.001
y(0)=-0.001
y(0)=−0.001这样一个很小的数，不然分母为零无法算出
y
˙
dot{y}
y˙。
机器携带符号计算是我不希望看到的，这很容易导致程序执行出错，所以我应该生成一个列表，其中保存了大量的常数 C C C的不同数值，我们将每个常数代入龙格库塔中可以获得一个序列（即各条最速线的轨迹），剩下问题是选择哪个常数？注意到我还有一个边界条件：终点 ( x 1 , y 1 ) (x_{1},y_{1}) (x1,y1)，所以我可以这样做：
- 取出每个序列的 x 1 x_{1} x1处对应的 y y y值，检验并选择最接近 y 1 y_{1} y1的那个 y y y，用它对应的线作为最速线，而它对应的常数就是对于我们这个问题的最优常数；
另外需要注意一个问题，我需简单分析一下常数 C C C的范围，这可以提高搜索的效率，观察 y ˙ = − C y − 1 dot{y}=-sqrt{frac{C}{y}-1} y˙=−yC−1 ，根号内的值应该大于0（因为质点是下滑的，导数应该小于0），所以常数 C C C的值应该小于质点下滑的最小纵坐标处，也就是终点处。
我现在把上述过程实现，并且拟定终点的边界条件为 ( 1.6 , − 1.0 ) (1.6,-1.0) (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
用已知的解析解验证数值解
对于终点为
(
1.6
,
−
1.0
)
(1.6,-1.0)
(1.6,−1.0)的情况，已知解析解的参数方程如下：
x
=
1
2
(
t
−
s
i
n
t
)
,
y
=
1
2
(
1
−
c
o
s
t
)
x=frac{1}{2}(t-sint),y=frac{1}{2}(1-cost)
x=21(t−sint),y=21(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
x=1.6
x=1.6
