静电学中，泊松方程的一般形式如下：

$$ \Delta \Phi=-\frac{\rho}{\varepsilon_0} $$

其中，$\varepsilon_0=8.854187817\times10^{12}\ \mathrm{F/m}$，称之为真空电导率，是一个近似值

<aside> 💡 同时，也 $\varepsilon_0$ 可以通过公式计算

$$ \varepsilon_0=\frac{1}{\mu_0c^2} $$

</aside>

其中 $c$ 是真空中的光速，$\mu_0$ 是真空磁导率

图中网格的边界为了不混淆，设 x 轴最大为 M，y 轴最大为 N

此时我们对泊松方程在 $x$ 与 $y$ 方向上分别离散，各个方向上的离散采用中心差分的方式。按照规则，当前的电子密度分布能计算出当前时刻的电势分布，故而这里各物理量的上标均为 $t$，为了推算简介，此处忽略。

$$ \begin{align} \nabla^2\Phi&=-\frac{\rho}{\varepsilon_{0}}\notag\\ \frac{\partial^2\Phi}{\partial x^2}+\frac{\partial^2\Phi}{\partial y^2} &= -\frac{\rho}{\varepsilon_{0}}\notag\\ \frac{\Phi_{\mathrm{N}} +\Phi_{\mathrm{S}} -2\Phi_{\mathrm{P}}}{\Delta x^2}+\frac{\Phi_{\mathrm{E}} +\Phi_{\mathrm{W}} -2\Phi_{\mathrm{P}}}{\Delta y^2}&= -\frac{\rho}{\varepsilon_{0}}\notag\\ \frac{\Phi_{\mathrm{N}} +\Phi_{\mathrm{S}} -2\Phi_{\mathrm{P}}}{h^2}+\frac{\Phi_{\mathrm{E}} +\Phi_{\mathrm{W}} -2\Phi_{\mathrm{P}}}{h^2}&= -\frac{\rho}{\varepsilon_{0}}\notag\\ \Phi_{\mathrm{N}} +\Phi_{\mathrm{S}} -2\Phi_{\mathrm{P}}+\Phi_{\mathrm{E}} +\Phi_{\mathrm{W}} -2\Phi_{\mathrm{P}}&= -\frac{h^2\rho}{\varepsilon_{0}}\notag\\ \Phi_{\mathrm{N}} +\Phi_{\mathrm{S}} -4\Phi_{\mathrm{P}}+\Phi_{\mathrm{E}} +\Phi_{\mathrm{W}} &= -\frac{h^2\rho}{\varepsilon_{0}}\notag\\ \Phi_{(i,j+1)} +\Phi_{(i,j-1)} -4\Phi_{(i,j)}+\Phi_{(i+1,j)} +\Phi_{(i-1,j)} &= -\frac{h^2\rho_{(i,j)}}{\varepsilon_{0}}\tag{1} \end{align} $$

现在，定义上述图 1 中网格一行之 $\Phi$ 所组成的向量。为了描述方便，直接默认矩阵是一个方阵，即 $\mathrm{N=M}$，且刨除掉首尾边界点

<aside> 🖊️ 节点下标如果是最大值为 $\mathrm N$，那么应当是 $\mathrm{0 \sim N-1}$

</aside>

$$ \Phi_j=(\Phi_{(1,j)},\Phi_{(2,j)},\cdots,\Phi_{(N-2,j)})^T,\ 1\leq j\leq \mathrm N-2\tag{2} $$

可以看到实际上五点差分只适合求解区域之内的点，如果是边界点就要用另外的方法。

根据式 (2) 的定义，实际上由 $\mathrm N-2$ 个 $\Phi_j$ 向量组成矩阵是一个 $\mathrm{(N-2)^2\times1}$ 的矩阵充当线性系统 $\mathrm Ax=b$ 中的 $x$。

<aside> 🖊️ $\mathrm{(N-2)\times(N-2)}$ 大小的标量矩阵压成一维向量

</aside>

而三对角矩阵必定是方阵，故而矩阵 A 的大小必定为 $(\mathrm N-2)^2\times(\mathrm N-2)^2$，进一步可以推敲出矩阵 b 的大小应当是 $\mathrm{(N-2)^2\times 1}$

那么同样，我们也可以定义上述图 1 中网格一行之 $\rho$ 所组成的向量

$$ \rho_j=(\rho_{(1,j)},\rho_{(2,j)},\cdots,\rho_{(N-2,j)})^T,\ 0\leq j\leq \mathrm N-2\tag{3} $$

