与咖啡杯有关的数学:焦散线方程

周末在家看了一天 YouTube 上 Numberphile 频道里面 Tadashi Tokieda 教授的视频,非常过瘾。Tokeida 教授的讲述非常生动、幽默,让你从一开始就被深深吸引。

其中有两个视频使用了杯子作为道具,第一个是 “听声辨形”,即通过倾听敲击杯子的音调变化来判断把手的位置:

第二个演示了杯子中小球的旋转方向会根据球的数目发生某种“相变”的现象:

其实与杯子有关的有趣的数学问题还有不少,还有一个可能是大家都已经非常熟悉的,即光线在杯子内壁反射形成的曲线:

从上图中可见光线在杯子内壁反射后形成了一个光斑 (术语叫焦散),焦散的形状是一条优美的曲线,叫做心脏线 (cardioid)。焦散与所有的反射光线是相切的。

形成焦散的原因是,从光源出发的光线在杯子内壁反射以后路径出现了交叉,其中某些位置经过的光线特别密集,所以亮度就更高。

焦散的具体形状其实和光源的位置有关:当光源是一个点光源,且位置恰好在杯子边缘上某一点时得到的才是严格意义下的心脏线;当光源是平行光源时 (可以视作无穷远处的点光源),得到的焦散是肾形线。心脏线和肾形线的形状非常像。当光源距离杯子有一段距离时,焦散的形状介于心脏线和肾形线之间。

于是趁着周末我动手推导了一番焦散的参数方程。这个推导其实网上的资料已经很多了,然而我很懒,看到那些微积分推导就头疼,也记不住它们,所以采用了一个另类点的方法,以下是具体步骤。


\(t\in[a,b]\)\(C_1(t)=(x_1(t),y_1(t))\)\(C_2(t)=(x_2(t), y_2(t))\) 是两条连续可微的曲线。对任何 \(t\),将两条曲线上对应的两点 \(C_1(t)\)\(C_2(t)\) 用线段 \(L\) 连起来,并让 \(t\) 变动,就得到了一组变化的线段 \(L(t)\)

如果存在另外的一条曲线 \(E(t)\),使得 \(E(t)\) 与线段族 \(L(t)\) 中的每一条都相切,就称 \(E(t)\)\(L(t)\) 的包线 (envelope),也叫做焦散线 (caustics)。

我们其实不难推导 \(E(t)\) 的参数方程:假设 \(E(t)\) 存在,固定一个参数 \(t_0\),则当 \(t\to t_0\) 时,切线 \(L(t)\) 会逐渐与 \(L(t_0)\) 重合,它俩的交点会趋近于 \(E(t)\) 上的一点 \(E(t_0)\)。所以记 \(P(t)\)\(L(t)\)\(L(t_0)\) 的交点,则 \[\lim_{t\to t_0}P(t) = E(t_0).\] 两个线段的交点 \(P(t)\) 的坐标可以解线性方程组得到\[P(t)=\begin{align*}\bigg(&\frac{(x_1(t_0)y_2(t_0)-y_1(t_0)x_2(t_0))(x_1(t)-x_2(t))-(x_1(t_0)-x_2(t_0))(x_1(t)y_2(t)-x_2(t)y_1(t))}{(x_1(t_0)-x_2(t_0))(y_1(t)-y_2(t))-(y_1(t_0)-y_2(t_0))(x_1(t)-x_2(t))},\\ &\frac{(x_1(t_0)y_2(t_0)-y_1(t_0)x_2(t_0))(y_1(t)-y_2(t))-(y_1(t_0)-y_2(t_0))(x_1(t)y_2(t)-x_2(t)y_1(t))}{(x_1(t_0)-x_2(t_0))(y_1(t)-y_2(t))-(y_1(t_0)-y_2(t_0))(x_1(t)-x_2(t))}\bigg)\end{align*}.\]

这是个 \(\frac{0}{0}\) 型的极限,使用洛必达法则可以目视写出答案来:

\[E(t_0)=\begin{align*} \bigg(&\frac{ (y_1-y_2)(x_1x_2' - x_1'x_2)-(x_1-x_2)(x_1y_2'-x_2y_1')}{(x_1-x_2)(y_1'-y_2') - (y_1-y_2)(x_1'-x_2')},\\ &\frac{(y_1-y_2)(y_1x_2' - x_1'y_2)-(x_1-x_2)(y_1y_2'-y_1'y_2)}{(x_1-x_2)(y_1'-y_2') - (y_1-y_2)(x_1'-x_2')}\bigg) \end{align*}.\]

这就是 \(E(t)\) 的参数表示。

对于杯子中光线所形成的焦散线方程,假设用单位圆周表示杯子,点光源位于 \(A\),不难看出红色和蓝色两端弧长度相等,所以反射光线是 \(C_1=(\cos t,\sin t)\)\(C_2=(\cos2t,\sin2t)\) 对应点的连线,其中 \(t\) 是从 \(A\) 出发顺时针转过的角度:

于是焦散线的参数方程可以通过将 \(C_1,C_2\) 代入上面 \(E(t)\) 的表达式化简得到。我写了个 Python 小脚本用 sympy 来算:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from sympy import *

t = symbols('t')
x1, y1 = cos(t), sin(t)
x2, y2 = cos(2 * t), sin(2 * t)
dx1, dy1 = diff(x1, t), diff(y1, t)
dx2, dy2 = diff(x2, t), diff(y2, t)
denumerator = simplify((x1 - x2) * (dy1 - dy2) - (y1 - y2) * (dx1 - dx2))
nume1 = (y1 - y2) * (x1 * dx2 - x2 * dx1) - (x1 - x2) * (x1 * dy2 - dy1 * x2)
nume2 = (y1 - y2) * (y1 * dx2 - y2 * dx1) - (x1 - x2) * (y1 * dy2 - dy1 * y2)
nume1, nume2 = simplify(nume1), simplify(nume2)
X = simplify(nume1 / denumerator)
Y = simplify(nume2 / denumerator)
print(X)
print(Y)

sympy 给出的结果是

1
2
(-3*cos(t) + cos(3*t) + 2)/(6*(cos(t) - 1))
-2*sin(t)**3/(3*cos(t) - 3)

这个表达式还不够简洁。我把它代入 Wolfram 的在线计算器里面进一步化简,就得到了 \[E(t) = \left(\frac{\cos2t + 2\cos t}{3}, \frac{\sin2t + 2\sin t}{3}\right).\] 这便是心脏线的参数方程。

 | 

当前网速较慢或者你使用的浏览器不支持博客特定功能,请尝试刷新或换用Chrome、Firefox等现代浏览器