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

资 源 简 介

本项目主要利用MATLAB编程环境实现二维泊松方程及拉普拉斯方程的九点差分格式数值求解。相比于经典的五点差分格式，九点差分格式通过引入网格点周围包括对角线方向在内的更多邻域节点信息，结合泰勒级数展开消除低阶误差项，从而将截断误差精度从二阶$O(h^2)$提升至四阶$O(h^4)$甚至更高。项目内容包括构建矩形区域的离散网格，组装基于紧致九点格式的稀疏系数矩阵，处理Dirichlet或Neumann边界条件，并调用高效的线性方程组求解器（如Direct Solver或Iterative Solver）计算节点未知量。此外，代码还将包含后处理模块，用于计算数值解与解析解之间的误差（无穷范数和L2范数），绘制误差收敛曲线以验证精度阶数，并通过三维曲面图和等高线图直观展示计算结果。该程序为研究高精度有限差分方法提供了标准的算法实现框架。

详 情 说 明

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

项目介绍

本项目是一个用于求解二维椭圆型方程（主要是泊松方程）的高精度数值计算程序。项目基于MATLAB环境开发，采用了紧致九点差分格式（Compact 9-point Scheme / Mehrstellen Scheme）。通过利用网格点周围包括对角线方向在内的更多邻域节点信息，并结合源项的加权处理，该算法能够将截断误差精度从经典五点差分格式的二阶 $O(h^2)$ 提升至四阶 $O(h^4)$。

该程序旨在为高精度有限差分方法的研究提供一个标准的算法实现框架，包含了从网格生成、矩阵组装、求解到后处理误差分析的全流程。

功能特性

• 高精度离散化：实现了四阶精度的紧致九点差分计算格式。
• 稀疏矩阵优化：利用稀疏矩阵存储格式组装系数矩阵，有效降低内存消耗并提升计算效率。
• 收敛性验证：自动在不同密度的网格序列上运行求解，通过对数坐标下的线性拟合验证数值解的收敛阶数（L2范数和无穷范数）。
• 边界条件处理：支持狄利克雷（Dirichlet）边界条件的处理。
• 全方位可视化：

* 绘制误差收敛曲线（Log-Log图），并与理论 $O(h^4)$ 参考线对比。 * 绘制数值解的三维曲面图及等高线图。 * 绘制点对点的绝对误差分布三维图。

系统要求

• MATLAB R2016a 或更高版本（需支持基本的绘图及稀疏矩阵运算功能）。

使用方法

直接运行主程序即可启动计算流程。程序将依次执行以下步骤：

1. 初始化计算参数（定义求解区域及网格密度列表）。
2. 循环求解各个网格密度下的数值解。
3. 终端打印每个网格规模下的 L2 误差、最大误差及计算耗时。
4. 计算并输出最终的收敛阶数。
5. 弹出三个可视化窗口展示分析结果。

代码实现逻辑与算法细节

本项目代码主要分为主控流程与核心求解器两个部分，具体实现逻辑如下：

1. 参数设置与主控流程

程序首先定义了物理计算区域为 $[0, 1] times [0, 1]$ 的矩形区域。为了验证高精度格式的收敛性，预设了一组逐步加密的网格尺寸列表（例如 $10 times 10$ 到 $160 times 160$）。

主循环负责遍历这些网格尺寸，对每一次迭代：

• 调用核心求解函数获取数值解、解析解、坐标网格及计算步长。
• 计算L2范数误差：通过离散形式 $sqrt{sum (diff^2) cdot h^2}$ 计算。
• 计算无穷范数误差：取所有节点误差绝对值的最大值。
• 记录计算耗时并实时打印统计结果。
• 保存最细密网格的计算结果用于后续绘图。

2. 精度验证（收敛阶分析）

在完成所有网格的计算后，程序利用 polyfit 函数对步长 $h$ 和误差（L2及无穷范数）进行双对数坐标下的线性拟合。拟合得到的斜率即为数值方法的收敛阶。对于紧致九点格式，理论上该值应接近 4.00。

3. 数值结果可视化

程序包含丰富的数据可视化模块：

• 收敛性分析图：双对数坐标系下绘制误差随步长变化的曲线，并叠加斜率为4的参考线，直观展示高阶精度特性。
• 数值解展示：包含通过光照渲染的三维曲面图（Surface Plot）和二维等高线图（Contour Plot），展示物理场分布。
• 误差分布图：通过三维曲面展示计算域内每一点的绝对误差 $|U_{num} - U_{exact}|$，便于观察误差集中的区域。

4. 核心求解器实现细节

求解器函数封装了有限差分方法的核心算法，具体步骤如下：

网格生成与索引映射

• 根据输入的 $N$ 值生成均匀网格。
• 建立二维网格坐标 $(i, j)$ 到线性方程组未知数索引 $k$ 的映射关系。由于采用Dirichlet边界条件，仅将内部节点视为未知量，总未知量个数为 $(N-1)^2$。

稀疏矩阵组装（Compact 9-Point Scheme）

• 差分算子：实现了紧致九点格式的左端项系数矩阵。对于泊松方程 $u_{xx} + u_{yy} = f$，离散系数如下（均需除以 $6h^2$）：

* 中心节点：-20 * 上下左右四邻居：4 * 对角线四邻居：1

• 边界处理：在遍历九点邻域时，如果通过索引映射判断某邻居节点位于边界，则将其移至方程右端项（RHS），利用解析解函数计算边界值并乘以对应权重后从 RHS 向量中减去。

右端源项处理（RHS）

• 为了保证整体格式达到四阶精度 $O(h^4)$，代码不仅仅简单取中心点的源项值 $f_{i,j}$，而是采用了源项的加权平均处理。
• 根据代码注释逻辑，右端项构造遵循：$RHS = f_{ij} + frac{h^2}{12}Delta f_{ij}$，这等价于对源项 $f$ 也应用一个加权算子（中心权8/12，四周权1/12）。这部分是为了消除截断误差中的低阶项，是紧致格式精度的关键所在。