2-39 一个二维稳态场，怎样用 5 点差分把拉普拉斯方程算成三维曲面？
MATLAB 二阶中心差分 + Jacobi 迭代 / 稀疏矩阵直接求解程序
一句话介绍： 这套 MATLAB 程序在 0～2 × 0～2 的二维区域上离散拉普拉斯方程，使用空间二阶中心差分构造 5 点差分格式，并分别通过固定次数 Jacobi 迭代和稀疏矩阵线性方程组求解内部节点，最后用 surf 输出二维解对应的三维曲面。
先看它解决什么问题
二维拉普拉斯方程描述的是一种典型稳态场问题：
∂²u/∂x² + ∂²u/∂y² = 0
它本身给出了区域内部应满足的关系，但并不会直接告诉我们每一个位置的数值。真正要得到可画出来的解，还需要先确定计算区域、网格和边界条件，再把连续偏微分方程转换成计算机能够处理的代数关系。
当前工程采用的主线非常清楚：
二维区域与边界条件 → 均匀网格 → 二阶中心差分 → 5 点离散 → 迭代求解或稀疏矩阵求解 → 三维曲面显示
这套程序的重点不是复杂物理建模，而是把“连续方程怎样变成数值解”完整展示出来。
1. 工程里真正有哪两条求解路线？
工程核心只有两个 MATLAB 脚本，两者求解的是同一类二维 Laplace 边值问题，但离散规模和求解方式不同。
文件 | 当前作用 |
|---|---|
Laplace_equation_2D.m | 60 × 60 网格，使用旧值数组 pn 做 10000 次 Jacobi 型显式迭代 |
LaplaceImplicit.m | 100 × 100 网格，组装稀疏系数矩阵后用 MATLAB 反斜杠运算直接求内部节点 |
license.txt | 工程随附的源码许可说明 |
两个脚本都不读取外部数据文件。输入直接写在脚本中，输出都是求解后的二维数值场及其三维曲面。
2. 当前程序求解的区域和边界条件是什么？
两个脚本都把 x、y 方向范围设置为 0～2。
迭代脚本使用：
nx = 60、ny = 60；
dx = 2/(nx-1)、dy = 2/(ny-1)。
稀疏矩阵脚本使用：
nx = 100、ny = 100；
同样由区域长度除以节点间隔数得到 dx、dy。
当前激活的边界条件可以概括为：
x = 0：u = 0
x = 2：u = y
y = 0：∂u/∂y = 0
y = 2：∂u/∂y = 0
也就是说，左右两侧采用 Dirichlet 边界，直接给定函数值；上下两侧采用零 Neumann 边界，要求法向梯度为 0。
因此右边界不是一个常数，而是随 y 从低到高变化。内部数值场要在这些约束之间形成平滑过渡。
3. 连续拉普拉斯方程怎样变成 5 点差分？
程序使用空间二阶中心差分近似二阶导数。
对一个内部网格点 u(i,j)，离散关系可以写成：
(u(i,j+1) - 2u(i,j) + u(i,j-1))/Δx² + (u(i+1,j) - 2u(i,j) + u(i-1,j))/Δy² = 0
当前两个方向都划分相同长度、相同节点数，所以 Δx = Δy。在这种等距网格下，中心点就等于四个直接相邻节点的平均：
u(i,j) = ¼·(u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1))
这就是程序所说的 5 点差分：一个中心点加上上、下、左、右四个邻点。
简单理解，内部每个点都不断被周围四个点“拉”到一个协调值，最终形成满足离散 Laplace 方程的稳态场。
4.Laplace_equation_2D.m 为什么属于 Jacobi 型迭代？
这个脚本先把 p 初始化为零，再施加四条边界。
随后每次循环都先执行：
pn = p
然后所有内部节点都只读取上一轮的 pn，一次性更新当前轮的 p。这正是 Jacobi 迭代最关键的特征：本轮所有新值都基于上一轮旧值，而不是立即使用刚刚更新的新值。
程序对内部区域 i = 2:ny-1、j = 2:nx-1 进行向量化更新，之后重新施加边界条件。
完整链路是：
初始化零场 → 写入边界 → 保存上一轮pn→ 更新全部内部节点 → 重施边界 → 重复 10000 次 → 绘图
它没有设置收敛误差阈值，而是固定执行 niter = 10000 次。因此这里的“结束条件”不是残差足够小，而只是迭代次数达到预设值。
5. 为什么每一轮都要重新写边界条件？
内部点更新时只应该由离散方程决定，而边界点属于外部给定约束。
如果边界在循环中也被普通差分更新，就会逐渐偏离原来的边值问题。当前脚本因此在每轮内部点计算完成后重新执行：
p(:,1)=0、p(:,nx)=y
并通过相邻行复 制 实现上下两侧零法向梯度：
p(1,:)=p(2,:)
p(ny,:)=p(ny-1,:)
对于零 Neumann 条件，这种离散处理相当于让边界值跟相邻内部值保持一致，从而使边界法向的一阶差分为零。
6.LaplaceImplicit.m 为什么不用 10000 次循环？
第二个脚本没有逐轮传播内部节点，而是把所有内部未知量一次写成一个大型线性方程组。
100 × 100 网格去掉四周边界后，内部区域为：
98 × 98 = 9604
因此程序最终求解的是 9604 个内部未知量。
它先分别构造 x、y 方向的二阶差分稀疏矩阵 Ax、Ay，再用 Kronecker 积组合成二维离散 Laplace 算子：
A = kron(Ay/dy², Iₓ) + kron(Iᵧ, Ax/dx²)
边界条件被整理到 bc 中，源项 S 当前初始化为全零，所以这里求解的仍然是纯 Laplace 方程，而不是带内部源项的 Poisson 方程。
最后通过：
S = A\S
一次求出全部内部节点，再 reshape 回二维数组。
源码注释把这条路线称为 implicit scheme。对于当前稳态 Laplace 问题，更准确地理解是：
它不是时间方向的隐式积分，而是把整个空间离散系统组装成稀疏线性方程组后直接求解。
7. 两种求解方式的核心差别是什么？
两个脚本使用的空间离散思想相同，真正不同的是“怎样把离散方程解出来”。
对比项 | Laplace_equation_2D.m | LaplaceImplicit.m |
|---|---|---|
网格 | 60 × 60 | 100 × 100 |
内部未知量 | 58 × 58 | 98 × 98 |
求解方式 | Jacobi 型固定次数迭代 | 稀疏线性系统直接求解 |
停止方式 | 固定 10000 次 | A\S 完成后结束 |
边界处理 | 每轮重新施加 | 先进入右端项，再回填边界 |
主要特点 | 过程直观，容易观察迭代思想 | 结构紧凑，直接利用 MATLAB 稀疏矩阵能力 |
所以这两个文件并不是两个不同物理模型，而是同一离散问题的两种数值求解路径。
8. 最终三维曲面表达的是什么？
迭代脚本最终调用：
surf(x,y,p,'EdgeColor','none')
稀疏矩阵脚本调用：
surf(x,y,u','EdgeColor','none')
后者之所以转置，是因为 u 在脚本中按 x 方向作为第一维存储，而 surf(x,y,Z) 的行列要与 y、x 坐标方向对应。
两者随后都使用 shading interp 去掉明显网格边线，让曲面显示得更连续。
需要注意：
shading interp只改变图形显示效果，不会增加数值网格分辨率，也不会提高差分求解精度。
从当前边界设置看，左侧边界固定为 0，右侧边界随 y 增大而增大，上下边界又保持零法向梯度，因此内部曲面会在这些条件之间形成连续、平滑的二维稳态分布。
9. 当前源码有哪些值得注意的数值边界？
这套程序结构很清楚，但如果准备在它的基础上继续扩展，有几个源码细节必须先知道。
迭代脚本没有自动收敛判断
Laplace_equation_2D.m 固定运行 10000 次，没有计算相邻两轮最大差值，也没有计算离散残差。
因此当网格、边界或区域改变后，10000 次不一定仍然是最合适的迭代次数。
当前显式更新公式依赖等距网格这一事实
源码当前设置 nx = ny = 60，而且 x、y 区间长度都为 2，因此 dx = dy，现有更新式与标准 5 点格式等价。
但代码中 dy² 与 i±1 方向、dx² 与 j±1 方向配对。由于 i 实际对应 y、j 对应 x，当以后让 dx ≠ dy 时，这一权重对应关系需要重新核对，不能直接把当前写法原样用于非等距网格。
角点最终服从后写入的 Neumann 处理
两个脚本都先写左右 Dirichlet 条件，再用相邻行复 制上下 Neumann 条件。
因此右上、右下角会被最后的 Neumann 赋值覆盖，不再严格等于 u = y 在两个端点的原始值。当前稀疏脚本的边界向量也按相邻 y 节点处理角点。
如果后续任务对角点值有严格定义，应单独确定角点到底服从哪一类边界。
10. 哪些参数最值得修改？
当前程序真正值得实验的主要是网格、迭代次数和边界，而不是绘图参数。
参数 / 设置 | 当前值 | 为什么值得改 | 调整后的主要影响 |
|---|---|---|---|
nx、ny | 60 × 60 或 100 × 100 | 决定空间离散分辨率 | 节点更多通常能更细致描述空间变化，但计算量和内存占用增加 |
niter | 10000 | 只影响迭代脚本的收敛程度与耗时 | 太少可能尚未稳定，过多则增加无必要计算 |
右边界 u = y | 线性随 y 变化 | 直接决定场在 x = 2 一侧的输入分布 | 修改后整个内部稳态场都会随之改变 |
上下边界类型 | 零 Neumann | 决定 y 方向边界是否允许法向梯度 | 改成 Dirichlet 或非零 Neumann 后，需要同步修改差分边界实现 |
如果要研究非等距网格，优先使用并核对 LaplaceImplicit.m 的矩阵离散关系，同时修正迭代脚本中 x、y 方向权重的对应关系。
11. 怎么运行这套程序？
运行 Laplace_equation_2D.m，即可执行 60 × 60 网格、10000 次 Jacobi 型迭代，并弹出三维求解曲面。
运行 LaplaceImplicit.m，即可建立 100 × 100 网格对应的稀疏离散系统，通过 A\S 求内部解，再输出三维曲面。
当前两个脚本都是自包含计算，不需要额外数据文件；所用的 sparse、kron、speye、反斜杠求解和 surf 均属于 MATLAB 基础数值与绘图能力。
如果只想快速理解有限差分的逐步传播过程，先看 Laplace_equation_2D.m；如果想理解二维差分算子怎样组装成全局稀疏线性系统，再看 LaplaceImplicit.m。
12. 一句话看懂这个项目
这是一个二维拉普拉斯方程有限差分求解项目：程序在 0～2 × 0～2 的矩形区域上设置左侧u=0、右侧u=y、上下零 Neumann 边界，用空间二阶中心差分把连续方程离散成 5 点格式；Laplace_equation_2D.m在 60 × 60 网格上通过旧值数组pn做 10000 次 Jacobi 型迭代，LaplaceImplicit.m则在 100 × 100 网格上把 9604 个内部未知量组装为稀疏线性系统后直接求解，最后都用surf输出三维稳态场曲面。当前实现适合作为二维 Laplace 方程、混合边界条件、5 点差分、迭代法和稀疏矩阵求解的数值计算示例；扩展到非等距网格时需要重新核对迭代式方向权重，角点边界优先级也应单独定义。
