超级会员免费看
数值传播与快速线性求解器
1. 数值传播中的数值扩散与色散
在数值传播的研究中，我们常常会遇到数值扩散和色散的问题。以Crank - Nicolson/中心差分法为例，它在不同时间步下的表现如图所示。以下是一段相关的代码：
for(int i=1; i<N-1; i++)
q[i-1] = uold[i] - 0.25*CFL*(uold[i+1]-uold[i-1])+
0.5*DN*(uold[i+1] - 2.0*uold[i] + uold[i-1]);
ThomasAlgorithm-per(N-1,c2,1.0e0+DN,c1,&unew[1],q);
unew[N-1] = unew[N-2];
delete[] q;
需要注意的是代码中“&”的使用，它表示取地址。在上述代码中，我们将
unew
数组第二个元素的地址传递给
ThomasAlgorithm - per
函数。当然，也可以使用指针算术来实现相同的功能，如下所示：
ThomasAlgorithm-per(N-1,c2,1.0e0+DN,c1,unew+1,q);
2. MPI非阻塞通信
2.1 阻塞与非阻塞通信
在之前的有限差分并行实现中，我们使用过
MPI Send
/
MPI Recv
或
MPI Sendrecv replace
命令来交换信息，这些都是阻塞通信，即它们在确保发送/接收缓冲区内容可以安全修改或使用之前不会返回。而MPI还提供了非阻塞版本的函数
MPI Isend
和
MPI Irecv
，“I”代表立即。这些函数允许进程先发布发送或接收信息的请求，之后再调用
MPI Wait
函数来完成发送/接收操作。
2.2 MPI Isend/MPI Irecv/MPI Wait的工作原理
假设进程0需要向进程1发送信息，但由于两个进程运行的特定算法，可能会出现同步不匹配的情况。此时，进程0发起
MPI Isend
请求（发布发送消息的意愿），然后继续执行不需要发送缓冲区内容的任务。当进程0无法继续执行而必须保证发送缓冲区内容可以修改时，调用
MPI Wait
等待事务完成。进程1类似，先通过
MPI Irecv
发布接收消息的意愿，当需要接收缓冲区内容才能继续执行时，调用
MPI Wait
等待事务完成。完成
MPI Wait
后，发送者可以修改发送缓冲区，接收者可以使用接收缓冲区中的数据。其通信过程如图所示：
graph LR
classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
A(Process 0):::process -->|MPI_Isend| B(Post Send):::process
C(Process 1):::process -->|MPI_Irecv| D(Post Recv):::process
B -->|MPI_Wait| E(Transaction Complete):::process
D -->|MPI_Wait| E
2.3 函数调用语法及参数解释
- MPI Isend
int MPI-Isend(
void* message /* in */,
int count /* in */,
MPI-Datatype datatype /* in */,
int dest /* in */,
int tag /* in */,
MPI-Comm comm /* in */,
MPI-Request* request /* out */)
- MPI Irecv
int MPI-Irecv(
void* message /* out */,
int count /* in */,
MPI-Datatype datatype /* in */,
int source /* in */,
int tag /* in */,
MPI-Comm comm /* in */,
MPI-Request* request /* out */)
- MPI Wait
int MPI-Wait(
MPI-Request* request /* in/out */,
MPI-Status* status /* out */)
参数解释如下表：
| 参数 | 含义 |
| ---- | ---- |
| message | 发送/接收缓冲区的起始地址 |
| count | 发送/接收缓冲区中的元素数量 |
| datatype | 发送缓冲区中元素的数据类型 |
| source | 发送数据的进程编号 |
| dest | 接收数据的进程编号 |
| tag | 消息标签 |
| comm | 通信器 |
| request | 通信请求 |
| status | 状态对象 |
2.4 使用示例
int mynode, totalnodes;
int datasize; // 要发送/接收的数据单元数量
int sender; // 发送进程的编号
int receiver; // 接收进程的编号
int tag; // 整数消息标签
MPI-Status status; // 用于保存状态信息的变量
MPI-Request request; // 用于维护isend/irecv信息的变量
MPI-Init(&argc,&argv);
MPI-Comm-size(MPI-COMM-WORLD, &totalnodes);
MPI-Comm-rank(MPI-COMM-WORLD, &mynode);
// 确定数据大小
double * databuffer = new double[datasize];
// 在发送/接收进程中填充sender, receiver, tag，并在发送进程中填充databuffer
if(mynode==sender)
MPI-Isend(databuffer,datasize,MPI-DOUBLE,receiver,tag,MPI-COMM-WORLD,&request);
if(mynode==receiver)
MPI-Irecv(databuffer,datasize,MPI-DOUBLE,sender,tag,MPI-COMM-WORLD,&request);
// 发送者/接收者可以执行不涉及databuffer数组的各种操作
MPI-Wait(&request,&status); // 同步以验证数据是否已发送
// 发送/接收完成
2.5 注意事项
-
发送者和接收者的消息数组通常应具有相同的类型，且大小至少为
datasize。 - 大多数情况下，发送类型和接收类型相同。
-
在
MPI Isend调用之后、MPI Wait调用之前，不应更改message的内容。 -
在
MPI Irecv调用之后、MPI Wait调用之前，不应使用message的内容。 -
MPI Send可以被MPI Irecv/MPI Wait接收，MPI Recv可以从MPI Isend/MPI Wait获取信息。 - 标签可以是0到32767之间的任意整数。
-
MPI Irecv可以使用通配符MPI ANY TAG来接收任意标签的消息，但MPI Isend必须指定具体的标签。 -
MPI Irecv可以使用通配符MPI ANY SOURCE来接收任意源的消息，但MPI Isend必须指定目标进程的编号。
3. 作业问题
以下是一系列相关的作业问题，涵盖了数值方法的稳定性分析、不同差分方案的实现以及非线性问题的求解等方面：
1. 使用冯·诺伊曼分析证明，蛙跳时间差分与中心差分结合求解$\frac{\partial \phi}{\partial t} + \frac{\partial \phi}{\partial x} = 0$（具有周期边界条件）会得到一个条件稳定的方案。能否从相应的等效微分方程得出相同的结论？
2. 使用冯·诺伊曼分析分析蛙跳时间差分与迎风差分结合求解$\frac{\partial \phi}{\partial t} + \frac{\partial \phi}{\partial x} = 0$（具有周期边界条件）的稳定性，并与从相应微分方程得出的启发式结论进行比较。
3. 使用三阶Adams - Bashforth/中心差分方案，参数类似于某组参数，但分别取$N = 30$和$N = 100$个点，求解一维线性对流方程，并将结果与二阶Adams - Bashforth/中心差分方案的结果进行比较。
4. 使用蛙跳/中心差分方案，$C = 0.2$，$N = 30$，求解一维对流方程，分别考虑周期和非周期边界条件，判断两种情况下的解是否稳定。
5. 数值求解一维对流 - 扩散方程$\phi_t + \phi_x = \frac{1}{100}\phi_{xx}$，初始条件为$\phi(x, t = 0) = e^{-2x^2}$，边界条件为$\phi(-3, t) = 0$，$\frac{\partial \phi}{\partial x}(L, t) = 0$（$L = 3$和$L = 10$）。使用中心差分进行空间导数计算，显式欧拉法进行时间积分。可以选择网格间距$\Delta x = 0.01$或$\Delta x = 0.1$，比较$t = 2$和$t = 5$时的解。
6. 非线性Lax - Friedrichs方案和Tadmor修正：求解无粘Burgers方程$\frac{\partial \phi}{\partial t} + \frac{\partial F(\phi)}{\partial x} = 0$，其中$F(\phi) = \frac{1}{2}\phi^2(x, t)$，初始条件为$\phi(x, 0) = \sin(\pi x)$，$x \in [-1, 1]$。该问题的解在$t_s \approx 0.31$时会出现激波间断，比较在激波形成前后（$t = 0.2$和$t = 0.4$）使用Lax - Friedrichs方案和Tadmor修正得到的数值解，具体计算不同网格点数（$N = 50$，$100$，$200$，$400$）下的$L_1$范数和$L_2$范数误差。
7. 将
EF FirstOrderUpwind
函数重写为使用
MPI Isend
和
MPI Irecv
的并行函数，验证代码并修改函数以求解$u_t + au_x = 0$（$a$为非零常数，可正可负），分析当波传播方向改变时MPI代码需要做哪些修改。
8. 修改
CrankNicolson CentralDifference
函数，使其使用欧拉向后法代替Crank - Nicolson法，重复某示例并与Crank - Nicolson法的结果进行比较。
9. 创建一个新的Crank - Nicolson C++函数，使用四阶紧致方案进行空间微分，重复某示例，分析空间精度提高带来的差异。
10. 高斯锥问题：考虑具有空间变化对流速度$u = +y$，$v = -x$的对流 - 扩散方程，初始条件为$u(x, y, 0) = e^{-\frac{(x - x_0)^2 + (y - y_0)^2}{2\lambda^2}}$。在由$100×100$个网格点组成的网格上使用二阶中心差分进行时间积分，时间为$t = 2\pi$。
- 采用欧拉方法和二阶Adams - Bashforth/Crank - Nicolson时间步长，取Courant数$C = 0.01$和扩散数$D = 0.5$，获得数值解。
- 使用基于显式中点规则的半拉格朗日方法进行向后积分，重复上述步骤。
- 将时间步长分别增大10倍和20倍，重复上述步骤，并与之前的结果进行比较，观察现象。
11. 重复上一个问题，但使用新的初始条件$u(x, y, t = 0) = \begin{cases}-16(r_0^2 - \frac{1}{16}) & \text{if } r_0 < \frac{1}{4} \ 0 & \text{elsewhere}\end{cases}$，其中$r_0^2 = (x - x_0)^2 + (y - y_0)^2$，$(x_0, y_0)$为锥的初始中心位置。比较欧拉方法和半拉格朗日方法在等高线图和$L_2$误差方面的结果。
4. 快速线性求解器
4.1 高斯消元法
高斯消元法是求解线性系统$Ax = b$的一种有效方法，托马斯算法是其在三对角系统中的特殊情况。高斯消元法基于线性系统的“叠加原理”，即可以用矩阵$A$的行和对应$b$值的线性组合形成的等效方程替换原系统的方程。具体做法是，从每一行中减去其倍数，使$A$对角线以下的元素变为零。
4.1.1 示例
考虑一个$3×3$的线性系统：
$\begin{cases}x_1 + \frac{1}{2}x_2 + \frac{1}{3}x_3 = 3 \ \frac{1}{2}x_1 + \frac{1}{3}x_2 + \frac{1}{4}x_3 = 2 \ \frac{1}{3}x_1 + \frac{1}{4}x_2 + \frac{1}{5}x_3 = 1\end{cases}$
求解该系统分两个阶段：
-
阶段1
：选择主元$a_{11} = 1$，计算乘数$\ell_{21}^{(1)} = \frac{a_{21}}{a_{11}} = \frac{1}{2}$，$\ell_{31}^{(1)} = \frac{a_{31}}{a_{11}} = \frac{1}{3}$。将第一个方程乘以$\ell_{21}^{(1)}$后从第二个方程中减去，乘以$\ell_{31}^{(1)}$后从第三个方程中减去，得到新的方程组：
$\begin{cases}x_1 + \frac{1}{2}x_2 + \frac{1}{3}x_3 = 3 \ \frac{1}{12}x_2 + \frac{1}{12}x_3 = \frac{1}{2} \ \frac{1}{12}x_2 + \frac{4}{45}x_3 = 0\end{cases}$
-
阶段2
：选择新的主元$a_{22}^{(1)} = \frac{1}{12}$，乘数$\ell_{32}^{(2)} = \frac{\frac{1}{12}}{\frac{1}{12}} = 1$。将第二个方程乘以$\ell_{32}^{(2)}$后从第三个方程中减去，得到：
$\begin{cases}x_1 + \frac{1}{2}x_2 + \frac{1}{3}x_3 = 3 \ \frac{1}{12}x_2 + \frac{1}{12}x_3 = \frac{1}{2} \ \frac{1}{180}x_3 = -\frac{1}{2}\end{cases}$
通过回代求解，从第三个方程开始，依次得到$x_3 = -90$，$x_2 = 96$，$x_1 = -15$。
在矩阵形式下，该系统可表示为$Ax = LUx = b$，其中$L$是下三角矩阵，$U$是上三角矩阵：
$L = \begin{bmatrix}1 & 0 & 0 \ \frac{1}{2} & 1 & 0 \ \frac{1}{3} & 1 & 1\end{bmatrix}$，$U = \begin{bmatrix}1 & \frac{1}{2} & \frac{1}{3} \ 0 & \frac{1}{12} & \frac{1}{12} \ 0 & 0 & \frac{1}{180}\end{bmatrix}$
4.2 LU分解
对于一般的$n×n$线性系统$Ax = b$，需要$(n - 1)$个消元阶段得到上三角系统。假设每一步的主元$a_{(ii)}^{(k)} \neq 0$，后续会讨论涉及行和/或列主元选择的算法来放松这个条件。
4.2.1 消元过程
-
第一阶段消元后得到：
$\begin{cases}a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 \ a_{22}^{(1)}x_2 + \cdots + a_{2n}^{(1)}x_n = b_2^{(1)} \ \cdots \ a_{n2}^{(1)}x_2 + \cdots + a_{nn}^{(1)}x_n = b_n^{(1)}\end{cases}$
其中中间系数$a_{ij}^{(1)} = a_{ij} - a_{1j}\ell_{i1}^{(1)}$，$\ell_{i1}^{(1)} = \frac{a_{i1}}{a_{11}}$，右侧项$b_i^{(1)} = b_i - b_1\frac{a_{i1}}{a_{11}}$。 - 类似地，第二阶段消元得到$a_{ij}^{(2)}$和$b_i^{(2)}$，以此类推，直到第$(n - 1)$阶段得到$a_{nn}^{(n - 1)}$和$b_n^{(n - 1)}$。
4.2.2 求解步骤
整个求解过程包括三个主要步骤：
1. LU分解：$A = L·U$
2. 前向求解$y$：$Ly = b$
3. 回代求解$x$：$Ux = y$
4.2.3 伪代码实现
for k = 1, n - 1
for i = k + 1, n
ℓik = aik / akk (assuming akk ≠ 0)
for j = k + 1, n
aij = aij - ℓik * akj
endfor
bi = bi - ℓik * bk
endfor
endfor
计算复杂度方面，LU分解的操作次数约为$\frac{2}{3}n^3$（考虑加法和乘法），若只考虑乘法约为$\frac{n^3}{3}$。回代求解的操作次数为$O(n^2)$。
4.3 相关注意事项
- 内存使用 ：在上述伪代码中，为了减少内存使用，对矩阵$A$和右侧向量$b$进行了覆盖操作。但如果后续还需要使用矩阵$A$，则应避免这种做法。
- 带宽矩阵 ：对于带宽为$m$的带状矩阵（$m \ll n$），LU分解的操作次数为$O(m^2n)$，回代求解的操作次数为$O(mn)$。
-
与其他方法比较
：与Cramer行列式法（复杂度为$O(n!)$）相比，高斯消元法（复杂度为$O(\frac{2n^3}{3})$）在计算效率上有巨大优势。例如，在使用1 Gflops的处理器时，不同$n$值下两种方法的计算时间对比如下：
| $n$ | Cramer法时间 | 高斯消元法时间 |
| ---- | ---- | ---- |
| 3 | $\approx 6$纳秒 | $\approx 18$纳秒 |
| 10 | $\approx 3$毫秒 | $\approx 0.6$微秒 |
| 20 | $\approx 80$年 | $\approx 5$微秒 | - 行列式计算 ：高斯消元法还可以高效地计算矩阵$A$的行列式，因为$\det(A) = \det(L)\det(U) = 1·[u_{11}·u_{22} \cdots u_{nn}]$，其中$u_{ii}$是上三角矩阵$U$的对角元素。
4.4 特殊矩阵情况
在实际应用中，会遇到各种特殊矩阵，它们的性质会影响高斯消元法和LU分解的计算过程和效率。
4.4.1 希尔伯特矩阵
希尔伯特矩阵是一种特殊矩阵，其元素$a_{ij} = \frac{1}{i + j - 1}$，第$n$行是向量$[\frac{1}{n}, \frac{1}{n + 1}, \frac{1}{n + 2}, \cdots, \frac{1}{2n - 1}]^T$。当$n$较大时，最后一行的元素比第一行的元素小约三个数量级，这种巨大的差异导致矩阵病态，条件数很大（例如，当$n \geq 5$时，条件数大于$10^5$）。在使用高斯消元法求解希尔伯特矩阵对应的线性系统时，会因为病态问题导致计算结果的误差较大，需要特别注意。
4.4.2 带状矩阵
对于带宽为$m$的带状矩阵（$m \ll n$），前面已经提到其LU分解的操作次数为$O(m^2n)$，回代求解的操作次数为$O(mn)$。这是因为带状矩阵的非零元素集中在主对角线附近的带状区域内，在消元过程中可以避免对大量零元素的计算，从而显著减少计算量。例如，三对角矩阵就是一种特殊的带状矩阵，托马斯算法就是针对三对角矩阵的高效求解方法，它是高斯消元法在三对角系统中的特殊情况。
4.5 与其他分解方法的比较
除了LU分解，还有其他矩阵分解方法，如Gram - Schmidt QR分解和Householder方法实现的QR分解，下面对它们进行比较。
4.5.1 Gram - Schmidt QR分解
Gram - Schmidt QR分解将矩阵$A$分解为一个正交矩阵$Q$和一个上三角矩阵$R$，即$A = QR$。该算法的计算复杂度为$O(2n^3)$（包括加法和乘法），是LU算法的三倍。这是因为Gram - Schmidt过程需要进行多次向量的正交化操作，涉及大量的内积和向量运算，导致计算量较大。
4.5.2 Householder方法实现的QR分解
Householder方法也可以实现矩阵的QR分解，其计算复杂度为$O(\frac{4}{3}n^3)$，是LU算法的约两倍。Householder方法通过构造一系列的Householder变换矩阵，逐步将矩阵$A$变换为上三角矩阵$R$，同时得到正交矩阵$Q$。与Gram - Schmidt方法相比，Householder方法在数值稳定性上更有优势，因为它减少了舍入误差的积累。
4.6 并行计算中的应用
在并行计算中，高斯消元法和MPI通信函数结合可以实现线性系统的并行求解。例如，使用MPI的广播命令
MPI Bcast
可以在高斯消元过程中高效地传播主元信息，同时通过
MPI Send
、
MPI Recv
、
MPI Allgather
和
MPI Allreduce
等函数进行数据的交换和同步。
4.6.1 MPI Bcast的应用
MPI Bcast
用于在一个通信器内将一个进程的数据广播到所有其他进程。在高斯消元的每一步中，主进程可以使用
MPI Bcast
将当前的主元信息广播给所有从进程，使得所有进程都能同步进行消元操作。
4.6.2 数据交换和同步
在消元过程中，不同进程可能需要交换部分数据以完成计算。例如，在计算$a_{ij} = a_{ij} - \ell_{ik} * a_{kj}$时，进程可能需要获取其他进程中的$a_{kj}$值。这时可以使用
MPI Send
和
MPI Recv
进行点对点的数据交换，或者使用
MPI Allgather
将所有进程的数据收集到每个进程中。
MPI Allreduce
则可以用于对所有进程的数据进行归约操作，如求和、求最大值等，以完成全局的计算。
4.7 总结
高斯消元法和LU分解是求解线性系统的重要方法，它们基于线性系统的叠加原理，通过逐步消元将原系统转化为上三角系统，然后进行回代求解。在实际应用中，需要考虑矩阵的性质（如病态性、带宽等）对计算结果和效率的影响，同时可以结合并行计算技术进一步提高求解效率。与其他矩阵分解方法相比，高斯消元法和LU分解在计算复杂度和数值稳定性上有各自的优势。在使用这些方法时，还需要注意内存的使用和数据的同步问题，以确保计算的正确性和高效性。
下面是一个简单的流程图，展示了高斯消元法的主要步骤：
graph TD
classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
A(输入线性系统 Ax = b):::process --> B(LU分解: A = L·U):::process
B --> C(前向求解 y: Ly = b):::process
C --> D(回代求解 x: Ux = y):::process
D --> E(输出解 x):::process
通过以上的介绍，我们对数值传播中的数值扩散与色散、MPI非阻塞通信、作业问题以及快速线性求解器（特别是高斯消元法和LU分解）有了更深入的了解。这些知识在科学计算、工程模拟等领域有着广泛的应用，希望读者能够通过实践进一步掌握这些方法和技术。
扫一扫
2381
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
