数值求解波动方程 [4]

2022-08-28 01:18 作者:nyasyamorina 0人读过 | 我要投稿

这篇里讲述了另外一种求得 PML 的方法以及简单实现二维波动方程的数值模拟.

另一个 PML 求解

在上一篇求一维 PML 时, PML 变换为 $%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D%5Cmapsto%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%5Ctilde%20x%7D%3D%5Cfrac%7B1%7D%7B1%2Bi%5Cfrac%7B%5Csigma%7D%7B%5Comega%7D%7D%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D$ , 以及在假设 $v%3De%5E%7Bi(kx-%5Comega%20t)%7D$ 时得到关系 $%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3D-i%5Comega%20v$ , 并且在最后恰好全部虚数单位 $i$ 都被抵消, 留下一个纯实数 (复数也可以) 的方程. 在部分情况下, 需要施加 PML 变换的波动方程比较复杂 (比如强迫振动或粘度等), 造成虚数单位无法恰好被消除, 留下一个需要在复数域上求解的方程. 也不是说数值求解复数方程不可以, 只不过考虑到内存和计算开销, 复数真的不太行.

实际上, 考虑傅里叶变换 $%5Cmathcal%20F%20f%3D%5Cmathcal%20F%5C%7Bf(t)%5C%7D(%5Comega)%3D%5Cint_%7B%5Cmathbb%20R%7Df(t)e%5E%7B-i%5Comega%20t%7Ddt%2C%5C%2C%5Comega%5Cin%5Cmathbb%20R$ , 无论 $v$ 如何都有 $%5Cmathcal%20F%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3D-i%5Comega%5Cmathcal%20Fv$ , 这与上面提到的关系是等价的. 另外再考虑拉普拉斯变换 $%5Cmathcal%20Lf%3D%5Cmathcal%20L%5C%7Bf(t)%5C%7D(s)%3D%5Cint_0%5E%7B%2B%5Cinfty%7Df(t)e%5E%7Bst%7Ddt%2C%5C%2Cs%5Cin%5Cmathbb%20C$ , 除了积分域和定义域不一样外, 它跟傅里叶变换是相同的. 在波动方程里, 不需要太关心 $t%3C0$ 部分 (这里只讨论正演而忽略反演), 那么把 $-i%5Comega%3Ds$ 替换掉所有出现的频率得到新的 PML 变换 $%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D%5Cmapsto%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%5Ctilde%20x%7D%3D%5Cfrac%7B1%7D%7B1%2B%5Cfrac%7B%5Csigma%7D%7Bs%7D%7D%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D$ 和关系 $%5Cmathcal%20L%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3Ds%5Cmathcal%20Lv$ . 下面继续以一维波动方程做例子:

对一维波动方程 $%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20x%5E2%7D$ 施加拉普拉斯变换得到 $%5Cmathcal%20L%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20t%5E2%7D%3D%5Cmathcal%20L%5Cleft(c%5E2%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20x%5E2%7D%5Cright)$ , 即 $s%5E2%5Cmathcal%20Lu%3Dc%5E2%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D%5Cmathcal%20Lu$ (偏微分算符和拉普拉斯变换的线性), 这个就是波动方程在 s 域的形式. 对 s 域的波动方程施加 PML 变换得 $s%5E2%5Cmathcal%20Lu%3D%5Cfrac%7Bc%5E2%7D%7B1%2B%5Cfrac%5Csigma%20s%7D%5Cfrac%5Cpartial%7B%5Cpartial%20x%7D%5Cleft(%5Cfrac1%7B1%2B%5Cfrac%5Csigma%20s%7D%5Cfrac%5Cpartial%7B%5Cpartial%20x%7D%5Cmathcal%20Lu%5Cright)$ , 两边同乘 $%5Cfrac%7Bs%2B%5Csigma%7D%7Bs%5E2%7D$ 得 $s%5Cmathcal%20L%20u%2B%5Csigma%5Cmathcal%20Lu%3Dc%5E2%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D%5Cleft(%5Cfrac1%7Bs%2B%5Csigma%7D%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D%5Cmathcal%20Lu%5Cright)$ . 引入新量 $v$ 使得 $%5Cmathcal%20Lv%3A%3D%5Cfrac1%7Bs%2B%5Csigma%7D%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D%5Cmathcal%20Lu$ , 即得到 $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Ds%5Cmathcal%20Lu%2B%5Csigma%5Cmathcal%20Lu%3Dc%5E2%5Cfrac%5Cpartial%7B%5Cpartial%20x%7D%5Cmathcal%20Lv%5C%5Cs%5Cmathcal%20Lv%2B%5Csigma%5Cmathcal%20Lv%3D%5Cfrac%5Cpartial%7B%5Cpartial%20x%7D%5Cmathcal%20Lu%5Cend%7Bmatrix%7D%5Cright.$ , 把方程组从 s 域转为时域得 $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D%2B%5Csigma%20u%3Dc%5E2%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20x%7D%5C%5C%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%2B%5Csigma%20v%3D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D%5Cend%7Bmatrix%7D%5Cright.$ , 可以看到这与之前得到的一模一样.

