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 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)

验证与结果展示

• 前轮轨迹非常简单：x=t, y=0
• 根据曳物线的解析表达式得到后轮的解析解为：x=L(t-tanht), y=Lsecht（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()

最后，求解一个相对复杂的运动并以动画形式展示结果。限于篇幅，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