多重网格算法：原理、实现与测试

### 多重网格算法：原理、实现与测试 #### 1. 多重网格算法基础 在求解偏微分方程时，我们常基于附近点的值 $u(x_p, y_q)$ 来近似二阶偏导数。例如，对于 $\frac{\partial^2 u}{\partial x^2}(x_i, y_j)$ 有以下近似： $\frac{\partial^2 u}{\partial x^2}(x_i, y_j) \approx \frac{u(x_{i + 1}, y_j) - 2u(x_i, y_j) + u(x_{i - 1}, y_j)}{h^2}$ $\frac{\partial^2 u}{\partial y^2}(x_i, y_j) \approx \frac{u(x_i, y_{j + 1}) - 2u(x_i, y_j) + u(x_i, y_{j - 1})}{h^2}$ 这些近似的误差分别约为 $\frac{\partial^4 u}{\partial x^4} \cdot \frac{h^2}{12}$ 和 $\frac{\partial^4 u}{\partial y^4} \cdot \frac{h^2}{12}$。由此，偏微分方程的精确解应满足： $\frac{u(x_{i + 1}, y_j) + u(x_{i - 1}, y_j) + u(x_i, y_{j + 1}) + u(x_i, y_{j - 1}) - 4u(x_i, y_j)}{h^2} \approx f(x_i, y_j)$ 为了计算 $u(x_i, y_j)$ 的近似值 $u_{ij}$，我们将上述近似替换为等式，得到方程组： 对于区域 $R$ 内的每个 $(x_i, y_j)$，有 $f(x_i, y_j) = \frac{u_{i + 1, j} + u_{i - 1, j} + u_{i, j + 1} + u_{i - 1, j - 1} - 4u_{i, j}}{h^2}$ 对于区域 $R$ 边界上的每个 $(x_i, y_j)$，有 $g(x_i, y_j) = u_{i, j}$ $h$ 越小，近似值 $u(x_i, y_j) \approx u_{ij}$ 越精确，但线性方程组的规模也越大。例如，若 $R$ 是一个 $1 \times 1$ 的正方形，取 $h = 0.01$，线性系统将有 $99^2 \approx 10^4$ 个未知数。对于三维问题，未知数数量可能轻松达到数百万甚至数十亿。 #### 2. 多重网格方法概述 多重网格方法使用一系列“网格”，每个网格都有其对应的矩阵。这些网格从最细网格（与待求解的线性系统矩阵相关）到最粗网格（与一个小矩阵相关），编号为 $0, 1, \cdots, N - 1$。由于网格数量随未知数数量呈对数增长，$N$ 通常不超过 20。与网格 $k$ 相关的矩阵是 $A_k$，其中 $A_0 = A$ 是我们要求解的线性系统中的矩阵。 对于每个网格，我们需要一个平滑算子 $R_k$ 作为 $A_k$ 的预条件子。实际上，有两种平滑算子：预平滑器 $R_k$ 和后平滑器 $S_k$。矩阵 $A_k$ 是一个 $n_k \times n_k$ 的方阵，且 $n_0 \geq n_1 \geq \cdots \geq n_{N - 1}$，其中 $n_0 = n$ 是原始系统中的未知数数量。为了连接这些网格，有插值矩阵 $I_{k}^{k + 1}$（$n_k \times n_{k + 1}$）和限制算子 $I_{k + 1}^{k}$（$n_{k + 1} \times n_k$）。以下是一个抽象的多重网格算法： ```cpp multigrid(k,bk)/* Solve Akxk = bk */ { if k = N solve AN xN = bN directly else { xk ← 0 repeat ν1 times: xk ← Rk(xk, bk) rk+1 ← I k+1 k (bk − Akxk) ek ← I k k+1 multigrid(k + 1,rk+1) xk ← xk + ek repeat ν2 times: xk ← Sk(xk, bk) } return xk } ``` 操作 $x_k \leftarrow R_k(x_k, b_k)$ 和 $x_k \leftarrow S_k(x_k, b_k)$ 可以是高斯 - 赛德尔方法、高斯 - 雅可比方法、SOR 方法、它们的反向版本，或其他线性迭代方法，如不完全 LU 分解（ILU）的迭代细化。这些操作通常可以用矩阵 $\tilde{R}_k$ 和 $\tilde{S}_k$ 表示： $R_k(x_k, b_k) = x_k + \tilde{R}_k(b_k - A_kx_k)$ $S_k(x_k, b_k) = x_k + \tilde{S}_k(b_k - A_kx_k)$ 为了保证整体多重网格方法实现对称预条件子，我们需要 $\tilde{S}_k = \tilde{R}_k^T$，所有 $A_k$ 对称，且 $I_{k + 1}^{k} = (I_{k}^{k + 1})^T$。这是一个 V 循环多重网格算法，使用相同的组件还可以构造其他变体，如 W 循环。 #### 3. 框架实现 我们使用 C++ 语言和模板来实现这个框架。主模板类是一个抽象类，所有组件（函数 $A_k$、$I_{k + 1}^{k}$、$I_{k}^{k + 1}$、$R_k$、$S_k$ 和支持函数）都可以在派生类中重载。这些“组件”由计算矩阵 - 向量积的函数表示，例如计算 $y = A_kx$。 主模板参数是向量和标量类型，这使我们可以使用相同的代码处理单精度、双精度和扩展精度（如果可用），以及不同的向量表示方式。我们选择使用印第安纳大学开放系统实验室开发的矩阵模板库（MTL）进行最终实现，其他可选库包括 MV++、SparseLib++、Blitz++、uBLAS 等。 在向量类上需要执行一些操作，有些操作可以由 STL 向量模板类支持，有些则不能。在设计算术运算时，我们选择使用输出参数的“三参数”函数，以提高效率。 我们还使用了类型绑定器技术。类型绑定器利用 C++ 类中可以定义类型的特性，使这些类型“局部”于类。以下是多重网格框架的类型绑定器的主要部分： ```cpp template< typename T, class VectorType, class MatrixType, class PermType > class TypeBinder { public: typedef TypeBinder< T, VectorType, MatrixType, PermType > Types; typedef T value type; // Scalar Type typedef unsigned int int type; // Integer Type typedef MatrixType Matrix; // Matrix (dense) typedef VectorType Vector; // Vector (dense) typedef PermType Perm; // Permutation or pivot typedef typename Vector::size type size type; // Grid class contains vectors for internal use ... }; //TypeBinder ``` 这个类型绑定器包含了创建多重网格类所需的所有类型。由于值类型 $T$ 是模板参数，我们可以轻松切换单精度、双精度或扩展精度，所有计算都会传播相应的精度级别。 类型绑定器还包含一个 Grid 类，用于存储每个级别的残差、右侧项和近似解： ```cpp struct Grid { Vector x, b, r; // Constructor of dimension dim Grid(size type dim) { Vector x(dim), b(dim), r(dim); } // Empty constructor Grid() { } // resize -- sets new dimension void resize(size type new dim) { x.resize(new dim); b.resize(new dim); r.resize(new dim); } }; // class Grid ``` 主多重网格模板类既是模板类又是抽象类。多重网格方法的主要“组件”是纯虚函数，必须在派生类中定义。定义多重网格 V 循环的例程不是虚函数，因此派生类只需定义多重网格方法的组件，基类会处理如何组合它们。 以下是多重网格模板类的大纲： ```cpp template < type ```