6.3 应用: Poisson 方程求解

Poisson方程是一类典型的偏微分方程,经过差分离散后,转化为一个线性方程组.我们以这个线性方程组为例,讨论Jacobi,G-S和SOR迭代法的相关性质.

6.3.1 一维 Poisson 方程

首先介绍一维 Poisson 方程的离散. 考虑如下带 Dirichlet 边界条件的一维 Poisson 方程

\left\{ \begin{array}{l} - \frac {\mathrm {d} ^ {2} u (x)}{\mathrm {d} x ^ {2}} = f (x), \quad 0 < x < 1, \\ u (0) = a, u (1) = b, \end{array} \right. \tag {6.16}

其中 $f(x)$ 是给定的函数， $u(x)$ 是需要计算的未知函数

差分离散

取步长 $h = \frac{1}{n + 1}$ ，网格节点为 $x_{i} = ih,i = 0,1,2,\ldots ,n + 1.$ 我们采用二阶中心差分来近似二阶导数，即

- \frac {\mathrm {d} ^ {2} u (x)}{\mathrm {d} x ^ {2}} \Bigg | _ {x _ {i}} = \frac {2 u (x _ {i}) - u (x _ {i - 1}) - u (x _ {i + 1})}{h ^ {2}} + O \left(h ^ {2} \left\| \frac {\mathrm {d} ^ {4} u}{\mathrm {d} x ^ {4}} \right\| _ {\infty}\right), i = 1, 2, \ldots , n.

将其代入 (6.9), 舍去高阶项后得到 Poisson 方程在 $x_{i}$ 点的近似离散方程

- u _ {i - 1} + 2 u _ {i} - u _ {i + 1} = h ^ {2} f _ {i},

其中 $f_{i} = f(x_{i}),u_{i}$ 为 $u(x_{i})$ 的近似.令 $i = 1,2,\ldots ,n,$ 则可得 $n$ 个线性方程.写成矩阵形式为

T _ {n} u = f, \tag {6.17}

其中

T _ {n} = \left[ \begin{array}{c c c c} 2 & - 1 & & \\ - 1 & \ddots & \ddots & \\ & \ddots & \ddots & - 1 \\ & & - 1 & 2 \end{array} \right] \triangleq \operatorname {t r i d i a g} (- 1, 2, - 1), \quad u = \left[ \begin{array}{c} u _ {1} \\ u _ {2} \\ \vdots \\ u _ {n - 1} \\ u _ {n} \end{array} \right], \quad f = \left[ \begin{array}{c} f _ {1} + u _ {0} \\ f _ {2} \\ \vdots \\ f _ {n - 1} \\ f _ {n} + u _ {n + 1} \end{array} \right]. \tag {6.18}

由边界条件可知 $u_0 = a, u_{n+1} = b$ .

系数矩阵 $T_{n}$ 的性质

易知 $T_{n}$ 是三对角对称矩阵，因此其特征值都是实数.事实上，我们有下面的结论

引理6.24 $T_{n}$ 的特征值和对应的特征向量分别为

\lambda_ {k} = 2 - 2 \cos \frac {k \pi}{n + 1},

z _ {k} = \sqrt {\frac {2}{n + 1}} \left[ \sin \frac {k \pi}{n + 1}, \sin \frac {2 k \pi}{n + 1}, \dots , \sin \frac {n k \pi}{n + 1} \right] ^ {\mathsf {T}}, \quad k = 1, 2, \dots , n,

即 $T_{n} = Z\Lambda Z^{\top}$ ，其中 $\Lambda = \mathrm{diag}(\lambda_1,\lambda_2,\ldots ,\lambda_n)$ 是对角矩阵， $Z = [z_{1},z_{2},\dots ,z_{n}]$ 是正交矩阵.

(留作课外自习, 直接代入验证即可)

由此可知, $T_{n}$ 是对称正定的. 更一般地, 我们有下面的结论

推论6.25 设 $T = \mathrm{tridiag}(a,b,c) \in \mathbb{R}^{n \times n}, ac > 0$ ，则 $T$ 的特征值为

\lambda_ {k} = b - 2 \sqrt {a c} \cos \frac {k \pi}{n + 1}, \quad k = 1, 2, \ldots , n,

对应的特征向量为 $z_{k}$ ，其第 $j$ 个分量为

z _ {k} (j) = \left(\frac {a}{c}\right) ^ {\frac {j}{2}} \sin \frac {j k \pi}{n + 1}.

特别地, 若 $a = c = 1$ , 则对应的单位特征向量为

z _ {k} = \sqrt {\frac {2}{n + 1}} \left[ \sin \frac {k \pi}{n + 1}, \sin \frac {2 k \pi}{n + 1}, \ldots , \sin \frac {n k \pi}{n + 1} \right] ^ {\mathsf {T}}.

(留作练习)

由引理6.24可知， $T_{n}$ 的最大特征值为

2 \left(1 - \cos \frac {n \pi}{n + 1}\right) = 4 \sin^ {2} \frac {n \pi}{2 (n + 1)} \approx 4,

最小特征值为

2 \left(1 - \cos \frac {\pi}{n + 1}\right) = 4 \sin^ {2} \frac {\pi}{2 (n + 1)} \approx \left(\frac {\pi}{n + 1}\right) ^ {2}.

因此, 当 $n$ 很大时, $T_{n}$ 的谱条件数约为

\kappa_ {2} (T _ {n}) \approx \frac {4 (n + 1) ^ {2}}{\pi^ {2}}.

6.3.2 二维Poisson方程

考虑二维 Poisson 方程的离散. 同样是带 Dirichlet 边界条件, 即

\left\{ \begin{array}{l} - \Delta u (x, y) = - \frac {\partial^ {2} u (x , y)}{\partial x ^ {2}} - \frac {\partial^ {2} u (x , y)}{\partial y ^ {2}} = f (x, y), \quad (x, y) \in \Omega , \\ u (x, y) = u _ {0} (x, y), \quad (x, y) \in \partial \Omega , \end{array} \right. \tag {6.19}

其中 $\Omega = [0,1]\times [0,1]$ 为求解区域， $\partial \Omega$ 表示 $\Omega$ 的边界

五点差分离散

为了简单起见, 我们在 $x$ -方向和 $y$ -方向取相同的步长 $h = \frac{1}{n + 1}$ , 网格节点为 $x_{i} = ih$ , $y_{j} = jh$ , $i, j = 0, 1, 2, \ldots, n + 1$ . 在 $x$ -方向和 $y$ -方向同时采用二阶中心差分近似, 可得

\begin{array}{l} \left. - \frac {\partial^ {2} u (x , y)}{\partial x ^ {2}} \right| _ {(x _ {i}, y _ {j})} \approx \frac {2 u (x _ {i} , y _ {j}) - u (x _ {i - 1} , y _ {j}) - u (x _ {i + 1} , y _ {j})}{h ^ {2}}, \\ - \frac {\partial^ {2} u (x , y)}{\partial y ^ {2}} \bigg | _ {(x _ {i}, y _ {j})} \approx \frac {2 u (x _ {i} , y _ {j}) - u (x _ {i} , y _ {j - 1}) - u (x _ {i} , y _ {j + 1})}{h ^ {2}}. \\ \end{array}

代入 (6.12) 后, 就可以得到二维 Poisson 方程在 $(x_{i}, y_{j})$ 点的近似离散方程

4 u _ {i, j} - u _ {i - 1, j} - u _ {i + 1, j} - u _ {i, j - 1} - u _ {i, j + 1} = h ^ {2} f _ {i, j},

其中 $f_{i,j} = f(x_i,y_j),u_{i,j}$ 为 $u(x_{i},y_{j})$ 的近似.这个离散格式可以用下面的图(左图)来描述

按自然顺序排列以上方程 (即先沿 $x$ -方向从左到右, 然后沿 $y$ 方向从下往上), 可得矩阵形式

T _ {N} u = h ^ {2} f, \tag {6.20}

其中

T _ {N} \triangleq I \otimes T _ {n} + T _ {n} \otimes I, \quad u = \left[ u _ {1, 1}, \dots , u _ {n, 1}, u _ {1, 2}, \dots , u _ {n, 2}, \dots , u _ {1, n}, \dots , u _ {n, n} \right] ^ {\top}.

这里 $\otimes$ 表示 Kronecker 乘积, $T_{n}$ 为一维 Poisson 方程离散后的系数矩阵, 即 (6.11). 上图 (右图) 为 $N = 7^{2}$ 时的矩阵示例图.

需要注意的是, 系数矩阵 $T_{N}$ 与网格点的排序有关, 不同的排序方式会得到不同的 $T_{N}$ .
类似地, 如果对三维 Poisson 方程进行中心差分离散, 则对应的系数矩阵为

T _ {n} \otimes I \otimes I + I \otimes T _ {n} \otimes I + I \otimes I \otimes T _ {n}.

系数矩阵 $T_{N}$ 的性质

定理6.26 设 $T_{n} = Z\Lambda Z^{\top}$ , 其中 $Z = [z_{1}, z_{2}, \ldots, z_{n}]$ 为正交阵, $\Lambda = \mathrm{diag}(\lambda_{1}, \lambda_{2}, \ldots, \lambda_{n})$ 为对角阵, 则 $T_{N}$ 的特征值分解为

T _ {N} = (Z \otimes Z) (I \otimes \Lambda + \Lambda \otimes I) (Z \otimes Z) ^ {\mathsf {T}},

即 $T_{N}$ 的特征值为 $\lambda_{i} + \lambda_{j}$ , 对应的特征向量为 $z_{i} \otimes z_{j}, i, j = 1,2,\ldots,n$ .

由于 $T_{N}$ 对称正定，所以其条件数为

\kappa (T _ {N}) = \frac {\lambda_ {\max} (T _ {N})}{\lambda_ {\min} (T _ {N})} = \frac {1 - \cos \frac {n \pi}{n + 1}}{1 - \cos \frac {\pi}{n + 1}} = \frac {\sin^ {2} \frac {n \pi}{2 (n + 1)}}{\sin^ {2} \frac {\pi}{2 (n + 1)}} \approx \frac {4 (n + 1) ^ {2}}{\pi^ {2}}.

故当 $n$ 越来越大时， $\kappa(T_N)$ 也越来越大，即 $T_N$ 越来越病态。

求解二维离散 Poisson 方程的 Jacobi 迭代法

算法6.5. 求解二维离散Poisson方程的Jacobi迭代法
1: Given an initial guess $u^{(0)} \in \mathbb{R}^N$
2: while not converge do
3: for $i = 1$ to $n$ do
4: for $j = 1$ to $n$ do
5: $u_{i,j}^{(k + 1)} = \frac{1}{4}\left(h^2 f_{i,j} + u_{i + 1,j}^{(k)} + u_{i - 1,j}^{(k)} + u_{i,j + 1}^{(k)} + u_{i,j - 1}^{(k)}\right)$
6: end for
7: end for
8: end while

求解二维离散 Poisson 方程的 G-S 迭代法

我们可以直接用Gauss-Seidel迭代法对基于自然顺序的离散Poisson方程(6.13)进行数值求解，即

算法6.6.求解二维离散Poisson方程的红黑排序G-S迭代法
1: Given an initial guess $u^{(0)} \in \mathbb{R}^N$
2: while not converge do
3: for $i = 1$ to $n$ do
4: for $j = 1$ to $n$ do
5: $u_{i,j}^{(k + 1)} = \frac{1}{4}\left(h^2 f_{i,j} + u_{i + 1,j}^{(k)} + u_{i - 1,j}^{(k + 1)} + u_{i,j + 1}^{(k)} + u_{i,j - 1}^{(k + 1)}\right)$
6: end for
7: end for
8: end while

但此时迭代解的更新必须按自然顺序进行 (即更新 $u_{i,j}^{(k + 1)}$ 时必须先更新 $u_{i - 1,j}^{(k + 1)}$ 和 $u_{i,j - 1}^{(k + 1)}$ 的值), 因此不适合并行计算.

下面我们介绍一种适合并行计算的排序方法：红黑排序，即将二维网格点依次做红黑记号，如右图所示.

在计算过程中, 对 $u_{i,j}^{(k + 1)}$ 的值进行更新时, 我们可以先更新红色节点上的值, 此时所使用的只是黑色节点的数据, 然后再更新黑色节点上的值, 这时使用的是红色节点的数据. 由于在更新红点时, 各个红点之间是相互独立的, 因此可以并行计算. 同样, 在更新黑点时, 各个黑点之间也是相互独立的, 因此也可以并行计算.

1: Given an initial guess $u^{(0)} \in \mathbb{R}^N$
2: while not converge do

3: for $(i,j)$ 为红色节点 do
4: $u_{i,j}^{(k + 1)} = \frac{1}{4}\left(h^2 f_{i,j} + u_{i + 1,j}^{(k)} + u_{i - 1,j}^{(k)} + u_{i,j + 1}^{(k)} + u_{i,j - 1}^{(k)}\right)$ % 使用旧数据
5: end for
6: for $(i,j)$ 为黑色节点 do
7: $u_{i,j}^{(k + 1)} = \frac{1}{4}\left(h^2 f_{i,j} + u_{i + 1,j}^{(k + 1)} + u_{i - 1,j}^{(k + 1)} + u_{i,j + 1}^{(k + 1)} + u_{i,j - 1}^{(k + 1)}\right)$ % 使用新数据
8: end for
9: end while

我们记基于红黑排序所得到的系数矩阵为 $T_N^{\mathrm{RB}}$ ，它可以写成 $2\times 2$ 分块形式，且两个对角块都是对角矩阵，如下图所示 $(N = 7^{2})$

事实上, $T_{N}^{\mathrm{RB}}$ 可以通过对 $T_{N}$ 进行行置换和列置换得到, 即存在置换矩阵 $P$ 使得

T _ {N} ^ {\mathrm {R B}} = P T _ {N} P ^ {\intercal}.

求解二维离散 Poisson 方程的 SOR 迭代法

算法6.7. 求解二维离散Poisson方程的SOR迭代法

1: Given an initial guess $u^{(0)} \in \mathbb{R}^N$ and a parameter $\omega$
2: while not converge do
3: for $i = 1$ to $n$ do
4: for $j = 1$ to $n$ do
5: $u_{i,j}^{(k + 1)} = (1 - \omega)u_{i,j}^{(k)} + \frac{\omega}{4}\left(h^2 f_{i,j} + u_{i + 1,j}^{(k)} + u_{i - 1,j}^{(k + 1)} + u_{i,j + 1}^{(k)} + u_{i,j - 1}^{(k + 1)}\right)$
6: end for
7: end for

8: end while

A上面的SOR算法是基于自然排序,基于红黑排序的SOR算法请读者自己完成.

例6.2 已知二维Poisson方程

\left\{ \begin{array}{l l} - \Delta u (x, y) = - 1, & (x, y) \in \Omega \\ u (x, y) = \frac {x ^ {2} + y ^ {2}}{4}, & (x, y) \in \partial \Omega \end{array} \right.

其中 $\Omega = (0,1)\times (0,1)$ . 该方程的解析解是 $u(x,y) = \frac{x^2 + y^2}{4}$ . 用五点差分格式离散后得到一个线性方程组.

(1) 分别用 Jacobi, G-S 和 SOR 迭代法计算该方程组的解.
(2) 分别用SOR和SSOR迭代法求解该方程组, 观察参数 $\omega$ 对算法收敛的影响.

解．(1)MATLAB程序参见Poisson_Jacobi_GS_SOR.m.定义近似解的相对误差：

\operatorname {r e l e r r} _ {k} \triangleq \frac {\left\| u ^ {(k)} - u _ {*} \right\| _ {2}}{u _ {*}},

其中 $u_{*}$ 表示精确解. 下图画出了 $n = 16$ (即 $N = 256$ ) 时三种方法的近似解相对误差的下降过程 (前 100 个迭代步).

图6.1.Jacobi,G-S和SOR的近似解相对误差的下降曲线

(2) MATLAB 程序参见 Poisson_SOR_omega.m 和 Poisson_SSOR_omega.m. 下图中画出了 $n = 8$ (即 $N = 64$ ) 时, SOR 和 SSOR 收敛结果与参数 $\omega$ 取值之间的关系.

图6.2.SOR和SSOR的收敛结果与 $\omega$ 取值的关系

6.3.3 收敛性分析

对于基本迭代法, 其收敛的充要条件是迭代矩阵的谱半径小于 1. 当谱半径不可求时, 我们可以根据迭代矩阵的范数来判断, 即如果迭代矩阵的某个范数小于 1, 则迭代法收敛.

本小节考虑求解二维离散 Poisson 方程的 Jacobi, G-S 和 SOR 的收敛性. 考虑这些方法的收敛性, 只需研究相应的迭代矩阵的谱半径即可. 对于二维离散 Poisson 方程, 系数矩阵为

A = T = I \otimes T _ {n} + T _ {n} \otimes I.

故 Jacobi 迭代法的迭代矩阵为

G _ {\mathrm {J}} = D ^ {- 1} (L + U) = (4 I) ^ {- 1} (4 I - T) = I - T / 4. \tag {6.21}

由于 $T$ 的特征值为

\lambda_ {i} + \lambda_ {j} = 2 \left(1 - \cos \frac {\pi i}{n + 1}\right) + 2 \left(1 - \cos \frac {\pi j}{n + 1}\right) = 4 - 2 \left(\cos \frac {\pi i}{n + 1} + \cos \frac {\pi j}{n + 1}\right),

所以 $G_{\mathrm{J}}$ 的特征值为

1 - \left(\lambda_ {i} + \lambda_ {j}\right) / 4 = \frac {1}{2} \left(\cos \frac {\pi i}{n + 1} + \cos \frac {\pi j}{n + 1}\right).

故

\rho (G _ {\mathrm {J}}) = \frac {1}{2} \max _ {i, j} \left\{\left| \cos \frac {\pi i}{n + 1} + \cos \frac {\pi j}{n + 1} \right| \right\} = \cos \frac {\pi}{n + 1} < 1,

即 Jacobi 迭代法是收敛的.

注意当 $n$ 越来越大时, $\kappa (T)\to \infty$ , 即 $T$ 越来越病态, 此时 $\rho (G_{\mathrm{J}})\rightarrow 1$ , 即 Jacobi 迭代法收敛越来越慢.

通常, 问题越病态就越难求解, 舍入误差对解的影响也越大.

关于G-S和SOR,我们有下面的结论

定理6.27 设 $G_{\mathrm{GS}}$ 和 $G_{\mathrm{SOR}}$ 分别表示求解基于红黑排序的二维Poisson方程的G-S和SOR的迭代矩阵, 则有

\rho (G _ {\mathrm {G S}}) = \rho (G _ {\mathrm {J}}) ^ {2} = \cos^ {2} \frac {\pi}{n + 1} < 1

\rho (G _ {\mathrm {S O R}}) = \frac {\cos^ {2} \frac {\pi}{n + 1}}{\left(1 + \sin \frac {\pi}{n + 1}\right) ^ {2}} < 1, \quad \omega = \frac {2}{1 + \sin \frac {\pi}{n + 1}}.

证明. 由于此时的系数矩阵 $T_{N}$ 具有性质A, 因此直接利用定理6.23即可.

在上述结论中, SOR 中的 $\omega$ 是最优参数, 即此时的 $\rho(G_{\mathrm{SOR}})$ 最小. 由 Taylor 公式可知, 当 $n$ 很大时, 有

\rho (G _ {\mathrm {J}}) = \cos \frac {\pi}{n + 1} \approx 1 - \frac {\pi^ {2}}{2 (n + 1) ^ {2}} = 1 - O \left(\frac {1}{n ^ {2}}\right),

\rho (G _ {\mathrm {S O R}}) = \frac {\cos^ {2} \frac {\pi}{n + 1}}{\left(1 + \sin \frac {\pi}{n + 1}\right) ^ {2}} \approx 1 - \frac {2 \pi}{n + 1} = 1 - O \left(\frac {1}{n}\right).

由于当 $n$ 很大时有

\left(1 - \frac {1}{n}\right) ^ {k} \approx 1 - \frac {k}{n} = 1 - \frac {k n}{n ^ {2}} \approx \left(1 - \frac {1}{n ^ {2}}\right) ^ {k n},

即SOR迭代 $k$ 步后误差的减小量与Jacobi迭代 $kn$ 步后误差减小量差不多.因此,对于二维离散Poisson方程，当SOR取最优参数时，收敛速度大约是Jacobi的 $n$ 倍

需要指出的是, 对于一般线性方程组, 上述结论不一定成立.

由于 $\rho (G_{\mathrm{GS}}) = \rho (G_{\mathrm{J}})^2$ ，因此,对于二维离散Poisson方程，G-S的收敛速度大约是Jacobi的2倍.

事实上, 当 $n$ 很大时, 这三个方法的收敛速度都很慢.

例6.3同例6.2,分别对不同的 $n$ ，观察Jacobi,Gauss-Seidel和SOR迭代法的收敛情况

解. 参见 MATLAB 程序 Poisson_Jacobi_GS_SOR_convergence.m. 下图中画出了 $N = 16,32,64,128$ 时, 这三种方法的相对误差下降情况.

事实上, 由于二维离散 Poisson 方程的系数矩阵 $T_{N}$ 是不可约对角占优的, 因此根据定理 6.8, Jacobi 和 G-S 都收敛. 另外, $T_{N}$ 是对称正定的, 故当 $0 < \omega < 2$ 时, SOR 迭代法收敛.

6.3.4 快速求解方法

如果已经知道矩阵 $A$ 的特征值分解 $A = X\Lambda X^{-1}$ , 则 $Ax = b$ 的解可表示为

x = A ^ {- 1} b = X \Lambda^ {- 1} X ^ {- 1} b.

如果 $A$ 是正规矩阵, 即 $X$ 是酉矩阵, 则

x = A ^ {- 1} b = X \Lambda^ {- 1} X ^ {*} b.

一般来说, 我们不会采用这种特征值分解的方法来解线性方程组, 因为计算特征值分解通常比解线性方程组更困难. 但在某些特殊情况下, 我们可以由此得到快速方法.

考虑二维离散 Poisson 方程

T u = h ^ {2} f, \tag {6.22}

其中

T = I \otimes T _ {n} + T _ {n} \otimes I, \quad T _ {n} = \operatorname {t r i d i a g} (- 1, 2, - 1) \in \mathbb {R} ^ {n \times n}.

由定理6.26可知

T = (Z \otimes Z) (I \otimes \Lambda + \Lambda \otimes I) (Z \otimes Z) ^ {\mathsf {T}},

其中 $Z = [z_{1},z_{2},\ldots ,z_{n}]$ 是正交矩阵.这里

z _ {k} = \sqrt {\frac {2}{n + 1}} \left[ \sin \frac {k \pi}{n + 1}, \sin \frac {2 k \pi}{n + 1}, \ldots , \sin \frac {n k \pi}{n + 1} \right] ^ {\mathsf {T}}, \quad k = 1, 2, \ldots , n.

所以，方程(6.15)的解为

u = T ^ {- 1} h ^ {2} f = \left[ (Z \otimes Z) (I \otimes \Lambda + \Lambda \otimes I) ^ {- 1} (Z \otimes Z) ^ {\intercal} \right] h ^ {2} f.

因此, 主要的运算是 $Z \otimes Z$ 与向量的乘积, 以及 $(Z \otimes Z)^{\mathsf{T}}$ 与向量的乘积. 而这些乘积可以通过快速 Sine 变换来实现.

快速Fourier变换

快速Fourier变换(FFT)是用来计算离散Fourier变换(DFT)矩阵与向量乘积的一种快速方法.

设 $x = [x_0, x_1 \ldots, x_{n-1}]^\top \in \mathbb{C}^n$ , 其DFT定义为 $y = \mathrm{DFT}(x) = [y_0, y_1 \ldots, y_{n-1}]^\top \in \mathbb{C}^n$ , 其中

y _ {k} = \sum_ {j = 0} ^ {n - 1} \omega_ {n} ^ {k j} x _ {j}, \quad k = 0, 1, \ldots , n - 1.

这里

\omega_ {n} = e ^ {- \frac {2 \pi \mathbf {i}}{n}} = \cos (2 \pi / n) - \mathbf {i} \sin (2 \pi / n)

是1的一个 $n$ 次本原根(primitive $n$ -th root of unity), $\mathbf{i}$ 是虚部单位.

这里说 $\omega_{n}$ 是primitive $n$ -th root of unity是指 $\omega_{n}^{n} = 1$ 且

\omega_ {n} ^ {k} \neq 1, \quad k = 1, 2, \dots , n - 1.

构造矩阵 $F_{n} = [f_{kj}]\in \mathbb{C}^{n\times n}$ ，其中 $f_{kj} = \omega_n^{kj} = e^{-\frac{2kj\pi\mathbf{i}}{n}} = \cos (2kj\pi /n) - \mathbf{i}\sin (2kj\pi /n)$ ，即

F _ {n} = \left[ \begin{array}{c c c c c} 1 & 1 & 1 & \dots & 1 \\ 1 & \omega_ {n} & \omega_ {n} ^ {2} & \dots & \omega_ {n} ^ {n - 1} \\ 1 & \omega_ {n} ^ {2} & \omega_ {n} ^ {4} & \dots & \omega_ {n} ^ {2 (n - 1)} \\ \vdots & \vdots & & \ddots & \vdots \\ 1 & \omega_ {n} ^ {n - 1} & \omega_ {n} ^ {2 (n - 1)} & \dots & \omega_ {n} ^ {(n - 1) ^ {2}} \end{array} \right],

则有

y = \mathrm {D F T} (x) = F _ {n} x.

我们称矩阵 $F_{n}$ 为DFT矩阵. 易知DFT矩阵具有以下性质:

(1) $F_{n}$ 是对称矩阵 (但不是 Hermitian);
(2) $F_{n}^{*}F_{n} = nI$ ，所以 $\frac{1}{\sqrt{n}} F_{n}$ 是酉矩阵

相应的离散Fourier反变换定义为 $x = \mathrm{IDFT}(y)$ ，其中

x _ {j} = \frac {1}{n} \sum_ {k = 0} ^ {n - 1} \omega_ {n} ^ {- j k} y _ {k}, \quad j = 0, 1, \ldots , n - 1.

写成矩阵形式为

x = \frac {1}{n} F ^ {*} y.

DFT和IDFT满足下面的性质：

\mathrm {I D F T} (\mathrm {D F T} (x)) = x,

\mathrm {D F T} (\mathrm {I D F T} (y)) = y.

在MATLAB中，计算DFT和IDFT的函数分别为fft和ifft,即： $y = \mathrm{fft}(x),x = \mathrm{ifft}(y)$ (测试代码见FFT_test.m)

离散Sine变换

离散Sine变换(DST)有多种定义, 我们这里只介绍与求解 Poisson 方程有关的一种定义, 其它定义可参见[17].

设 $x = [x_{1},x_{2},\ldots ,x_{n}]^{\mathsf{T}}\in \mathbb{R}^{n}$ ，其离散Sine变换定义为 $y = \mathrm{DST}(x) = [y_1,y_2,\dots ,y_n]^{\mathsf{T}}\in \mathbb{R}^{n}$ 其中

y _ {k} = \sum_ {j = 1} ^ {n} x _ {j} \sin \left(\frac {k j \pi}{n + 1}\right), \quad k = 1, 2, \ldots , n.

对应的离散Sine反变换记为IDST,即 $x = \mathrm{IDST}(\mathrm{y})$ ，其中

x _ {j} = \frac {2}{n + 1} \sum_ {k = 1} ^ {n} y _ {k} \sin \left(\frac {j k \pi}{n + 1}\right), \quad j = 1, 2, \dots , n.

DST 和 IDST 满足下面的性质:

\operatorname {I D S T} (\operatorname {D S T} (x)) = x,

\mathrm {D S T} (\mathrm {I D S T} (y)) = y.

在MATLAB中，计算DST和IDST的函数分别为dst和idst,即： $y = \mathrm{dst}(x),x = \mathrm{idst}(y)$ (测试代码见DST_test.m)

Possion方程与DST

我们首先考虑矩阵 $Z$ 与一个任意给定向量 $b$ 的乘积.设 $y = Zb$ ，则

y _ {k} = \sum_ {j = 1} ^ {n} Z (k, j) b _ {j} = \sqrt {\frac {2}{n + 1}} \sum_ {j = 1} ^ {n} b _ {j} \sin \left(\frac {k j \pi}{n + 1}\right) = \sqrt {\frac {2}{n + 1}} \left[ \mathrm {D S T} (b) \right] _ {k},

其中 $[\mathrm{DST}(b)]_k$ 表示第 $k$ 个分量. 因此, 乘积 $y = Zb$ 可以通过 DST 来实现. 类似地, 乘积 $y = Z^{\mathsf{T}}b = Z^{-1}b$ 可以通过离散 Sine 反变换 IDST 实现, 即

y = Z ^ {\mathsf {T}} b = Z ^ {- 1} b = \left(\sqrt {\frac {2}{n + 1}}\right) ^ {- 1} \mathrm {I D S T} (b).

所以对于一维离散 Poisson 方程, 其解为

u = T _ {n} ^ {- 1} (h ^ {2} f) = (Z \Lambda^ {- 1} Z ^ {\intercal}) (h ^ {2} f) = h ^ {2} Z \Lambda^ {- 1} Z ^ {\intercal} f = h ^ {2} \mathrm {D S T} (\Lambda^ {- 1} \mathrm {I D S T} (b)).

可以通过公式 $T_{n} = Z\Lambda Z^{\top}$ 来计算一维离散 Poisson 矩阵的特征值: 记 $[\alpha_{1}, \alpha_{2}, \ldots, \alpha_{n}]^{\top}$ 和 $[\beta_{1}, \beta_{2}, \ldots, \beta_{n}]^{\top}$ 分别为 $Z^{\top}T_{n}$ 和 $Z^{\top}$ 的第一列, 则 $T_{n}$ 的特征值为

\lambda_ {i} = \frac {\alpha_ {i}}{\beta_ {i}}.

对应的MATLAB代码为:Lam=idst([2,-1,zeros(1,n-2)]')./idst_eye(n,1))

而对于二维离散 Poisson 方程, 我们需要计算 $(Z \otimes Z)b$ 和 $(Z^{\mathsf{T}} \otimes Z^{\mathsf{T}})b$ . 它们对应的是二维离散 Sine 变换和二维离散 Sine 反变换.

设 $b = [b_1^{\mathsf{T}}, b_2^{\mathsf{T}}, \ldots, b_n^{\mathsf{T}}]^{\mathsf{T}} \in \mathbb{R}^{n^2}$ , 其中 $b_k \in \mathbb{R}^{n \times n}$ . 令 $B = [b_1, b_2, \ldots, b_n] \in \mathbb{R}^{n \times n}$ , 则由 Kronecker 乘积的性质可知

(Z \otimes Z) b = (Z \otimes Z) \operatorname {v e c} (B) = \operatorname {v e c} (Z B Z ^ {\intercal}) = \operatorname {v e c} \left(\left(Z (Z B) ^ {\intercal}\right) ^ {\intercal}\right).

因此, 我们仍然可以使用 DST 来计算 $(Z \otimes Z)b$ . 类似地, 我们可以使用 IDST 来计算 $(Z^{\mathsf{T}} \otimes Z^{\mathsf{T}})b$ .

算法6.8.二维离散Poisson方程的快速方法

1: 计算 $b = h^{2}f$
2: $B = \mathrm{reshape}(b, n, n)$
3: $B_{1} = (Z^{\mathsf{T}}B)^{\mathsf{T}} = (\mathrm{IDST}(B))^{\mathsf{T}}$
4: $B_{2} = (Z^{\mathsf{T}}B_{1})^{\mathsf{T}} = (\mathrm{IDST}(B_{1}))^{\mathsf{T}}$
5: $b_{1} = (I \otimes \Lambda + \Lambda \otimes I)^{-1}\mathrm{vec}(B_{2})$
6: $B_{3} = \mathrm{reshape}(b_{1},n,n)$
7: $B_4 = (ZB_3)^\top = (\mathrm{DST}(B_3))^\top$
8: $B_{5} = (ZB_{4})^{\mathsf{T}} = (\mathrm{DST}(B_{4}))^{\mathsf{T}}$
9: $u = \text{reshape}(B_5, n^2, 1)$

MATLAB 程序见 Poisson_DST.m

6.3.5 求解方法小结

由于 Poisson 方程的特殊结构和性质, 除了常规求解方法以外, 人们还设计出了一些特殊的快速算法. 下表列出了求解二维 Poisson 方程的不同方法的比较 [30], 这里假定网格剖分为 $n \times n$ , 并记 $N = n^2$ .

6.3_应用_Poisson_方程求解