问题重述和分析
所谓最速降线问题是:设 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)2dx问题变成求
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)2dx为了解决上述问题,我们取出主部,使用欧拉方程来做,主部即:
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}{\partial\dot{y}})=0 ∂y∂f−dxd(∂y˙∂f)=0所以下面,我们获取 ∂ f ∂ y \frac{\partial f}{\partial y} ∂y∂f和 ∂ f ∂ y ˙ \frac{\partial f}{\partial\dot{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}{\partial\dot{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+y\dot{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+y\dot{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<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
)
(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