思考:两个椭圆片能粘合成一个立体吗?
By 苏剑林 | 2019-07-21 | 59074位读者 |前两周又在群里看到一个颇为有趣的问题:两个同样大小的椭圆片可以沿着它们的长轴弯曲,沿着边缘线粘贴,能完美地贴合成一个封闭立体吗?问题来源于知乎《两个椭圆片可否以柱面弯曲边缘完美贴合?》。
问题可以用只言片语表达清楚,甚至普通读者都能理解,而问题本身是有一定难度的,这就符合了一个漂亮的问题的条件,所以也就吸引了笔者陆陆续续思考了好多天,最终在昨天算是给出了这类问题通用的列方程思路和数值求解方案,而今天则完成了理论证明,确认两个相同椭圆片总是可以完美贴合。
弧长参数方程 #
作为准备,我们先求出椭圆的弧长参数方程,如下图所示,当从$(a,0)$出发的一段椭圆弧长为$s$时,对应的$y=h(s)$是多少?
理论上,这个问题不难解决。因为椭圆的标准方程为:
\begin{equation}\frac{x^2}{a^2}+\frac{y^2}{b^2}=1\end{equation}
也就是$x=a\sqrt{1-y^2/b^2}$,这样椭圆弧长就可以表示为:
\begin{equation}s=\int_0^y \sqrt{1+\left(\frac{dx}{dy}\right)^2}dy=\int_0^y \sqrt{\frac{a^2 y^2}{b^4-b^2 y^2}+1}dy\label{eq:bs}\end{equation}
这个积分可以用“第二类不完全椭圆积分”来表示。在$[0,b]$区间内,$s$关于$b$是单调递增的,所以反函数必然存在。作为数值计算来说,我们不需要知道这个反函数是多少,我们只需要通过数值积分$\eqref{eq:bs}$算出一堆$(y,s)$的对应关系,将$y,s$交换位置,就相当于得到了关系$y=h(s)$(通过逐点法给出一个函数)。
粘贴边缘方程 #
现在来看原始的问题:假设如果能粘合的情况下,边缘曲线方程是啥?我们按照下图(左图)所示来建立坐标系:
在上图中,只考虑第一卦限,蓝色线就是我们希望求参数方程的曲线,其中空间曲线与原来椭圆平面曲线的对应关系,已经用了跟前一图(右图)同一种颜色来表示。$P=(x(s),y(s),z(s))$为曲线上的一点,假设弧长$OP$为$s$,这段弧长实际上也就是椭圆上原来的一段弧长。而前面我们已经讨论过了,根据弧长可以求出$y$坐标
\begin{equation}y=h(s)\end{equation}
然后由于对称(互补)关系,$z$坐标事实上是弧长为$l-s$时对应的$y$坐标,即
\begin{equation}z=h(l-s)\end{equation}
其中$l$是整个椭圆周长的1/4,也就是蓝色曲线的总长度。
现在$y,z$都有了,还缺$x$,因为$s$是弧长参数,所以必然满足$dx^2+dy^2+dz^2=ds^2$,即
\begin{equation}\left(\frac{dx(s)}{ds}\right)^2+\left(\frac{dh(s)}{ds}\right)^2+\left(\frac{dh(l-s)}{ds}\right)^2=1\end{equation}
这样就可以把$x(s)$积分出来:
\begin{equation}x(s) = \int_0^s \sqrt{1-\left(\frac{dh(s)}{ds}\right)^2-\left(\frac{dh(l-s)}{ds}\right)^2}ds\label{eq:f}\end{equation}
而
\begin{equation}1-\left(\frac{dh(s)}{ds}\right)^2-\left(\frac{dh(l-s)}{ds}\right)^2 \geq 0\label{eq:c}\end{equation}
是否恒成立,就是是否能粘合的判别条件。
证明可以粘合 #
下面我们从理论上证明:两个椭圆片是可以完美贴合的,也就是验证$\eqref{eq:c}$恒成立。
首先,我们需要证明一个引理:
在椭圆$x^2/a^2 + y^2/b^2=1(a \geq b > 0)$中有
\begin{equation}\left.\left(\frac{dy}{dx}\right)^2\right|_{s=s_0}\leq \left.\left(\frac{dx}{dy}\right)^2\right|_{s=l-s_0}\end{equation}
也就是弧长为$s_0$处的斜率绝对值,不超过弧长为$l-s_0$处的斜率倒数的绝对值。
证明过程如下(并不难,手动画图,比划一下就明白了):
利用椭圆常见的参数方程
\begin{equation}x=a\cos t,\quad y=b\sin t\end{equation}
注意这个参数$t$增加的方向,跟$s$增加的方向是是一样的,都是从$(a,0)$出发逆时针方向,因此$s$与$t$也是单调递增关系。可以得到
\begin{equation}\left.\left(\frac{dy}{dx}\right)^2\right|_{t=t_0} = \frac{b^2\cos^2 t_0}{a^2\sin^2 t_0}\leq \frac{a^2\sin^2 (\pi/2-t_0)}{b^2\cos^2 (\pi/2-t_0)}=\left.\left(\frac{dx}{dy}\right)^2\right|_{t=\pi/2-t_0}\end{equation}
也就是说,$t=t_0$处的斜率,肯定不超过$t=\pi/2-t_0$处的斜率倒数的绝对值,并且很显然,$t_0$越小时$\left.\left(\frac{dx}{dy}\right)^2\right|_{t=\pi/2-t_0}$越大。然后我们考虑
\begin{equation}\frac{ds}{dt}=\sqrt{\left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2}=dt\sqrt{a^2 - (a^2-b^2)\cos^2 t}\end{equation}
注意,这是一个$t$的单调递增函数(都只考虑$t\in[0,\pi/2]$区间),这意味着$s$随着$t$的增加而增加,并且增加的速度也越来越快,这意味着
\begin{equation}s_{[0,t_0]}\leq s_{[\pi/2-t_0,\pi/2]}\end{equation}
也就是说$[0,t_0]$区间的弧长,不超过$[\pi/2-t_0,\pi/2]$区间的弧长,这又意味着如果$s_{[0,t_0]} = s_{[\pi/2-t_1,\pi/2]}$,则$t_1 \leq t_0$。换言之,如果$s=s_0$对应$t=t_0$、$s=l-s_0$对应$t=t_1$的话,则$t_1 \leq t_0$,根据刚才的讨论:
\begin{equation}\left.\left(\frac{dy}{dx}\right)^2\right|_{s=s_0}=\left.\left(\frac{dy}{dx}\right)^2\right|_{t=t_0} \leq \left.\left(\frac{dx}{dy}\right)^2\right|_{t=\pi/2-t_1}=\left.\left(\frac{dx}{dy}\right)^2\right|_{s=l-s_0}\end{equation}
有了这个引理之后,我们去证明:
\begin{equation}\left(\frac{dh(s)}{ds}\right)^2+\left(\frac{dh(l-s)}{ds}\right)^2\leq 1\end{equation}
这里$h(s)$本来就是原来椭圆的$y$坐标,那么假设对应的$x$坐标为$g(s)$,则:
\begin{equation}\left(\frac{dh(s)}{ds}\right)^2+\left(\frac{dg(s)}{ds}\right)^2 = 1\end{equation}
那么等价于证明
\begin{equation}\left(\frac{dh(l-s)}{ds}\right)^2\leq\left(\frac{dg(s)}{ds}\right)^2\end{equation}
两边取倒数,然后利用$ds^2=dg^2+dh^2$,事实上就是要证明
\begin{equation}\left.\left(\frac{dg}{dh}\right)^2\right|_{s=l-s_0}\geq\left.\left(\frac{dh}{dg}\right)^2\right|_{s=s_0}\end{equation}
这就是刚才的引理。
数值求解结果 #
从数学上来看,$\eqref{eq:f}$已经称得上所求曲线的显式解(解析解)了,当然显式解不意味着是初等函数解,也不意味着这样的解答容易算。虽然微分几何中大多数结果都是以弧长为参数的方程,但那都是理论上的,事实上除了圆和直线,就没几条曲线的弧长参数方程是简单的。
为了数值验证条件$\eqref{eq:c}$以及求解$\eqref{eq:f}$,我发现现成的数学软件(如Mathematica)并不好用,(对于笔者来说)最简单的方法还是直接通过Python写一段代码来模拟求解。以$a=2,b=1$为例,它通过了$\eqref{eq:c}$的数值验证,并且得到如下的解:
代码如下:
#! -*- coding: utf-8 -*-
import numpy as np
a = 2
b = 1
h = 1e-6 # 步长
def Int(y, x):
"""自定义数值积分函数
"""
assert len(y) == len(x)
dx = np.diff(x)
_x = (x[:-1] + x[1:]) / 2
_y = (y[:-1] + y[1:]) / 2
s = np.cumsum(_y * dx)
s = np.interp(x, _x, s, left=0)
return s, x
def D(y, x):
"""自定义数值微分函数
"""
assert len(y) == len(x)
dy = np.diff(y) + 1e-10
dx = np.diff(x) + 1e-10
dy_dx = dy / dx
_x = (x[:-1] + x[1:]) / 2
dy_dx = np.interp(x, _x, dy_dx)
return dy_dx, x
def s2y():
"""弧长到椭圆y轴的映射(即h(s)。)
"""
y = np.arange(0, b, h)
ds = np.sqrt(1 + a**2 * y**2 / b**2 / (b**2 - y**2))
return Int(ds, y)[::-1]
y, s = s2y()
l = s[-1]
z = np.interp(l - s, s, y)
dx2 = 1 - D(y, s)[0]**2 - D(z, s)[0]**2
assert (dx2 > - 10 * h).all() # 判断是否能粘合,即在10倍步长的误差内是否恒大于等于0
dx = np.sqrt(np.clip(dx2, 0, a))
x = Int(dx, s)[0]
# 采样1/10的点去绘图
x = x[::10]
y = y[::10]
z = z[::10]
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
plt.clf()
fig = plt.figure(figsize=(8, 8))
ax = fig.gca(projection='3d')
ax.plot(x, y, z, color='blue')
ax.plot(x, - y, z, color='blue')
ax.plot(x, y, - z, color='green')
ax.plot(x, - y, - z, color='green')
ax.plot(x, np.zeros_like(x), z, color='red')
ax.plot(x, np.zeros_like(x), - z, color='red')
for i in range(0, len(x), 2000):
_ = ax.plot([x[i], x[i]], [y[i], - y[i]], [z[i], z[i]], color='blue')
_ = ax.plot([x[i], x[i]], [y[i], - y[i]], [- z[i], - z[i]], color='green')
for i in [-1, -1000]:
_ = ax.plot([x[i], x[i]], [y[i], - y[i]], [z[i], z[i]], color='blue')
if i != -1:
_ = ax.plot([x[i], x[i]], [y[i], - y[i]], [- z[i], - z[i]], color='green')
ax.set_xlim(0, a)
ax.set_ylim(-b * 1.2, b * 1.2)
ax.set_zlim(-b * 1.2, b * 1.2)
plt.show()
小结 #
本文分析了两个椭圆片的粘合问题,其核心思想是用弧长参数化曲线方程,虽然文章是以椭圆为例,但显然可以推广到一般的两个平面曲线面的弯曲粘合问题,只需要分别求出各自的弧长参数方程(即$h(l-s)$换成另一条曲线的弧长参数化方程),然后沿用$\eqref{eq:f}$式积分。
总算基本解决了一个比较有意思的问题,来证明自己的数学兴趣和能力还在~
转载到请包括本文地址:https://kexue.fm/archives/6818
更详细的转载事宜请参考:《科学空间FAQ》
如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。
如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!
如果您需要引用本文,请参考:
苏剑林. (Jul. 21, 2019). 《思考:两个椭圆片能粘合成一个立体吗? 》[Blog post]. Retrieved from https://kexue.fm/archives/6818
@online{kexuefm-6818,
title={思考:两个椭圆片能粘合成一个立体吗?},
author={苏剑林},
year={2019},
month={Jul},
url={\url{https://kexue.fm/archives/6818}},
}
July 22nd, 2019
小苏同学干的不错,搞出一个判别式来了!
记得有一个曲率是弧长的线性函数的曲线,其弧长参数方程是一个简明的积分式。
谢谢大佬。大佬到此捧场,受宠若惊呀。
除此之外好像悬链线的弧长参数方程也不算复杂,反正这类例子也不多。
那么多测地线方程,没几个真正有解析解出来的。
July 24th, 2019
你是我见过的理论、编程与讲解结合的最好的学者了
谢谢~^_^
July 31st, 2019
问题来了,围成的空间体积是多少
目测只能算个数值结果~
基本计算公式是
$$4\int_0^l h(s)h(l-s)\sqrt{1-h'(s)^2-h'(l-s)^2}ds$$
August 19th, 2019
剑林,我想纠正个小问题: 在弧长参数方程里面, 应该把逆函数写成反函数, 函数是可逆的(reversible), 而没有逆函数.
OK,接受并感谢你这个建议,已修正~
November 26th, 2019
请问楼主“粘贴边缘方程”一节中两幅图是用什么制图软件画的?
这个就是用powerpoint做的而已~
November 27th, 2019
最后一段说得很对,不但椭圆可以,一类曲线都可以,比如:$(\frac{x}{a})^p+(\frac{y}{b})^p=1$,其中$1\leq p\leq 2$.
错啦哈,请忽略!$p=1$或$p=2$,前者就是一个菱形而已。