使用有限差分矩阵运算符的泊松方程的数值解

date

Apr 10, 2025

slug

datapossoi

status

Published

1. 引言

泊松方程是一种椭圆偏微分方程（PDE），在物理学和工程学的许多领域中无处不在。它的两个常见用途包括模拟静电势和重力势。泊松方程的某些形式也出现在流体流动 [ 1] 和传热问题 [ 2， 3] 中。通常，由于复杂的几何结构，对于实际感兴趣的情况，这些问题的解析解是不可能的。对泊松方程进行数值求解是研究这些系统的唯一方法。因此，泊松方程的数值求解方法具有重要意义 [ 2， 4， 5]。

有几种方法可以对泊松方程进行数值求解 [ 6]。有限差分法（FDM）是最简单和流行的方法之一 [ 7， 8， 9， 10]。该方法涉及将连续导数运算符替换为采用矩阵形式的近似离散有限差分运算符 [ 11]。该方法很直观，求解过程涉及成熟的线性代数方法。它存在一些缺点，例如准确表示复杂的几何图形和边界。有限元法（FEM）等方法更擅长处理任意几何图形。然而，由于 FDM 的易于实现和直观性，它仍然是一种非常流行的求解 PDE 的方法。

有几篇参考文献讨论了使用 FDM 求解泊松方程 [ 12]。虽然这些论文涵盖了许多细节，但它们通常没有提出为任意几何形状和任意边界条件构建数值系统（即系统矩阵形成）的系统方法。因此，对于 Poisson 方程的每种变体，线性代数系统都必须从头开始重建。本文的目标是介绍一种系统的方法，用于求解任何任意几何和边界条件的泊松方程，同时对数值求解器/计算机代码进行最少的修改。这将允许人们研究物理系统的各种配置，而无需重新制定数值配方。

设置 FDM 求解器有一些基本步骤。第一步是解空间的离散化。定义了空间中将求解方程的一系列网格点。在每一个点上，泊松方程都使用差分方程以离散形式写成。在传统方法中，通过修改边界点处的差分方程，在此步骤中加入边界条件。此步骤取决于几何结构和边界条件的类型。经过此调整后，我们剩下了一个线性方程组，从中可以形成全秩矩阵方程。此矩阵称为系统矩阵。然后对系统进行求解。请注意，当几何结构或边界条件发生变化时，系统矩阵会发生变化，因此必须重建。大多数参考文献没有为读者提供系统的方法来实现这一点。在许多情况下，需要研究同一问题的不同几何变化。因此，人们经常需要为每种配置修改其数值求解器的繁琐任务。

在本文中，系统矩阵定义为两个矩阵的组合。一个矩阵独立于问题几何，因此在所有情况下都保持不变。这称为导数运算符矩阵。对于 Poisson 方程，这是二阶导数运算符。另一个矩阵被定义为取决于几何和边界条件。这称为边界运算符矩阵。修改边界算子矩阵相对简单直观。结果表明，可以通过适当地编译这些矩阵中的行来构造系统矩阵。请注意，当几何和/或边界条件发生变化时，只有边界算子矩阵和系统矩阵构造过程发生变化。因此，可以通过对相同的数值代码进行细微调整来研究不同的几何变化（和/或边界条件）。这是所提出的方法的关键优势。

在本文中，详细介绍了如何实现灵活的泊松方程求解器的完整数值方法。此外，还讨论了该算法的 Python 实现，并通过 GitHub 存储库 [ 13] 提供给读者。此外，还提供了一些示例数值结果，这些结果展示了求解器的功能。

本文的其余部分组织如下：第 2 节定义了泊松方程和通常遇到的边界条件类型，第 3 节讨论了一维（1D）和二维（2D）情况的有限差分公式;第 4 节讨论了使用 Python 的计算机实现算法;第 5 节介绍了一些示例问题的数值结果;第 6 节中进行了总结性评论。

2. 泊松方程和边界条件

泊松方程由三个二阶空间导数组成。笛卡尔坐标中方程的一般形式为：

Δ𝑢(𝑥,𝑦,𝑧)=𝑓(𝑥,𝑦,𝑧),onΩ.

哪里

Δ=∇2=∂2∂𝑥2+∂2∂𝑦2+∂2∂𝑧2,

是笛卡尔坐标中的拉普拉斯算子，x、y、z 是独立的笛卡尔空间坐标， Ω 并且是解域（即要求解方程的 x、y、z 值的集合）。因变量 u 是要求解的函数，称为源函数，通常给出。通常，它是空间坐标的函数。在的特殊 𝑓=0 情况下，方程简化为拉普拉斯方程。

与任何微分方程一样，在尝试求解之前，必须定义边界条件。该问题的完整定义包括方程（ 1），其中定义了一个源函数，并在 u 或其一阶导数（即 ∂𝑢∂𝑥 、 ∂𝑢∂𝑦 和 ∂𝑢∂𝑧 ）上强制执行边界条件。边界条件也可以同时对 U 及其导数强制执行。解域 Ω 的边界表示为 ∂Ω 。边界可以细分为 M 个任意段或区域， ∂Ω1,∂Ω2,⋯,∂Ω𝑀 每个段或区域都有自己的边界条件。通常，会遇到四种常见的边界条件类型：

• 狄利克雷边界条件，其形式为： 𝑢|∂Ω𝑖=𝑢𝑖

• Neumann 边界条件，其形式为： ∂𝑢∂𝜂|∂Ω𝑖=𝑔𝑖

• 混合边界条件，形式为： 𝑢|∂Ω𝑖=𝑢𝑖 ， ∂𝑢∂𝜂|∂Ω𝑗=𝑔𝑗

• Robin 边界条件，其形式为： (𝐴𝑢+𝐵∂𝑢∂𝜂)|∂Ω𝑖=𝑔𝑖

这里 𝜂∈{𝑥,𝑦,𝑧} 表示笛卡尔空间变量之一， 𝑖,𝑗∈{1,2,⋯,𝑀} ， A， B 是常数，而 𝑢𝑖 是 𝑔𝑗 给定的量。请注意，对于诺依曼边界条件，由于约束施加在导数上，因此最多只能确定一个常数。因此，需要额外的约束才能获得唯一的解。此约束可以采用以下形式 [ 14， 15]：

边界区域也可以在仿真域本身内部指定，而不必在外周定义（即，其中一些 ∂Ω𝑖 ∂Ω1,∂Ω2,⋯,∂Ω𝑀 可能不一定位于的 Ω 外边缘）。这在狄利克雷边界中更为常见 [ 11]。这里介绍的解方案不区分内部和外部边界区域。该方法适用于 any arbitrary ∂Ω𝑖 for all 。

3有限差分离散化和矩阵运算符

方程（ 1）中的所有量都是连续变量。当以数值方式求解时，方程仅在离散点处求解。第一步是将变量和方程从连续域转换为离散域。为了简单起见，首先考虑 1D 情况。可以从一维情况构建更高维的情况。

3.1. 一维泊松方程

仅考虑一个空间变量 x，方程（ 1）可以写成：