尽管两个方法得到的结果都是一样的, 但在 s 域上处理 PML 变换是非常通用的, 至少不需要假设一些条件.

二维波动方程

二维波动方程为 $%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5Cleft(%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20x_1%5E2%7D%2B%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20x_2%5E2%7D%5Cright)$ , 其中 $u(x_1%2Cx_2%2Ct)$ 为所求波的解. 习惯上使用 $x%2Cy$ 作为二维笛卡尔坐标, 但是 unicode 没有下标 y, 所以这里使用 $x_1%2Cx_2$ 作为坐标了 (方便以后的编程实现). 跟一维波动方程类似, 二维波动方程离散化为

$%5Cfrac%7Bu_%7Bx_1%2Cx_2%2Ct%2B1%7D%2Bu_%7Bx_1%2Cx_2%2Ct-1%7D-2u_%7Bx_1%2Cx_2%2Ct%7D%7D%7B%5CDelta%20t%5E2%7D%3Dc%5E2%5Cleft(%5Cfrac%7Bu_%7Bx_1%2B1%2Cx_2%2Ct%7D%2Bu_%7Bx_1-1%2Cx_2%2Ct%7D-2u_%7Bx_1%2Cx_2%2Ct%7D%7D%7B%5CDelta%20x_1%5E2%7D%2B%5Cfrac%7Bu_%7Bx_1%2Cx_2%2B1%2Ct%7D%2Bu_%7Bx_1%2Cx_2-1%2Ct%7D-2u_%7Bx_1%2Cx_2%2Ct%7D%7D%7B%5CDelta%20x_2%5E2%7D%5Cright)$

有一说一, 确实有点长了. 整理后得到

$u_%7Bx_1%2Cx_2%2Ct%2B1%7D%3Dc_1(u_%7Bx_1%2B1%2Cx_2%2Ct%7D%2Bu_%7Bx_1-1%2Cx_2%2Ct%7D)%2Bc_2(u_%7Bx_1%2Cx_2%2B1%2Ct%7D%2Bu_%7Bx_1%2Cx_2-1%2Ct%7D)%2Bc_3u_%7Bx_1%2Cx_2%2Ct%7D-u_%7Bx_1%2Cx_2%2Ct-1%7D$

其中 $c_1%3D%5Cfrac%7Bc%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x_1%5E2%7D%3B%5C%3Bc_2%3D%5Cfrac%7Bc%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x_2%5E2%7D%3B%5C%3Bc_3%3D2(1-c_1-c_2)$ . 在离散化后还需要截断计算区域, 截断方式与一维的一样, 但只不过现在不管是截断一个维度, 而是两个. 但是值得留意的是, 边界围成的形状不一定只能是正方形长方形, 但奇形怪状的边界很多时候需要额外计算边界处的离散化, 所以这里只展示正方形的边界了.

二维实现

在模拟二维波动方程里, 有一个比较大的问题就是计算速度, 因为网格大小会随着空间精度和大小的增大呈平方上升, 所以与其放在 cpu 上面跑, 不如放在 gpu 上面.

julia 的 cuda 支持在 https://cuda.juliagpu.org/stable 有详细记述. 一般来说有显卡驱动时直接在 REPL 里 `]add CUDA` 就可以安装了, 并且 julia 函数可以编译为 cuda 代码, 而不用像 python 那样在字符串里写 C++ CUDA 函数. 下面是 Wave2D 相关的实现, 在代码里上下边界使用了 Neumann 边界, 而左右边界则是循环边界.

不过在渲染 gif 的时候无论如何都不能把文件大小降下来, 明明隔壁 mp4 都超级小了, 怎么 gif 就那么大 (mp4: 211K, gif: 13.8M). 所以下面贴图片算了, 但是代码是直接渲染视频的.

观察到图中, 波包在空间扩散时, 在后面留下了不为零的"尾迹", 并且尾迹不会消失 (可以参考《数理方程》) , 亦即不符合 Huggens 原理, 这种现象称为波的弥散或有后效现象, 并且这是偶数维度空间特有的. 下图是沿轴方向的切片, 可以更直观地看到这个现象:

在上面的代码实现里, 已经包括了不均匀介质, 下图为空间中心有一个半径为 0.1 的实心玻璃圆.

渲染的代码稍微有点长, 所以这里单独拿出来贴一份:

数值稳定性

因为二维波动方程有三个维度 (两个空间一个时间), 那么应该讨论三个维度的数值稳定性. 据自己的观察经验, 当 $%5Cleft.%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D%5Cright%7C_%7Bt%3D0%7D%3D0$ 时, $c_3%3C0$ 会造成空间数值不稳定, 否则 $c_3%3C1$ 是空间数值不稳定 (? 我也不知道为何). 在空间维度上, $c_1%2Cc_2$ 是衡量时间不稳定的指标, 与一维波动方程类似, 这两个数值在 $%3C1$ 时会在相应维度上产生时间不稳定, 但必须确保全局上不能有空间不稳定, 即 $c_3%3D2(1-c_1-c_2)%5Cgeq0$ , 所以任何维度都不能不产生时间不稳定.