Python 求解自行车前后轮轨迹问题 - 文章 - 开发者社区

picture.image

♚

作者：crazyhat，Python及科学计算爱好者。

数月前偶遇一道自行车相关的趣味数学题：根据下图[1]所示自行车前、后轮轨迹，判断自行车的前进方向，是从左往右还是从右往左？

picture.image

自行车趣题

不曾想过，平常无奇的自行车的两轮之间还蕴藏着恒定的规律。在赞叹答案思路清奇之余，也萌生了基于前轮轨迹用Python画出相应后轮轨迹的想法；于是便有了本文。

数学模型

文章开头的问题最早出自《Which Way did the Bicycle Go? … and Other Intriguing Mathematical Mysteries》一书，封面用的即是这道看似无法求解的问题。福尔摩斯探案集系列中也设计了类似的案件，可惜给出的却不是正确的解释。在揭晓答案之前，不妨先花几分钟尝试一下能否得到思路。

picture.image

封面

这实在是一道精彩的几何问题。考虑一下自行车的结构和特点：前轮可以自由控制行进方向，而后轮因为焊死在车架上而只能沿着车身方向行进。因此，后轮轨迹曲线上任一点的切线方向即为车身方向，并且该方向与前轮轨迹曲线的交点即为自行车前后轮轮毂中心的距离，而这保持一个定长（车身长度）。用数学语言描述为：

作后轮曲线上任一点 P的切线，得到与前轮轨迹曲线的交点Q，线段PQ保持定长 。

所以解题方法为，分别假设前后轮轨迹，然后按照上述结论进行验证即可。本题答案为自行车从左向右行驶，具体参考下图：

picture.image

自行车趣题解法

picture.image

自行车趣题解法

好了，数学趣题到此结束。而本文才刚刚开始——已知任意给定的前轮轨迹，如何求解相应的后轮轨迹？为方便起见，我们用参数方程来描述曲线上任一点的坐标(x,y)，例如已知前轮轨迹x=α(t), y=β(t)，求解后轮轨迹x=x(t), y=y(t)。

由于数学并不是本文的重点，略去具体的推导步骤（高中数学知识即可）得到下图中红色所示的x'和y'的方程。这两个方程便是冥冥之中的“规律”，默默地支配着我们习以为常的自行车后轮的运动轨迹。

picture.image

数学模型

常微分方程数值解

与我们初等数学阶段学习的方程的不同之处在于，这是关于未知数 (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 odeint  
  
  
class 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 component  
            fy  : function object of front wheel track y component  
    '''   
  
    def \_\_init\_\_(self, fx, fy):  
        # front wheel track  
        self.front_track_x = fx  
        self.front_track_y = fy  
  
        # first derivative of front wheel track on parameter t  
        self.dfx = grad(fx)  
        self.dfy = grad(fy)  
  
        # solved back track represented with t, x, y  
        self.t, self.X, self.Y = None, None, None  # back track  
        self.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 = Y  
    k1 = 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 t  
        P0  : initial position of back wheel (x, y)  
    '''  
  
    # initial point of back wheel  
    P0 = np.array(P0)  
  
    # initial point of front wheel is defined by parametric equations  
    t0, t1 = span  
    Q0 = np.array([self.front_track_x(t0), self.front_track_y(t0)])  
  
    # frame length is defined by P0 and Q0  
    self.L = np.sum((P0-Q0)**2)**0.5  
  
    # solving  
    self.t = np.linspace(t0, t1, num)  
    res = odeint(self.governing_equation, P0, self.t, tfirst=True)  
  
    # solved back track  
    self.X, self.Y = res[:, 0], res[:, 1]  
  
    # front wheel track  
    self.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=L*secht （ L 为车长）。


        
import numpy as np  
import matplotlib.pyplot as plt  
from bicycle_track import BicycleTrack  
  
L = 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()

picture.image

验证曳物线

最后，求解一个相对复杂的运动并以动画形式展示结果。限于篇幅，matplotlib动画过程具体参考文末汇总代码。

picture.image

轨迹动画

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

[2]scipy.integrate.odeint

[3] HIPS/autograd

[4] 曳物线

赞赏作者

picture.image

Python中文社区作为一个去中心化的全球技术社区，以成为全球20万Python中文开发者的精神部落为愿景，目前覆盖各大主流媒体和协作平台，与阿里、腾讯、百度、微软、亚马逊、开源中国、CSDN等业界知名公司和技术社区建立了广泛的联系，拥有来自十多个国家和地区数万名登记会员，会员来自以工信部、清华大学、北京大学、北京邮电大学、中国人民银行、中科院、中金、华为、BAT、谷歌、微软等为代表的政府机关、科研单位、金融机构以及海内外知名公司，全平台近20万开发者关注。

picture.image

▼ 点击成为 社区注册会员 「在看」 一下，一起PY