拉普拉斯贝尔特拉米算子

拉普拉斯贝尔特拉米算子（Laplace–Beltrami operator）是在具有黎曼度量的空间上，对标量函数先取梯度、再取散度得到的二阶微分算子。它把平面上的 $f_{x x} + f_{y y}$ 推广到曲面及更一般的黎曼流形；比较相邻位置时，长度、面积和方向都由度量决定。

本文采用 $Δ_{g} = {div}_{g} {grad}_{g}$ 的符号。按这一约定，圆周上的余弦函数被送到自身的负倍数，热方程写成 $u_{t} = κ Δ_{g} u$ ，其中扩散系数 $κ > 0$ 。有些文献把负号并入算子，使用 $- Δ_{g}$ ；比较特征值时须同时核对这个约定。

从圆周上的温度变化开始

半径为 $R > 0$ 的圆周可写成 $X (θ) = (R \cos θ, R \sin θ)$ 。角度 $θ$ 用弧度表示，沿圆周走过的弧长为 $s = R θ$ 。设温度只随圆周位置变化，记为 $f (θ)$ 。相同的角度变化在大圆上对应更长的路程，所以沿弧长求导时 $\frac{d}{d s} = \frac{1}{R} \frac{d}{d θ}, Δ_{g} f = \frac{d^{2} f}{d s^{2}} = \frac{1}{R^{2}} \frac{d^{2} f}{d θ^{2}} .$ 这里讨论的是圆周本身这个一维空间，不是圆盘内的二维温度场。例如 $f (θ) = \cos (m θ), Δ_{g} f = - \frac{m^{2}}{R^{2}} f (m = 0, 1, 2, \dots) .$ $m$ 必须是整数，才能使函数在绕行一周后与原值及各阶导数相接。峰谷越密， $m$ 越大，二阶变化越强；保持相同的角度分布而将半径加倍，算子的结果缩小为原来的四分之一。

梯度、体积元与坐标公式

设 $M$ 是无边界的光滑 $n$ 维流形， $g$ 是其光滑正定黎曼度量。以下经典公式用于 $f \in C^{2} (M)$ 。在局部坐标 $x^{1}, \dots, x^{n}$ 中，用 $\partial_{i} = \partial / \partial x^{i}$ 表示坐标方向，记 $g_{i j} = g (\partial_{i}, \partial_{j}), G = det (g_{i j}) > 0, (g^{i j}) = (g_{i j})^{- 1} .$ 上标 $i j$ 在这里表示逆矩阵的元素，并非对 $g_{i j}$ 逐项取倒数。度量给出的体积元为 $d V_{g} = \sqrt{G} d x^{1} \dots d x^{n} .$ 二维情形就是面积元。它能在局部写出并拼合为全局体积密度，不要求流形可定向。

梯度 ${grad}_{g} f$ 是满足 $g ({grad}_{g} f, V) = V (f)$ 的唯一切向量场； $V (f)$ 表示沿向量 $V$ 的方向导数。若将梯度写成 $\sum_{i} a^{i} \partial_{i}$ ，分别取 $V = \partial_{j}$ ，便有 $\sum_{i} g_{i j} a^{i} = \partial_{j} f$ 。解这个线性方程组得到 $a^{i} = \sum_{j} g^{i j} \partial_{j} f .$ 这一步说明，函数的偏导数与梯度的坐标分量通常不同。

对一个 $C^{1}$ 切向量场 $V = \sum_{i} V^{i} \partial_{i}$ ，散度表示单位几何体积内的净流出率。它的坐标表达式为 ${div}_{g} V = \frac{1}{\sqrt{G}} \sum_{i} \partial_{i} (\sqrt{G} V^{i}) .$ 把梯度的分量代进去，即得 $Δ_{g} f = \frac{1}{\sqrt{G}} \sum_{i, j = 1}^{n} \partial_{i} (\sqrt{G} g^{i j} \partial_{j} f) .$ 度量同时出现在逆矩阵与体积权重中：前者决定方向导数如何变成梯度，后者决定流出量按什么体积平均。只把普通二阶偏导数相加，一般不能得到曲面上的这个算子。^[1]

