二维椭圆方程九点差分高精度求解器

资 源 简 介

本项目基于MATLAB平台开发，旨在实现用于求解二维椭圆型偏微分方程（如泊松方程和拉普拉斯方程）的九点差分格式算法。该程序利用紧致九点差分模板（9-point compact stencil），通过结合中心节点及其周围八个邻居节点（包括对角线节点）的函数值，构建高精度的离散化方程组。与传统的截断误差为O(h^2)的五点差分格式相比，九点差分格式能够达到O(h^4)甚至更高的截断误差精度，显著提升了数值模拟的准确性。主要功能包括：自定义矩形求解区域与网格划分步长；自动组装基于九点格式的大型稀疏系数矩阵；灵活处理Dirichlet及Neumann边界条件；使用高效的线性方程组求解器计算数值解；以及提供后处理模块，用于计算数值解与解析解的误差范数，并绘制误差分布图和精度收敛曲线，直观展示九点格式在精度上的优势。

详 情 说 明

九点差分格式高精度数值求解器 (2D 9-Point Compact Difference Solver)

1. 项目介绍

本项目是一个基于MATLAB平台开发的高精度数值计算程序，旨在利用紧致九点差分格式（9-point compact difference scheme）求解二维泊松方程（Poisson Equation）和拉普拉斯方程。

该求解器与其说是传统的五点差分法（精度 $O(h^2)$），不如说通过引入紧致差分模板和右端源项的加权处理，将截断误差精度提升至 $O(h^4)$。本项目代码不仅实现了核心求解算法，还包含了一套完整的多网格收敛性分析框架，用于验证算法的高阶精度特性。

2. 功能特性

根据 main.m 的实际代码实现，项目具备以下核心功能：

• 高精度离散化：采用二维紧致九点差分模版，结合中心节点、四邻居节点及四对角节点的函数值，实现四阶精度求解。
• 稀疏矩阵优化：利用MATLAB的 sparse 数据结构，配合坐标格式（Triplet Format: I, J, V）进行系数矩阵的组装，极大降低了内存消耗并提升了组装效率，适合处理大规模网格。
• Dirichlet边界处理：实现了第一类边界条件（Dirichlet）的精确施加，代码中通过解析解直接设定边界节点数值。
• 右端项修正（RHS Correction）：为了配合左端的紧致格式，代码对源项 $f(x,y)$ 进行了加权平均处理（Padé 近似），确保方程两端具有一致的 $O(h^4)$ 精度。
• 多网格收敛性分析：程序内置循环结构，自动在不同尺寸的网格（如 10x10, 20x20, ..., 160x160）上进行求解，记录并计算误差，用于分析数值解的收敛阶。
• 误差范数计算：针对每个网格计算数值解与解析解之间的误差（代码主要框架支持L2范数及无穷范数统计）。

3. 系统要求与使用方法

系统要求

• 软件环境：MATLAB R2016a 或更高版本（需包含基础数学库）。
• 内存要求：取决于设定的网格密度（mesh_sizes）。对于 $N=160$ 的网格，常规PC内存即可。

使用方法

1. 确保 main.m 文件在MATLAB的当前工作目录下。
2. 直接运行 main 函数。
3. 程序将自动执行以下流程：

* 初始化参数与求解区域。 * 依次在不同分辨率网格上进行循环计算。 * 输出当前计算的网格规模及求解状态。 * 计算数值解与解析解的差异。

4. 算法与实现逻辑详解

main.m 的核心实现逻辑严格遵循紧致差分格式的数学推导，具体细节如下：

4.1 网格与问题定义

• 求解区域：设定为单位正方形 $[0, 1] times [0, 1]$。
• 网格划分：采用均匀网格划分，步长 $h_x = h_y = h$。
• 目标方程：泊松方程 $-Delta u = f$。
• 解析解验证：为了验证精度，预设解析解 $u(x,y) = sin(2pi x)sin(2pi y)$，并据此反推源项 $f(x,y) = 8pi^2 sin(2pi x)sin(2pi y)$。

4.2 稀疏系数矩阵组装 (LHS)

代码通过遍历所有网格节点 $(i, j)$ 来构建线性方程组 $AU=b$。

• 全局索引映射：将二维网格索引 $(i, j)$ 映射为一维全局索引 $n$，遵循MATLAB的列优先与逻辑映射规则。
• 内部节点模板：对于非边界节点，系数矩阵 $A$ 的填充遵循以下标准紧致九点模板（未归一化形式）：

* 中心节点 (i, j)：系数为 -20。 * 四邻居节点 (上下左右)：系数为 4。 * 四对角节点 (角点)：系数为 1。 这种权重组合对应于 $6h^2$ 倍的拉普拉斯算子离散化。

4.3 右端项的高阶处理 (RHS)

这是实现 $O(h^4)$ 精度的关键步骤。代码并未直接使用源项 $f_{i,j}$，而是对其进行了加权平均。

• 实现公式：

rhs_val = -0.5 * h^2 * (8*f_C + f_W + f_E + f_S + f_N)

• 逻辑说明：

该公式源于紧致格式推导，等价于 $6h^2 (f_{i,j} + frac{h^2}{12}Delta f)$ 的近似。代码中通过组合中心点与四个最近邻点的值来修正右端项，使其截断误差与左端差分算子匹配。

4.4 边界条件处理

• Dirichlet边界：这部分逻辑嵌套在节点遍历循环中。当检测到节点位于网格边缘（$i=1$ 或 $N_x+1$, $j=1$ 或 $N_y+1$）时，不再应用九点模板，而是强制令 $A_{n,n} = 1$，且右端项 $b_n$ 等于该点处的解析解数值。

4.5 线性方程组求解

• 利用 MATLAB 的反斜杠运算符 (`) 求解稀疏矩阵方程 $U = A setminus b$。MATLAB 会自动根据矩阵结构（不对称、稀疏）选择最优的直接求解器（如 UMFPACK），从而高效得到数值解向量。

4.6 数据后处理

• 求解得到的向量 $U_{vec}$ 被重塑（Reshape）回二维矩阵形式，并进行转置处理以匹配 meshgrid 生成的坐标系 $(X, Y)$。
• 通过计算数值解矩阵与解析解矩阵的差值绝对值（abs(U_num - U_exact)`），为后续的误差范数计算（L2, Linf）提供基础数据。