Python代码中的数学之美：从自由落体到爬虫悖论，十分钟开启数学思维

数学思维，就是用数学的方式去解决问题，就象吃饭用筷子、喝水用杯子一样，自然而然又理所当然。数学思维并非知识的积累，而是一种由特定思维习惯蕴育而成的能力——这种特定习惯的养成，往往是从解决看似简单的问题开始。本文正是从最简单的自由落体运动抽象出数学分析的一般性原理，进而用来回答蠕虫悖论问题。即使对数学抱有恐惧心理的同学，若能静下心来花上十分钟阅读一遍，也一定可以从中体会到数学之美以及数学分析这个工具的威力。

将长度1米的橡皮筋的一端固定，小明拉着另一端以1米/秒的速度向前奔跑。与此同时，一条毛毛虫从橡皮筋固定的一端开始以1厘米/秒的速度追赶小明。假定橡皮筋可以无限均匀拉伸，小明和毛毛虫都可以长生不老，那么，毛毛虫最终能够追得上小明吗？

这就是被称为“蠕虫悖论”的经典问题，据说是由古希腊数学家芝诺（Zeno of Elea）提出来的——我对此表示怀疑：四千多年前的古希腊人应该不会制造橡皮筋吧？古希腊儿童也许会跳皮筋，但跳的肯定不是橡皮筋。

不过，考古不是我们的话题，搞明白毛毛虫最终能否追得上小明才是当务之急。乍一看，这是一个典型的追击问题，小学生最拿手。可细细琢磨，每秒钟爬1厘米的毛毛虫似乎永远也不可能追上每秒钟跑100厘米的小明。

抱着脑袋想了很久，稍微有了一点眉目：毛毛虫爬过的距离和橡皮筋总长之比kkk，一定是时间ttt的函数，记作：
k=f(t)k = f(t)k=f(t)

因为橡皮筋的拉伸是均匀的，当毛毛虫处于某个位置时，如果停止爬行，kkk就保持不变；只要向前爬，kkk就一定是单向递增的；当kkk等于1时，就意味着毛毛虫追上了小明，此时对应的ttt就是毛毛虫追上小明所花费的时间。

可是，f(t)f(t)f(t)又该如何表示呢？既然是追击问题，我们不妨梳理一下从小到大积累起来的、和速度路程相关的知识点，也许用最简单的思想就能解决这个看似复杂的问题。

小学阶段，我们学会了解决以匀速运动为前提的速度和路程问题。假设物体在ttt时间内移动的距离为sss，则物体移动的平均速度vvv记作：

v=stv = \frac{s}{t}v=ts

中学阶段，我们理解了加速度的概念。假设物体运动的初始速度是v0v_0v0，加速度为aaa，经过ttt时间后的物体运动速度vvv记作：

v=v0+atv = v_0 + atv=v0+at

对于自由落体而言，加速度一般用ggg表示（约等于9.8m/s29.8m/s^29.8m/s2），其初始速度为0，自由落体的瞬时速度vvv记作：

v=gtv = gtv=gt

在ttt时间内，自由落体下落的距离sss记作：

s=12gt2s = \frac{1}{2}gt^2s=21gt2

这里，自由落体的瞬时速度vvv和下落距离sss都是时间ttt的函数，二者之间可以互相推导：对距离函数sss求导就是瞬时速度函数vvv，对瞬时速度函数vvv积分，就是距离函数sss。

v=s′=(12gt2)′=gtv = s' = ( \frac{1}{2}gt^2)' = gtv=s′=(21gt2)′=gt
s=∫vdt=∫gtdt=12gt2+Cs = \int{vdt} = \int{gtdt} = \frac{1}{2}gt^2 + Cs=∫vdt=∫gtdt=21gt2+C

上式中的常数C表示在自由落体开始降落前的下降距离，显然，C=0。如果理解这两个公式有点吃力，那么请牢记以下三点，这是应用数学分析的方法解决实际问题的三大法宝。

导数是变化率
导数是速度
导数是斜率

导数是变化率！忽然间灵光乍现：既然毛毛虫爬过的距离和橡皮筋总长之比kkk函数不容易直接写出来，何不尝试一下这个函数的导数k′k'k′——也就是在任意的ttt时刻单位时间内毛毛虫的爬行距离与橡皮筋的长度之比？

不难理解，任意ttt时刻的橡皮筋长度为1+t1+t1+t；单位时间内毛毛虫的爬行距离是1厘米，换算成标准单位是0.01米。MyGod，没想到k′k'k′就这么轻易地写出来了：

k′=0.011+tk' = \frac{0.01}{1+t}k′=1+t0.01

接下来，只要对k′k'k′积分，就可以得到函数kkk了。积分，就得需要求导公式。也许有些同学已经忘记了如何求导，没关系，只要理解导数的意义，求导公式可以很容易地在网上找到。下面是常用六大基本初等函数的求导公司，其中三角函数和反三角函数的求导公式，我只记住了两个。