公式中的体积权重

可以用分部积分法检验散度公式，而不靠画面猜测它。取支撑紧含在一个坐标邻域内的光滑函数 $φ$ ，故积分时没有边界项。逐个坐标积分得到 $- \int_{M} V (φ) d V_{g} = - \int \sum_{i} V^{i} \partial_{i} φ \sqrt{G} d x = \int φ \sum_{i} \partial_{i} (\sqrt{G} V^{i}) d x .$ 将最后一个积分改回 $d V_{g}$ ，就得到 $\int_{M} φ {div}_{g} V d V_{g} = - \int_{M} V (φ) d V_{g} .$ 任意这样的 $φ$ 都满足此式，因此右侧确定散度的值。左侧若少掉 $1 / \sqrt{G}$ ，就不能成立。这种表达还说明散度与坐标的选择无关：方向导数、几何体积和积分本身都不依赖坐标。

平面极坐标是一个可直接复核的例子。在 $r > 0$ 的坐标邻域内， $g = (\begin{matrix} 1 & 0 \\ 0 & r^{2} \end{matrix}), \sqrt{G} = r, Δ f = \frac{1}{r} \partial_{r} (r \partial_{r} f) + \frac{1}{r^{2}} \partial_{θ}^{2} f = f_{r r} + \frac{1}{r} f_{r} + \frac{1}{r^{2}} f_{θ θ} .$ 其中 $f_{r} / r$ 来自半径增大时扇形面积与圆弧长度的变化。比如同一个函数 $f = x^{2} + y^{2} = r^{2}$ ，在直角坐标中得到 $2 + 2 = 4$ ；极坐标中得到 $2 + 2 r / r = 4$ 。若删去中间的一阶导数项，就只剩下 $2$ 。^[2]

半径一与二处同样角宽和径向厚度的两块环扇形，外圈圆弧更长 — 固定角宽和径向厚度，外圈的环扇形面积更大。极坐标中的体积权重为 r；图中的有限小块面积还包含径向厚度的二阶项。

图中的环扇形取 $Δ θ = π / 6$ 、 $Δ r = 1 / 4$ 。内半径为 $r$ 时，其精确面积是 $\frac{1}{2} ((r + Δ r)^{2} - r^{2}) Δ θ = (r Δ r + \frac{(Δ r)^{2}}{2}) Δ θ .$ 在 $r = 1$ 与 $r = 2$ 处分别为 $3 π / 64$ 、 $17 π / 192$ 。微分面积元 $r d r d θ$ 描述一阶极限，不应把它当成任意有限小块的精确面积。

圆柱上的完整计算

考虑半径 $R$ 的圆柱侧面 $X (θ, z) = (R \cos θ, R \sin θ, z) .$ 沿两个坐标方向求导得 $X_{θ} = (- R \sin θ, R \cos θ, 0)$ 、 $X_{z} = (0, 0, 1)$ 。它们垂直，长度分别为 $R$ 与 $1$ ，所以 $(g_{i j}) = (\begin{matrix} R^{2} & 0 \\ 0 & 1 \end{matrix}), \sqrt{G} = R, Δ_{g} f = \frac{1}{R^{2}} f_{θ θ} + f_{z z} .$ 取 $R = 2$ ，考察圆柱上的标量函数 $f (θ, z) = \cos (2 θ) + z^{2} .$ 这里使用无量纲教学坐标；若把函数当成物理温度模型，各项还应带有相应单位的系数。两个二阶偏导分别为 $- 4 \cos (2 θ)$ 与 $2$ ，故 $Δ_{g} f = - \cos (2 θ) + 2.$ 在 $θ = π / 3, z = 1$ 处， $f = 1 / 2$ 、 $Δ_{g} f = 5 / 2$ 。函数值与算子值是不同的量：前者是该点的读数，后者描述周围变化的综合。