<aside> 💡 这边定义向量实际上对应于图中的一行，但是在矩阵中还是按照列主序排列

</aside>

综合式 (2) (3)，能将式 (1) 改写为

$$ D\Phi_{j+1}+C\Phi_{j}+D\Phi_{j-1} = -\frac{h^2}{\varepsilon_{0}}\rho_j\ ,\ 1\leq j\leq \mathrm N-2\tag{4} $$

令 $-\frac{h^2}{\varepsilon_{0}}\rho_j=f_j$，同时可以将式 (4) 写成矩阵的形式

$$ \left(\begin{array}{ccccc}C & D & & & \\D & C & D & & \\& \ddots & \ddots & \ddots & \\& & D & C & D \\& & & D & C\end{array}\right){\mathrm{(N-2)^2\times (N-2)^2}}\left(\begin{array}{c}\boldsymbol{\Phi}1 \\\boldsymbol{\Phi}2 \\\vdots \\\boldsymbol{\Phi}{M-2} \\\boldsymbol{\Phi}{M-1}\end{array}\right){\mathrm{(N-2)^2}\times1}=\left(\begin{array}{c}\boldsymbol{f}_1-\boldsymbol{D} \boldsymbol{\Phi}0 \\\boldsymbol{f}2 \\\vdots \\\boldsymbol{f}{M-2} \\\boldsymbol{f}{M-1}-\boldsymbol{D} \boldsymbol{\Phi}M\end{array}\right){\mathrm{(N-2)^2}\times1}\tag{5} $$

其中， $C$ 与 $D$ 同样也是一个矩阵，实际上应该写成 $C_{\mathrm{(N-2)\times (N-2)}}$ 以及 $D_{\mathrm{(N-2)\times (N-2)}}$。且 C 是一个三对角矩阵，而 D 是一个对角阵。所以，到此位置我们就可以建立关于二维电势泊松方程的差分格式。

作为补充，可以看看矩阵 C 和 D 到底长什么样。

根据式 (1)，可以得到

$$ \begin{align} \Phi_{(i,j+1)} +\Phi_{(i,j-1)} -4\Phi_{(i,j)}+\Phi_{(i+1,j)} +\Phi_{(i-1,j)} &= -\frac{h^2\rho_{(i,j)}}{\varepsilon_{0}}\tag{1}\\ -\frac{\varepsilon_{0}}{h^2}\Phi_{(i,j+1)} -\frac{\varepsilon_{0}}{h^2}\Phi_{(i,j-1)} +\frac{4\varepsilon_{0}}{h^2}\Phi_{(i,j)}-\frac{\varepsilon_{0}}{h^2}\Phi_{(i+1,j)}-\frac{\varepsilon_{0}}{h^2}\Phi_{(i-1,j)} &= \rho_{(i,j)} \tag{6} \end{align} $$

根据对照关系，可以知道 C 矩阵作为三对角线矩阵，其主对角线附近的三元素分别为

$$ C_{\mathrm{elements}}=[-\frac{\varepsilon_{0}}{h^2},\frac{4\varepsilon_{0}}{h^2},-\frac{\varepsilon_{0}}{h^2}] $$

<aside> ❓ 为什么 $C$ 和 $D$ 都是一个 ${\mathrm{(N-2)\times (N-2)}}$ 的方阵？

根据矩阵乘法，就看式 (5) 的第一行，可以得到

$$ C\Phi_1+D\Phi_2=f_1-D\Phi_0 $$

首先，各个 $\Phi$ 的向量尺寸一致，都应该是 $\mathrm{(N-2)\times1}$

矩阵 $f$ 本身也是一个 $\mathrm{(N-2)\times1}$ 的向量

式子右侧的结果最终必定是 $\mathrm{(N-2)\times1}$，所以可以推出矩阵 D 是一个 ${\mathrm{(N-2)\times (N-2)}}$ 的方阵，同理可得 C 也是

</aside>

<aside> 💡 对于 $f$ 组成的向量，本身是只包含内点信息，但是通过计算之后加入了边界信息，即虽然一个 $f$ 是有 $\mathrm{M-2}$ 个元素，但实际上是完整的信息

</aside>

而矩阵 D 作为对角线矩阵，其主对角线上的元素为

$$ D_{\mathrm{elements}}=[-\frac{\varepsilon_{0}}{h^2}] $$

<aside> ❗ 上述描述的情况都是结构化网格在 $x$ 与 $y$ 方向是等距的情况

</aside>