🗒️微分算子及其离散化4：拉普拉斯
1 定义
前面介绍了梯度和散度。梯度作用于标量场，得到向量场；散度作用于向量场，得到标量场。如果先对一个标量函数求梯度，再对这个梯度场求散度，就得到了拉普拉斯算子（Laplacian）：
拉普拉斯算子常记作 或 。这里的 表示梯度的散度，不能理解为梯度向量的长度平方。
以二维标量函数 为例：
因此，拉普拉斯算子作用于标量场，输出的仍然是标量场。三维情况只需再加上 。下文也会把这几个二阶偏导数简写为 、、。
这些二阶导数相加，在几何上表示什么？首先从一维函数开始看。
2 一维：比较中心与两侧
一维函数只有一个方向，此时 。一阶导数描述斜率，二阶导数描述斜率如何变化。
首先考虑 ，在 附近取间距 的三个点：
如下图，把左右两个点连起来。由于两点到中心的水平距离相等，连线在 处的高度，就是左右函数值的平均数：
它比中心的函数值高 。对于 ，左右两点都变为 ，中心仍为 ，平均值与中心的差就变为 。
图1：x²与2x²，邻点连线中点的高度分别比中心高1和2
图中黄色点是三个采样点，空心点是邻点连线的中点，红色线段标出它与中心的高度差。两个函数在原点的斜率都是 ，但向左右走时， 的斜率变化得更快。
如果把函数换成 ，左右平均值变为 ，比中心低 。再看斜直线 ：左、中、右的值分别是 ，左右平均值恰好等于中心值。
图2：−x²的邻点平均低于中心，斜直线的邻点平均等于中心
所以，函数沿着一个方向增加，并不代表二阶导数为正。直线可以一直增加，但它的斜率不变，二阶导数仍为 。
从两段斜率推导三点差分
上面只比较了高度。要得到二阶导数，还需要考虑点之间的距离。
在均匀网格上，记 ，相邻点间距为 。首先连接中心与右邻点，这段割线的斜率为：
把它作为右侧半格位置 的导数近似。中心与左邻点同理：
这里两段斜率都按 增大的方向计算：左段也是右端的值减左端的值。
图3：先计算左右两段斜率，再在两个半格位置之间做差分
这两个斜率所代表的位置相距 。用右侧斜率减左侧斜率，再除以这个距离，就得到了斜率的变化率：
这就是一维拉普拉斯的三点差分格式。分母里的两个 分别来自求斜率和求斜率的变化率。
记这个离散结果为 。如果左右邻点的平均值为 ，上式也可以写成：
回到 的四个例子，得到的结果依次为 ，恰好等于各自的二阶导数。对一般函数，有限间距上的差分只是近似；当函数足够光滑且 时，它趋向 。
3 二维：把两个方向的结果相加
对于 ，可以把函数值画成曲面的高度。取定一点 ，首先固定 ，沿 方向切出一条曲线。这条曲线在该点的二阶导数，就是 。
再固定 ，沿 方向切出另一条曲线，得到 。两个结果相加：
以原点为例，先看 。沿 方向的截线是 ，沿 方向的截线是 。两条截线都向上弯曲，各自贡献 ，因此 。
如果换成 ，两个方向都向下弯曲，各自贡献 ，因此 。
可以在下面切换函数、旋转观察，调整中心位置与步长，并查看连续值和离散值。数值参数以演示面板为准。
静态 SVG 参考｜图4：向上的碗面和向下的碗面，两个方向的二阶导数同号
图中蓝色是沿 方向的截线，红色是沿 方向的截线，黄色点是原点。灰色网线表示曲面的其余部分。
接下来考虑一个斜平面 。虽然沿着两个坐标方向走，函数值都在增加，但两条截线都是直线，所以 ，拉普拉斯也为 。
再看 。沿 方向的截线向上弯曲，贡献 ；沿 方向的截线向下弯曲，贡献 。相加后仍然是 ：
静态 SVG 参考｜图5：斜平面与鞍面都满足拉普拉斯为零，但截线的形状不同
因此，拉普拉斯为零并不表示曲面是平的。它只要求这两个方向的二阶导数之和为零，允许一正一负相互抵消。
同样，不能只看拉普拉斯的正负来判断局部极值。例如 ，拉普拉斯为 ；但从原点沿 方向走，函数值上升，沿 方向走，函数值下降。原点仍然是鞍点，并不是局部最小值。
4 从梯度场的流出量理解
前面从二阶导数观察了曲面的形状。也可以回到定义 ：先画出梯度场，再考察这个向量场在一个小区域内的净流出量。
注意， 的定义域是二维的 平面。上一节的曲面只是用高度显示函数值。梯度 是定义在这个二维平面上的向量，下面的方块也取在 平面上。
以 为中心，取一个边长为 的方块。它的右、左、上、下边中点，分别记为 ：
图6：二维定义域中的控制方块及四条边的外法向
记梯度场为 。首先看右边界：只有垂直于边界的 分量穿过这条边，平行于边界的 分量不贡献流出量。若方块很小，用边中点的值近似整条边上的值，右边的流出量就是 。
左边的外侧是 减小的方向。因此左边流出量为 ：当 时，向量朝右，实际是流入方块，要记负号。
同理，上边流出量为 ，下边为 。四边相加得到整个方块的净流出量：
这个量会随着方块的大小而变。要得到单位面积的净流出量，需要除以方块面积 ：
右边的两个分式，就是 沿 方向、 沿 方向的变化率。当 时，得到中心处的散度，也就是拉普拉斯：
这里 表示四条边上梯度场的实际净通量；边中点取值给出的是它的近似。由此也能看出，拉普拉斯描述的是净流出量的面积密度，不能直接把一个有限方块的总流出量当成某点的拉普拉斯值。
再看原点附近的三个例子：
图7：同一个控制方块上，谷底梯度向外、峰顶向内、鞍面两个方向抵消
对于 ，梯度为 。在右边中点 处，，所以右边流出量为 ；左边中点 ，流出量同样为 。上下两边也是如此，总量为 ，除以面积后得到 。
对于 ，梯度指向中心，四边都贡献负的流出量，得到 。对于 ，左右两边向外，上下两边向内，正负贡献相互抵消，得到 。
图中画的是梯度场。非零梯度指向函数值增长最快的方向；它还不是热流，热流的方向会在后面单独讨论。
5 离散化：从交错网格得到五点格式
现在把上面的过程放到网格上。考虑二维均匀方格，横纵方向的格子中心间距都是 。把标量 存在格子中心，当前格子及其四个邻居的值记为：
这里 分别表示右、左、上、下的邻居中心。首先要做的，是由这些中心值算出格子各条边上的梯度分量。
先把梯度放在面中心
沿用散度一篇中的交错网格（MAC grid）： 放在左右两边的中点， 放在上下两边的中点。二维里的“面中心”，就是一条边的中点。图中将右、左、上、下四个面上的分量简写为 。
图8：标量存于格子中心，梯度的法向分量存于对应面中心
当前中心的编号为 。右面的编号是 ，因为它位于当前中心与右邻居的正中间。用这两个中心值做差：
这是右面中心处 的近似。注意，虽然右面到当前中心的距离是 ，但用来求差的两个采样点是 和 ，它们相距 ，所以分母为 。
左面同理。仍然按从左到右的方向计算斜率：
图9：中心之间相距h，左右两个面中心之间也相距h
对于上下两面，按从下到上的方向计算：
这里用等号定义离散梯度 ；它对连续梯度的近似含义并没有消失。
再把散度算回格子中心
已经有了四个面上的分量，就可以使用上一节的通量公式。左右两个面的中点相距 ，上下也一样。因此中心处的离散散度为：
首先把左右两项代入：
上下两项同理，得到 。再把两个方向相加：
这就是二维拉普拉斯的五点差分格式：用到当前中心和上下左右四个邻居，四个邻居的系数都是 ，中心的系数是 。
式子里的 也可以从推导中直接数出来。每个方向求一次二阶差分，都产生 ；两个方向相加，就是 。
对于足够光滑的函数，这个离散结果近似中心处的连续拉普拉斯。在内部点使用上述均匀网格格式时，空间误差为 ：网格间距减半，主导误差项约变为原来的四分之一。
为什么要说明存储位置
连续的 可以直接写在一起，但离散时，两次求差分别用哪些点，会影响最后的格式。
例如，在一维中，如果先用同位中心差分 ，再对 使用同样的中心差分，得到的是：
它用到了相距两个格子的点，是另一套更宽的差分格式。上面的三点、五点格式来自中心与面中心的交错：中心值求面上梯度，面上梯度求中心散度。
单独打开曲面交互实验，调整中心位置与步长，比较连续值与离散值。
6 邻域平均：周围平均值减去中心
五点公式还可以换一种写法。记四个邻居的平均值为：
将邻居之和替换为 ，得到：
因此，在这个离散格式下，拉普拉斯的正负取决于四邻平均值与中心值的差。系数 来自四个邻居， 来自前面的两次差分。
先取 ，直接代入几组数：
- 中心为 ，四邻都为 。平均值 ，所以 。
- 中心为 ，四邻都为 。平均值仍为 ，所以 。
- 中心为 ，左、右、上、下分别为 。平均值为 ，所以 。
图10：四组五点数据，对比中心值与四邻平均值
第三组里，中心有一个更低的邻居，也有一个更高的邻居，但两边的差相互抵消。所以结果为零，并不要求每个邻居都等于中心。
图中最后一组来自前面的鞍面 。在原点采样，左右值都是 ，上下值都是 ，四邻平均值为 ，恰好等于中心值。沿两个方向观察时，一条截线向上弯，另一条向下弯；改看四邻数值时，正负差仍然抵消。
这里的等式是有限间距 下，离散算子 的精确代数关系。对一般连续函数，不能把 与同样的四邻差直接画等号；只有当网格逐渐变细时，这个离散量才在光滑性条件下趋向连续拉普拉斯。
也可以回到曲面交互实验，换一个函数继续比较连续值与离散值。
7 热传导：由温差得到温度变化
设 是一块均匀材料的温度。梯度 指向温度升高的方向，而热量由高温处传向低温处。根据傅里叶定律，热流密度为：
其中 是导热系数。负号表示热流与温度梯度方向相反。
图11：中心80℃、四邻20℃时，温度梯度朝向中心，热流离开中心
接下来考虑一个小区域的热量收支。若流出的热量多于流入的热量，区域就会降温。因此，在没有内部热源时，能量守恒给出：
这里 是密度， 是比热容， 表示单位体积升高单位温度所需的热量。把热流表达式代入；如果材料均匀、各向同性，且 都为常数，就得到：
这就是热方程。其中 称为热扩散率。
注意这里有两个负号：热流相对于梯度有一个负号，温度变化相对于热流的净流出量又有一个负号。因此温度变化最终与 对应。在这个无源模型中， 表示升温， 表示降温。
手算一次更新
将空间上的拉普拉斯换成五点格式，再用一个很短的时间间隔 近似温度变化。记 为第 个时刻的中心温度，则：
整理后得到显式更新公式：
记 。若 用米、 用秒，则 的单位为平方米每秒， 没有单位。
例如取 、、，便有 。设当前中心为 ，四邻都是 ，一次更新后：
中心高于四邻平均，离散拉普拉斯为负，因此中心温度降低。用平均值 来写，同一步更新就是：
在这个例子中，中心向四邻平均温度移动了差值的 。更新整张网格时，每个新值都必须由同一时刻的旧值计算，不能一边覆盖旧数组，一边拿刚更新的值继续计算。
时间步不能任意取大
把更新式中的中心项合并：
对于这里使用的二维等间距五点格式与显式时间推进，通常要求：
在这个范围内，中心和四邻的权重都非负，而且总和为 。新温度是这五个旧温度的加权平均，不会一步越过它们的最大值或最小值。若把时间步取得过大，中心权重会变成负数，就可能产生不符合扩散直觉的振荡甚至数值发散。这个限制属于这套显式算法，换用其他时间推进方法时，条件也会不同。
下面可以涂抹初温、选择边界，使用播放、暂停或单步观察温度变化。数值参数以演示面板为准。
以下静态参考保留原文的固定参数。
下面再看一整块网格的变化。使用 个采样点，中央 个点初始为 ，其余为 ，四条边界始终固定在 。网格间距、热扩散率和时间步仍取上面的数值。
静态 SVG 参考｜图12：固定20℃边界下的二维热扩散，三个时刻共用20至80℃色标
热区逐渐扩散，最高温度下降。三幅图使用同一条色标，所以颜色可以直接比较。这里边界与外部的低温环境交换热量，材料内部的总热量会减少；不能把它当成绝热边界下的热量守恒例子。
8 调和：温度不再变化时
在前面的无源热方程中，如果温度已经不再随时间变化，就有 ，因而区域内部满足 。
更一般地，若函数在一个区域内处处满足：
就称它在这个区域内是调和函数（harmonic function）。只在某一点算出 ，还不能称整个函数在某个区域内调和。
对于五点离散格式，内部格子满足 ，就等价于：
也就是每个内部中心值都等于四邻平均值。这是一组同时成立的条件：一个格子的值会影响它的邻居，邻居又受更远的格子影响，最后还需要边界上的条件。
两端固定温度的一根杆
先看更简单的一维情况。沿均匀细杆等间距取五个点，两端温度始终固定为 和 ，假设杆内部没有热源，侧面也不散热。稳态时，三个内部点都满足：
把上式移项，得到 。这说明相邻两点之间的温差相等。两端总共相差 ，分成四段，每段就是 ，因此五个点依次为：
图13：固定两端温度时，杆上的线性稳态仍然存在热流
例如中间的 ，正好等于左右两点 和 的平均值。沿杆的温度梯度不为零，热量仍然从 的一端流向 的一端；只是每一小段的流入量与流出量相等，所以温度保持不变。
边界值决定了这根杆的稳态温度分布。如果两端都固定为 ，最终得到的就是处处为 的常数解。两端不同，得到的则是上图这样的线性分布。
连续函数的平均值性质
连续调和函数也有平均值性质。在二维中，如果一个闭圆盘完全位于调和区域内部，圆心值等于圆周上所有值的平均：
这里平均的是圆周上的一整圈，而前面的离散公式只取四个邻居。两者平均的对象不同，不能把连续调和函数的圆周平均性质直接换成任意有限间距的四邻平均等式。
前面的 就是一个调和函数，因为它在整个平面上都满足 。这个例子也说明，调和状态允许不同方向的变化相互抵消，并不要求函数值处处相同。
9 三维与矩形网格
三维的七点格式
三维情况下，拉普拉斯多了一个 方向的二阶导数。在三个方向的网格间距都为 时，沿用相同的差分方式：
合并后，就是六个轴向邻居之和减去 倍中心值，再除以 。它一共使用七个点，因此称为七点差分格式。若六邻平均值为 ，也可以写成：
图14：三维七点的系数与二维矩形网格的方向间距
横纵间距不相等时
如果二维网格是矩形，横向间距为 ，纵向间距为 ，就要分别使用各自的距离：
此时两个方向的系数不同，不能再直接用四邻的普通平均。将中心项合并，可以写成：
其中 是带权平均值：
左右两个邻居各取权重 ，上下两个邻居各取权重 ，再除以四个权重的总和。若 较小，左右邻居的系数就较大。当 时，这个带权平均又回到了前面的四邻普通平均。
上述公式都针对内部点。到了边界，还要按照给定的边界温度、法向热流等条件来处理。实际计算时，先确定各个量存放的位置、网格间距和边界条件，再按每个方向构造差分，就能得到所需的离散拉普拉斯。
参考
- John K. Hunter：Notes on Partial Differential Equations，第2章：拉普拉斯方程与调和函数的平均值性质。
- MIT：Heat Equation, Discrete and Continuous：由热流与能量守恒推导热方程。
附件：正文与原生 SVG 源文件
静态稿与原生 SVG 源文件备份：压缩包包含 Markdown 正文、14 幅 SVG 和审阅摘要。
离线 HTML 静态稿备份：公式、字体和 14 幅 SVG 已内嵌，可下载后离线打开。整页浏览器排版尚未完成实测，验证范围见压缩包中的审阅摘要。
Loading...