把圆柱展开为平面，令 $s = 2 θ$ ，同一个函数变成 $f (s, z) = \cos s + z^{2}$ ，直接计算 $f_{s s} + f_{z z}$ 仍得 $- \cos s + 2$ 。圆柱在三维中弯曲，但这种局部展开保留长度，算子也随之保留。全局上则须让圆周方向周期相接，不能将展开图左右边当成两个自由边界。

变式。 保持同一个角度函数 $\cos (2 θ) + z^{2}$ ，改用 $R = 4$ 。此时 $Δ_{g} f = - \frac{1}{4} \cos (2 θ) + 2, Δ_{g} f (π / 3, 1) = \frac{17}{8} .$ 圆周方向的变化被更长的实际距离摊开，轴向的 $2$ 则不变。

球面上的高度函数

半径 $R$ 的球面用余纬 $θ \in (0, π)$ 、经度 $ϕ$ 参数化： $X (θ, ϕ) = (R \sin θ \cos ϕ, R \sin θ \sin ϕ, R \cos θ) .$ 这里 $θ = 0$ 对应北极，不能与从赤道量起的纬度混用。两个坐标方向的长度为 $R$ 和 $R \sin θ$ ，因而 $Δ_{g} f = \frac{1}{R^{2} \sin θ} \partial_{θ} (\sin θ \partial_{θ} f) + \frac{1}{R^{2} \sin^{2} θ} \partial_{ϕ}^{2} f .$ 对无量纲高度 $f = \cos θ = z / R$ ，经度方向的导数为零；另一方面 $\partial_{θ} f = - \sin θ, \partial_{θ} (\sin θ \partial_{θ} f) = \partial_{θ} (- \sin^{2} θ) = - 2 \sin θ \cos θ .$ 故 $Δ_{g} f = - \frac{2}{R^{2}} f .$ 例如 $R = 2$ ，在 $θ = π / 3$ 处， $f = 1 / 2$ 、 $Δ_{g} f = - 1 / 4$ 。球面公式的分母在极点为零，源于经纬坐标在那里失效；高度函数和算子值本身光滑，以上等式由连续性延伸到两极。

余纬从零到π时高度cosθ与半径二球面上的算子值负二分之一cosθ的曲线 — 半径为二时，球面归一化高度 f=cosθ 经算子作用后变为 −f/2。赤道处两者为零，北半球与南半球的符号相反。

分部积分、符号与边界条件

在紧致、无边界的黎曼流形上，对光滑实函数 $u, v$ ，散度的乘积法则给出 ${div}_{g} (v {grad}_{g} u) = g ({grad}_{g} v, {grad}_{g} u) + v Δ_{g} u .$ 对全空间积分，散度定理使左侧为零，因此 $\int_{M} v Δ_{g} u d V_{g} = - \int_{M} g ({grad}_{g} v, {grad}_{g} u) d V_{g} .$ 取 $v = u$ ，右侧非正。这说明 $Δ_{g}$ 的二次型非正；若 $Δ_{g} u = λ u$ 且 $u$ 非零，则 $λ \leq 0$ 。常函数对应零特征值。若空间还连通，光滑函数满足 $Δ_{g} u = 0$ ，上述等式迫使梯度处处为零，沿任意路径函数值不变，所以函数只能是常数。

带有光滑边界时会多出 $\int_{\partial M} v \partial_{ν} u d S_{g}, \partial_{ν} u = g ({grad}_{g} u, ν),$ 其中 $ν$ 是位于流形切空间内、垂直其边界并指向外侧的单位向量。它不是嵌入曲面在三维空间中的法向。指定函数值的 Dirichlet 条件、指定外向导数的 Neumann 条件，会产生不同问题。局部微分表达式相同，不意味着其边界问题或谱相同。^[1]在非紧致空间上，积分也须具备紧支撑或足够衰减等条件，不能直接删去无穷远处的边界贡献。

扩散中的尺度与频率

