Python 求解自行车前后轮轨迹问题
作者:crazyhat,Python及科学计算爱好者。
数月前偶遇一道自行车相关的趣味数学题:根据下图[1]所示自行车前、后轮轨迹,判断自行车的前进方向,是从左往右还是从右往左?
自行车趣题
不曾想过,平常无奇的自行车的两轮之间还蕴藏着恒定的规律。在赞叹答案思路清奇之余,也萌生了基于前轮轨迹用Python画出相应后轮轨迹的想法;于是便有了本文。
数学模型
文章开头的问题最早出自《Which Way did the Bicycle Go? … and Other Intriguing Mathematical Mysteries》一书,封面用的即是这道看似无法求解的问题。福尔摩斯探案集系列中也设计了类似的案件,可惜给出的却不是正确的解释。在揭晓答案之前,不妨先花几分钟尝试一下能否得到思路。
封面
这实在是一道精彩的几何问题。考虑一下自行车的结构和特点:前轮可以自由控制行进方向,而后轮因为焊死在车架上而只能沿着车身方向行进。因此,后轮轨迹曲线上任一点的切线方向即为车身方向,并且该方向与前轮轨迹曲线的交点即为自行车前后轮轮毂中心的距离,而这保持一个定长(车身长度)。用数学语言描述为:
作后轮曲线上任一点P
的切线,得到与前轮轨迹曲线的交点Q
,线段PQ
保持定长。
所以解题方法为,分别假设前后轮轨迹,然后按照上述结论进行验证即可。本题答案为自行车从左向右行驶,具体参考下图:
自行车趣题解法
好了,数学趣题到此结束。而本文才刚刚开始——已知任意给定的前轮轨迹,如何求解相应的后轮轨迹?为方便起见,我们用参数方程来描述曲线上任一点的坐标(x,y)
,例如已知前轮轨迹x=α(t), y=β(t)
,求解后轮轨迹x=x(t), y=y(t)
。
由于数学并不是本文的重点,略去具体的推导步骤(高中数学知识即可)得到下图中红色所示的x'
和y'
的方程。这两个方程便是冥冥之中的“规律”,默默地支配着我们习以为常的自行车后轮的运动轨迹。
数学模型
常微分方程数值解
与我们初等数学阶段学习的方程的不同之处在于,这是关于未知数(x, y)
的一阶导数的方程。实际上,大自然正是以微分方程的形式控制着一切的运动和变化,无论是日月星辰的运行,还是手中一杯茶的冷却。所幸的是,本文要处理的问题相对来说非常简单,没有高阶导数、多自变量、非线性的耦合,而仅仅只是一阶常微分方程组
。
然而,即便这么一个“简单”的问题,我们大多数情况下依然不能解析求解后轮曲线的代数方程;所以本文采用的是数值求解的方法。欧拉在求解形如y'(x)=f(x, y)
的微分方程的数值解问题上给出了一些标准解法。本文待求解的方程组中α, β
都是参数t
的函数,因此可以归结为x'(t)=f(t, x, y)
和y'(t)=g(t, x, y)
。其中t
为自变量,x
和y
是待求解的量。
为了套用标准的y'(x)=f(x, y)
,我们将上述方程组改写为向量形式Y'=F(t, Y)
。其中,
Y=(x,y)
表示待求解微分方程组的各个分量组成的向量,也就是后轮曲线上的x
和y
坐标;F=(f,g)
表示待求解微分方程组的各个控制方程。
作为练习,我们当然可以自己实现欧拉方法或者改进的算法(例如Runge-Kutta
方法)来求解上述问题,但这不是本文的重点。所以我们直接采用Python科学计算库scipy
中的odeint()
方法来求解常微分方程。经过漫长的铺垫,终于来到Python部分了!
odeint()
函数声明如下,除了最重要的前三个参数外,其余统一用**kw
表示,具体参考官方文档[2]。
scipy.integrate.odeint(func, y0, t, **kw)
func
表示待求解的函数对象,默认参数列表callable(Y, t, …)
,即待求解变量在前,自变量t
在后(可以设置tfirst=True
改变为callable(t, Y, …)
的形式);如果求解微分方程组,则函数返回值为各个微分方程返回值组成的向量。y0
表示给定的初始条件。t
表示自变量的求解区间或者给定序列;注意t[0]
必须正好对应初始解y0
。
函数返回值为二维numpy
数组,每一行对应给定的自变量序列t
的一个值,每一列对应待求解变量y
的一个分量。
Python求解方案
前两节分别建立了数学模型和给出了求解方法,至此理论上这个问题已经解决了,不过我们的Python实践才刚刚开始。考虑到微分方程的解依赖于初始条件,也就是说即便前轮走同样的轨迹,后轮轨迹也会因为刚启动自行车时车身姿态(后轮位置)的不同而不同。当初想到这个问题,也正是想探究以不同姿态启动时,后轮是如何追随前轮轨迹的。所以,下面以前轮轨迹参数方程为构造函数参数开始自行车轨迹问题的面向对象求解。
import autograd.numpy as np
from autograd import grad
from scipy.integrate import odeintclass BicycleTrack:'''Solving and plot back wheel track according to front wheel track and initial position of the bicycle.The track of front wheel is given by:x = fx(t)y = fy(t)Arguments:fx : function object of front wheel track x componentfy : function object of front wheel track y component''' def __init__(self, fx, fy):# front wheel trackself.front_track_x = fxself.front_track_y = fy# first derivative of front wheel track on parameter tself.dfx = grad(fx)self.dfy = grad(fy)# solved back track represented with t, x, yself.t, self.X, self.Y = None, None, None # back trackself.FX, self.FY = None, None # front track
从这个问题的控制方程可知,我们需要求解前轮轨迹的一阶导数,这里借助第三方库autograd
来实现[3]。例如上述代码中self.dfx = grad(fx)
即为得到前轮横坐标参数方程的函数对象fx
的一阶导函数对象dfx
。此外,为了避免使用该库过程中出现未知的问题,平时导入numpy
库的代码改为import autograd.numpy as np
。
接下来定义待求解的微分方程组:
def governing_equation(self, t, Y):''' ODEs of Bicycle Track Problem '''x, y = Yk1 = np.array([self.dfx(t), self.dfy(t)])k2 = np.array([x-self.front_track_x(t), y-self.front_track_y(t)]) return np.sum(k1*k2) * k2 / self.L**2
本文采用了
callable(t, Y, …)
形式的参数列表(仅仅因为个人习惯),因此后面被odeint
调用时需要指明tfirst=True
;self.L
为自行车车身长度(前后轮毂中心距离),它在给定求解区间和初始条件时被确定。
最后,调用odeint
求解数值解:
def solve(self, span, P0, num=100):''' span: solving range of parameter tP0 : initial position of back wheel (x, y)'''# initial point of back wheelP0 = np.array(P0)# initial point of front wheel is defined by parametric equationst0, t1 = spanQ0 = np.array([self.front_track_x(t0), self.front_track_y(t0)])# frame length is defined by P0 and Q0self.L = np.sum((P0-Q0)**2)**0.5# solvingself.t = np.linspace(t0, t1, num)res = odeint(self.governing_equation, P0, self.t, tfirst=True)# solved back trackself.X, self.Y = res[:, 0], res[:, 1]# front wheel trackself.FX, self.FY = self.front_track_x(self.t), self.front_track_y(self.t)
self.FX
、self.FY
、self.X
和self.Y
分别保存为求解后的前后轮横纵坐标值,可以直接用于后续的可视化。
验证与结果展示
曳物线
[4]为例进行验证。曳物线可以理解为:自行车初始沿着y
轴停放且保持前轮正好在原点,然后前轮沿着x
轴正方向行进,此时后轮的轨迹即为曳物线。
- 前轮轨迹非常简单:
x=t, y=0
- 根据曳物线的解析表达式得到后轮的解析解为:
x=L(t-tanht), y=Lsecht
(L
为车长)。
import numpy as np
import matplotlib.pyplot as plt
from bicycle_track import BicycleTrackL = 1 # length# numeric solution
BT = BicycleTrack(lambda t: t, lambda t: 0*t)
BT.solve([0, 5], np.array([0, L]), 50) # analytical solution
t = BT.t
x = L*(t-np.sinh(t)/np.cosh(t))
y = L/np.cosh(t)# plot
plt.plot(BT.X, BT.Y, 'ro', label='numeric solution')
plt.plot(x, y, 'deepskyblue', label='analytical solution')
plt.legend()
plt.show()
最后,求解一个相对复杂的运动并以动画形式展示结果。限于篇幅,matplotlib
动画过程具体参考文末汇总代码。
轨迹动画
BicycleTrack
类及主要测试代码汇总如下:https://github.com/dothinking/blog/blob/master/src/bicycle_track/bicycle_track.py
[1] [1.9 Which Way Did the Bicycle Go? – and other bicycle ponderings](https://gdaymath.com/lessons/gmp/1-9-which-way-did-the-bicycle-go-and-other-bicycle-ponderings/)
[2][scipy.integrate.odeint](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html)
[3] [HIPS/autograd](https://github.com/HIPS/autograd)
[4] [曳物线](https://www.lfhacks.com/t/tractrix)
赞 赏 作 者
Python中文社区作为一个去中心化的全球技术社区,以成为全球20万Python中文开发者的精神部落为愿景,目前覆盖各大主流媒体和协作平台,与阿里、腾讯、百度、微软、亚马逊、开源中国、CSDN等业界知名公司和技术社区建立了广泛的联系,拥有来自十多个国家和地区数万名登记会员,会员来自以工信部、清华大学、北京大学、北京邮电大学、中国人民银行、中科院、中金、华为、BAT、谷歌、微软等为代表的政府机关、科研单位、金融机构以及海内外知名公司,全平台近20万开发者关注。
▼ 点击成为社区注册会员 「在看」一下,一起PY
Python 求解自行车前后轮轨迹问题相关推荐
- python 文件名随自变量变化_Python 求解自行车前后轮轨迹问题
原标题:Python 求解自行车前后轮轨迹问题 转载自:python中文社区 作者:crazyhat,Python及科学计算爱好者. 数月前偶遇一道自行车相关的趣味数学题:根据下图[1]所示自行车前. ...
- 汽车转向前后轮轨迹matlab程序,车前进后退方向的口诀,动画图解车前后轮转弯轨迹...
我们在学习驾驶前,最好先掌握一些车辆行驶原理,比方说汽车转弯轨迹.前进和后退打方向的关系等,只有具备这些基础后,驾驶才得心应手,下面提供汽车前轮和后轮转向轨迹动画图解. 车前进后退方向的口诀 口诀一: ...
- 自行车为什么前轮和后轮受到的摩擦力相反呢 自行车前轮后轮转动方向一样 自行车运动原理...
后轮是主动轮是引擎(脚链子)给与的驱动使得后龙有向后传的趋势从而产生向前驱动的地面静摩擦(若地面完全光滑则后轮只会打滑不会再前进了). https://zhidao.baidu.com/questio ...
- mysql求回购率_用户行为分析——回购率、复购率(SQL、Python求解)
有一个多月没有用Python了,有些生疏o(╥﹏╥)o.通过秦路老师的一道题目,分别使用sql和python求解,顺便复习下python点,重点关注[复购率].[回购率]的解法 ☞秦路老师视频讲解(使 ...
- python求解迷宫问题,配js实现的走迷宫动画,动起来才有意思~
前言 继昨天手动实现了走迷宫问题,虽然是实现了,但是看到被我画成乱七八糟的草稿纸,总是觉得不爽,不仔细看,又得把自己给走迷糊了,于是自己使用js实现了一下,效果还不错!先看一下展示效果吧!(文末配有j ...
- 偏微分方程数值解法python_基于python求解偏微分方程的有限差分法资料
基于python求解偏微分方程的有限差分法资料 Computer Era No. 11 2016 0 引言 在数学中, 偏微分方程是包含多变量和它们的偏 导数在内的微分方程.偏微分方程通常被用来求解 ...
- Python 求解非零和博弈的纳什均衡策略——以虚构的两个企业之间的商品价格博弈为例
系列文章目录 目录 系列文章目录 文章目录 前言 一.模型的建立 二.算法的步骤 三.代码的实现 四.运行的结果 总结 前言 博弈论是研究具有斗争或竞争性质现象的数学理论和方法,它考虑游戏中的个体的预 ...
- 贝叶斯定理及其python求解
文章目录 1. 什么是著名的贝叶斯定理: **2.在实际问题中利用贝叶斯定理求解后验概率的一般步骤:** 3. 典型例子:曲奇饼问题 1. 什么是著名的贝叶斯定理: 就是这么一个不起眼的公式,却能解决 ...
- python如何求解微分方程_常微分方程数值解:Python求解
这里对使用python求解常微分方程提供两种思路,一种是自己编程实现欧拉法,改进欧拉法或者四阶龙格库塔,这样有助于理解上述三种数值计算方法的原理:一种是调用python已有的库,不再重复造轮子. 本文 ...
最新文章
- 面向对象编程(OOP)----BLUE大师JS课堂笔记(二)
- asyncdata 获取参数_载入页面初始数据(asyncData)《 Nuxt.js:异步数据 》
- 可视化生信分析利器 Galaxy 之 Docker 部署
- Python 3.X 练习集100题 02
- 【树莓派】更新系统镜像下载地址,可能是最简单粗暴的树莓派搭建个人网站教程...
- python tkinter 背景色改变不了_python – 在Tkinter中动态更改小部件背景颜色
- sklearn模型评选择与评估
- 计算机考研考编程,计算机考研面试------编程语言
- GraphQL一些hello world级别的例子
- 2017-12-04HTML table布局
- sql语句的执行过程和优化
- Xcode7.0.1:升级Xcode7上传AppStore失败问题
- 因设备需求超供应预期 摩托罗拉折叠机Razr推迟在美上市时间
- TCP/IP网络协议栈面试经典题目
- Node.js nodemon
- python爬虫面试自我介绍范文_走过路过不容错过,Python爬虫面试总结
- Cocos2d lua 破解方案集合
- JPA结合querydsl使用
- 融云 SDK 5.0.0 功能迭代
- 关于sfc /scannow后主题文件的重置