问题重述和分析

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

将A点取为坐标原点，B点取为(x1,y1)(x_{1},y_{1})(x1,y1)，如下图所示：

根据能量守恒定律，质点在曲线y(x)y(x)y(x)上任意取一点处的速度满足：12m[(dxdt)2+(dydt)2]=mgy\frac{1}{2}m[(\frac{dx}{dt})^{2}+(\frac{dy}{dt})^{2}]=mgy21m[(dtdx)2+(dtdy)2]=mgy可以获得：dt=dx2+dy22gy=1+(dydx)22gydxdt=\sqrt{\frac{dx^{2}+dy^{2}}{2gy}}=\sqrt{\frac{1+(\frac{dy}{dx})^{2}}{2gy}}dxdt=2gydx2+dy2=2gy1+(dxdy)2dx问题变成求y=y(x)y=y(x)y=y(x)，使得泛函最小，同时满足边界条件：y(0)=0,y(x1)=y1y(0)=0,y(x_{1})=y_{1}y(0)=0,y(x1)=y1，泛函为：J(y(x))=∫0x11+(dydx)22gydxJ(y(x))=\int_{0}^{x_{1}}\sqrt{\frac{1+(\frac{dy}{dx})^{2}}{2gy}}dxJ(y(x))=∫0x12gy1+(dxdy)2dx为了解决上述问题，我们取出主部，使用欧拉方程来做，主部即：1+(dydx)22gy=1+y˙22gy\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 pltmechanics.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)

对于主部fff，回忆欧拉方程为：∂f∂y−ddx(∂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()

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

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

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

solve(eq4, diff(y))

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

四阶龙格库塔计算数值解

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

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

其中，h=xn+1−xnh=x_{n+1}-x_{n}h=xn+1−xn代表步长，xn+12=xn+12hx_{n+\frac{1}{2}}=x_{n}+\frac{1}{2}hxn+21=xn+21h，以及f(xn,yn)f(x_{n},y_{n})f(xn,yn)代表导数：f(xn,yn)=f(xn,y(xn))=y(xn+1)−y(xn)hf(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˙=−Cy−1\dot{y}=-\sqrt{\frac{C}{y}-1}y˙=−yC−1，所以我其实获得了f(xn,yn)f(x_{n},y_{n})f(xn,yn)，然后边界条件之一有y(0)=0y(0)=0y(0)=0，因此我可以从起点开始逐步计算，从而得到一个带有符号常数CCC的轨迹，这里需要注意：我应该确保有y(0)=−0.001y(0)=-0.001y(0)=−0.001这样一个很小的数，不然分母为零无法算出y˙\dot{y}y˙。

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

取出每个序列的x1x_{1}x1处对应的yyy值，检验并选择最接近y1y_{1}y1的那个yyy，用它对应的线作为最速线，而它对应的常数就是对于我们这个问题的最优常数；

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

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

xdst,ydst=1.6,-1.0def dyexpr(y,const):# y'的表达,const为常数return -((const/y)-1.0)**0.5def rungekutta(C):"""四阶龙格库塔C为常数,且C<y_min"""x = 0.0y = -0.001h = 0.001xs = [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 += hxs.append(round(x,3))y += h * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0ys.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.0res_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=12(t−sint),y=12(1−cost)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.6x=1.6x=1.6

Python求解最速降线问题相关推荐

python 物理学中的应用_利用python求解物理学中的双弹簧质能系统详解
前言本文主要给大家介绍了关于利用python求解物理学中双弹簧质能系统的相关内容,分享出来供大家参考学习,下面话不多说了,来一起看看详细的介绍吧. 物理的模型如下: 在这个系统里有两个物体,它们的质 ...
python代码物理_利用python求解物理学中的双弹簧质能系统详解
前言本文主要给大家介绍了关于利用python求解物理学中双弹簧质能系统的相关内容,分享出来供大家参考学习,下面话不多说了,来一起看看详细的介绍吧. 物理的模型如下: 在这个系统里有两个物体,它们的质 ...
最大公约数python语言算法_使用Python求解最大公约数的实现方法
这篇文章主要介绍了使用Python求解最大公约数的实现方法,包括用Python表示欧几里得算法和Stein算法的求解原理. 1. 欧几里德算法欧几里德算法又称辗转相除法, 用于计算两个整数a, b的 ...
利用python求解节点介数和边介数
利用python求解节点介数和边介数利用networkx里面的函数betweenness_centrality(G)来求解节点介数和函数edge_betweenness_centrality(G)来 ...
python最大公约数计算_使用Python求解最大公约数的实现方法
1. 欧几里德算法欧几里德算法又称辗转相除法, 用于计算两个整数a, b的最大公约数.其计算原理依赖于下面的定理: 定理: gcd(a, b) = gcd(b, a mod b) 证明: a可以表示 ...
mysql求回购率_用户行为分析——回购率、复购率（SQL、Python求解）
有一个多月没有用Python了,有些生疏o(╥﹏╥)o.通过秦路老师的一道题目,分别使用sql和python求解,顺便复习下python点,重点关注[复购率].[回购率]的解法 ☞秦路老师视频讲解(使 ...
用Python求解线性规划问题
线性规划简介及数学模型表示线性规划简介一个典型的线性规划问题线性规划模型的三要素线性规划模型的数学表示图解法和单纯形法图解法单纯形法使用python求解简单线性规划模型编程思路求解案例例1:使用sci ...
python求解迷宫问题，配js实现的走迷宫动画，动起来才有意思~
前言继昨天手动实现了走迷宫问题,虽然是实现了,但是看到被我画成乱七八糟的草稿纸,总是觉得不爽,不仔细看,又得把自己给走迷糊了,于是自己使用js实现了一下,效果还不错!先看一下展示效果吧!(文末配有j ...
人工智能 --- Python求解线性和非线性规划问题
基于jupyter notebook的Python编程一.线性规划问题求解 1.Excel中大M法与Excel的"规划求解"包对实际问题的求解比较实际例题: 求解以下约束条件的 ...

Python求解最速降线问题

目录

问题重述和分析

欧拉方程化简

四阶龙格库塔计算数值解

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

Python求解最速降线问题相关推荐

最新文章

热门文章