若理想介质只沿空间内部扩散，且扩散系数恒为 $κ > 0$ ，可用 $u_{t} = κ Δ_{g} u$ 描述相对于恒定基准的温差。对半径 $R$ 的圆周，初始温差为 $\cos (m θ)$ 时，试取 $u = A (t) \cos (m θ)$ ，代入得到 $A^{'} (t) = - \frac{κ m^{2}}{R^{2}} A (t), A (0) = 1, u (θ, t) = e^{- κ m^{2} t / R^{2}} \cos (m θ) .$ 这是直接代回方程可验证的解，不需要先掌握完整的谱理论。 $m = 0$ 的常量部分保持不变；高频分量更快衰减。 $κ$ 的单位为长度平方除以时间，故指数 $κ m^{2} t / R^{2}$ 无量纲。

同一时间下圆周一阶与三阶余弦温差的衰减振幅，以及半径加倍的一阶温差 — 取 κ=1、t=1/4：R=1 的一阶与三阶分量振幅分别为 exp(−1/4) 和 exp(−9/4)，R=2 的一阶分量振幅为 exp(−1/16)。图中横轴是角度，各曲线对应各自的初始余弦分量。

在紧致无边界空间上，若 $u$ 是足够光滑的热方程解，积分公式还给出 $\frac{d}{d t} \int_{M} u d V_{g} = 0, \frac{d}{d t} \frac{1}{2} \int_{M} u^{2} d V_{g} = - κ \int_{M} | {grad}_{g} u |_{g}^{2} d V_{g} \leq 0.$ 前式表示总量守恒，后式说明平方幅度不会因扩散而增长。它们也为数值结果提供检查依据。

圆周上的离散计算

用 $N \geq 3$ 个等角节点 $θ_{j} = j h$ ， $h = 2 π / N$ ，下标按模 $N$ 周期循环。将圆周二阶导数换成中心差分，再将时间导数换成向前差分，得到 $u_{j}^{n + 1} = u_{j}^{n} + α (u_{j - 1}^{n} - 2 u_{j}^{n} + u_{j + 1}^{n}), α = \frac{κ τ}{R^{2} h^{2}} .$ 其中 $τ$ 为时间步长。空间差分对足够光滑函数的截断误差为 $O (h^{2})$ ，时间差分为 $O (τ)$ ：将两侧节点作 Taylor 展开，奇数阶项抵消，留下 $u_{θ θ} + h^{2} u_{θ θ θ θ} / 12 + \dots$ ；时间展开则留下 $u_{t} + τ u_{t t} / 2 + \dots$ 。这些局部误差阶不等于任意网格上已给出了整体误差数值。

将更新式写成 $u_{j}^{n + 1} = α u_{j - 1}^{n} + (1 - 2 α) u_{j}^{n} + α u_{j + 1}^{n},$ 就能检查一个明确的稳定条件： $0 \leq α \leq 1 / 2$ 时，三个权重非负且和为一，新值落在旧值的最小值和最大值之间。按全部节点求和，周期循环使左右邻项的总和都等于原总和，所以离散平均值不变。

四节点圆周的周期相邻关系与节点零在四分之一权重下的更新 — N=4、α=1/4 时，节点零接收自身的二分之一和两个相邻节点各四分之一；节点三与节点零相邻。每步要先保存全部旧值，不能用刚更新的值接着计算。

取 $N = 4$ 、 $α = 1 / 4$ 和初始向量 $(1, 0, - 1, 0)$ 。节点零第一步为 $(0 + 2 \times 1 + 0) / 4 = 1 / 2$ ，节点一为 $(1 + 2 \times 0 - 1) / 4 = 0$ ，其他两点同理。连续三步得到 $(1, 0, - 1, 0) ⟼ (1 / 2, 0, - 1 / 2, 0) ⟼ (1 / 4, 0, - 1 / 4, 0) ⟼ (1 / 8, 0, - 1 / 8, 0) .$ 这是离散格式的精确结果，不是连续方程解的精确采样。例如 $R = κ = 1$ 时，一步对应 $τ = π^{2} / 16$ ，连续一阶模态的振幅应为 $e^{- π^{2} / 16}$ ，约 $0.540$ ，并非 $1 / 2$ 。细化空间时须同时缩小时间步长，使 $α$ 保持在稳定范围内。

以下独立实现采用整数或 Fraction 输入，精确演示有限网格更新；不接收浮点数，也不直接求任意曲面上的算子。函数返回最终节点值，保留输入序列不变。

Python 3

from fractions import Fraction


def circle_heat_steps(values, alpha, steps):
    def rational(value):
        if isinstance(value, bool) or not isinstance(value, (int, Fraction)):
            raise TypeError("Use integers or Fraction values")
        return Fraction(value)

    if isinstance(steps, bool) or not isinstance(steps, int):
        raise TypeError("steps must be an integer")
    if steps < 0:
        raise ValueError("steps must be nonnegative")
    a = rational(alpha)
    if not 0 <= a <= Fraction(1, 2):
        raise ValueError("alpha must lie in [0, 1/2]")
    u = [rational(value) for value in values]
    if len(u) < 3:
        raise ValueError("At least three periodic nodes are required")
    size = len(u)
    for _ in range(steps):
        old = u
        u = [(1 - 2 * a) * old[j]
             + a * (old[(j - 1) % size] + old[(j + 1) % size])
             for j in range(size)]
    return u


result = circle_heat_steps([1, 0, -1, 0], Fraction(1, 4), 3)
print([str(value) for value in result])
# ['1/8', '0', '-1/8', '0']

每一轮都只读取 old 的值，所以实现的正是同步更新式。由凸组合性质和周期求和，上述范围约束及平均值守恒在每一轮继续成立；归纳即得全部步数的结论。 $α = 0$ 或零步时原值保持不变，常数输入也保持不变。执行 $K$ 步需要 $O (N (K + 1))$ 次有理数基本运算和 $O (N)$ 个额外有理数存储位置；分子分母的位数可随步数增加，故这里不是忽略位长后仍成立的固定比特成本。

越界变式。 若不加限制而取 $α = 3 / 4$ ，初值 $(1, - 1, 1, - 1)$ 在一步后变成 $(- 2, 2, - 2, 2)$ 。因为两侧邻值都与本点相反，每个分量都乘以 $1 - 4 α = - 2$ ；继续迭代会放大振荡。代码因此拒绝该步长。

在一般曲面上，节点分布、面积权重和几何近似都要参与离散化，不能直接把圆周三点格式搬到任意邻接图上。有限元方法会从分部积分式建立质量矩阵和刚度矩阵，并同时处理函数近似与曲面几何近似造成的误差。^[1]本文的圆周代码只用于看清权重、尺度和时间步长之间的关系。

参考文献

↑ ^1.0 ^1.1 ^1.2 Andrea Bonito、Alan Demlow、Ricardo H. Nochetto：Finite Element Methods for the Laplace-Beltrami Operator，arXiv:1906.02786v1，2019，§2.1–2.2，尤其式 (8)、(10)、(12)、(13) 与引理 1。该文从参数曲面的度量、面积元构造切向微分算子。本文对一般黎曼度量的推导及教学例子独立展开。
↑ Victor Ivrii：Partial Differential Equations, §6.3: Laplace operator in different coordinates，University of Toronto，式 (6.3.3)、(6.3.8)、(6.3.9)，极坐标、球面角向算子及度量坐标公式。访问日期：2026年10月25日。

[bdn-1] 1.0 ^1.1 ^1.2 Andrea Bonito、Alan Demlow、Ricardo H. Nochetto：Finite Element Methods for the Laplace-Beltrami Operator，arXiv:1906.02786v1，2019，§2.1–2.2，尤其式 (8)、(10)、(12)、(13) 与引理 1。该文从参数曲面的度量、面积元构造切向微分算子。本文对一般黎曼度量的推导及教学例子独立展开。

[ivrii-2] Victor Ivrii：Partial Differential Equations, §6.3: Laplace operator in different coordinates，University of Toronto，式 (6.3.3)、(6.3.8)、(6.3.9)，极坐标、球面角向算子及度量坐标公式。访问日期：2026年10月25日。

[1]

[2]