常数函数
(C)′=0(C)' = 0(C)′=0
幂函数
(xn)′=nxn−1(x^n)' = nx^{n-1}(xn)′=nxn−1
指数函数
(ax)′=axln⁡a(a^x)' = a^x\ln{a}(ax)′=axlna
(ex)′=ex(e^x)' = e^x(ex)′=ex
对数函数
(ln⁡x)′=1x(\ln{x})' = \frac{1}{x}(lnx)′=x1
三角函数
(sin⁡x)′=cos⁡x(\sin{x})' = \cos{x}(sinx)′=cosx
(cos⁡x)′=−sin⁡x(\cos{x})' = -\sin{x}(cosx)′=−sinx
反三角函数
(arcsin⁡x)′=11−x2(\arcsin{x})' = \frac{1}{\sqrt{1-x^2}}(arcsinx)′=1−x21
(arccos⁡x)′=−11−x2(\arccos{x})' = -\frac{1}{\sqrt{1-x^2}}(arccosx)′=−1−x21

显而易见，k′k'k′的模样最接近对数函数的导数形式，可以很轻松地写出积分式：

k=∫k′dt=∫0.011+tdt=0.01×ln⁡(1+t)+Ck = \int{k'dt} = \int{\frac{0.01}{1+t}dt} = 0.01\times\ln{(1+t)} + Ck=∫k′dt=∫1+t0.01dt=0.01×ln(1+t)+C

在毛毛虫开始爬行前kkk为0，因此上式中的CCC也等于0。有了kkk函数，蠕虫悖论问题就变成了求解如下方程的问题：

0.01×ln⁡(1+t)−1=00.01\times\ln{(1+t)} - 1 = 00.01×ln(1+t)−1=0

至此，蠕虫悖论问题总算走对了方向，剩下的工作都是水到渠成、手到擒来的事了。对于这个方程，至少有3种方式可以求解。

1. 调包侠：使用SciPy模块求解非线性方程

SciPy的优化与求根子模块scipy.optimize提供了非线性方程的求解方案。使用SciPy求解非线性方程要求将方程表示为 f(x)=0f(x)=0f(x)=0 或f(x,y,…)=0f(x,y,…)=0f(x,y,…)=0的形式，这里的函数fff就是求根函数的参数。如果函数fff只有一个未知数，则称为标量函数，可用optimize.root_scalar()函数求根；如果函数fff包含多个未知数，则称为向量函数，可用optimize.root()函数求根。

对标量函数求根的函数optimize.root_scalar()，除了需要传入标量函数这个必要参数，还接受method（求根方法）、bracket（求根区间）、fprime（标量函数是否可导）等参数。大多数情况下，只需要给出 bracket 参数即可。对于本例，猜测方程的根在10410^4104到105010^{50}1050之间。

>>> import numpy as np
>>> from scipy import optimize
>>> def k(t):return 0.01*np.log(1+t) - 1>>> result = optimize.root_scalar(k, bracket=[1e4, 1e50]) # bracket指定求根区间
>>> result.root # 方程的根
2.688117141816133e+43
>>> k(result.root) # 检验方程根，越接近0，近似程度越高
0.0
>>> result.iterations # 迭代求解的次数
32

结果显示，经过32次迭代，root_scalar()返回的方程根为2.688×10432.688\times10^{43}2.688×1043，带入原方程检测，误差为0。也就是说，毛毛虫最终一定可以追上小明，只是时间有点长，大约需要2.688×10432.688\times10^{43}2.688×1043秒钟，相当于8.524×10358.524\times10^{35}8.524×1035年。

2. 经典派：使用牛顿逼近法求解非线性方程

如果不喜欢调包，这个方程也可以用我在《Python代码中的数学之美：用牛顿逼近法计算2的算术平方根》一文中介绍的牛顿逼近法来求解。

>>> import numpy as np
>>> def k(t):return 0.01*np.log(1+t) - 1>>> def dk(t):return 0.01/(1+t)>>> def newton_raphson_method():# 牛顿-拉弗森方法tiny = 1e-15 # 当k(t_n)小于tiny时，迭代结束i, t = 0, 10000 # 迭代计数器和初始值while True:i += 1 # 迭代计数器加1t = t-k(t)/dk(t) # 计算下一个tif abs(k(t)) < tiny: # 满足迭代结束条件return(i, t)>>> newton_raphson_method()
(32, 2.688117141816148e+43)

如果以10410^4104作为初始值的话，同样经过32次迭代，使用牛顿逼近法求得的方程根为2.688×10432.688\times10^{43}2.688×1043，精度和上面的方法相比略有差异，小到可以忽略不计。

3. 太极派：四两拨千斤

不管是调包侠还是经典派，本质上都是反复迭代渐次逼近的方法，所得到的方程根是近似解。实际上，这个方程可以用纯数学的方式轻松得到精确解。以下是求解过程

0.01×ln⁡(1+t)−1=00.01\times\ln{(1+t)} - 1 = 00.01×ln(1+t)−1=0
⇒ln⁡(1+t)=100\Rightarrow\ln{(1+t)} = 100⇒ln(1+t)=100
⇒eln⁡(1+t)=e100\Rightarrow e^{\ln{(1+t)}} = e^{100}⇒eln(1+t)=e100
⇒1+t=e100\Rightarrow 1+t = e^{100}⇒1+t=e100
⇒t=e100−1\Rightarrow t = e^{100} - 1⇒t=e100−1

最终求得t=e100−1t=e^{100}-1t=e100−1，这个数字写出来如下：

>>> import numpy as np
>>> np.power(np.e, 100) - 1
2.6881171418161212e+43

这个精确解和前两种方式计算出来的近似解，尾数的前13位相同，三者之间的差异几乎可以忽略。