数值求解波动方程 [1]

2022-08-16 16:13 作者:nyasyamorina 0人读过 | 我要投稿

关于这个系列的专栏, 与其说是教程, 不如说更像是笔记, 因为数值求解的内容实在在太多太多了, 我也只是学到什么写什么而已, 希望可以帮到有需要的人?

波动方程

波动方程是一个用于描述像弦振动, 薄膜振动, 声波传播和地震波等现象的方程, 甚至麦克斯韦方程组也可以退化为波动方程, 从而使用波动方程去研究电磁行为.

一般地, 波动方程的形式为 $%5Cfrac%7B%5Cpartial%5E2%20u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5CDelta%20u$ , 其中 $%5CDelta%3D%5Cfrac%7B%5Cpartial%5E2%7D%7B%5Cpartial%20x_1%5E2%7D%2B%5Cfrac%7B%5Cpartial%5E2%7D%7B%5Cpartial%20x_2%5E2%7D%2B%5Ccdots$ 为拉普拉斯算子, $c$ 为波在介质里传播的速度, 又叫相速度. 那么求解函数 $u(x_1%2Cx_2%2C%5Ccdots%2Ct)$ 即为所求波的解. 关于波动方程的推导, 数学解等, 在《数理方程》里就有, 而且也不是这个专栏的重点. 因为波动方程是二阶微分方程, 所以在给出初值条件 (一般来说是 $t%3D0$ 的 $u$ 和 $%5Cpartial%20u%2F%5Cpartial%20t$ ) 时可以求出相应的特解. 另外波动方程具有很多种变形: 边值条件, 强迫振动, 非均匀介质等.

数值求解波动方程就是在已知初值条件时使用计算机求得相应的特解, 主要实现方法有有限元法和有限差分法等. 相较于有限元法, 有限差分法更容易, 所以这系列的专栏是使用有限差分法进行波动方程的数值求解.

有限差分法

由于计算机无法处理无限的数据, 所以当进行数值模拟时必须对模型进行离散化 (避免无限小的间隔) 和截断 (避免无限大的区间).

由导数的定义 $%5Cfrac%7Bdf%7D%7Bdx%7D(x)%3D%5Clim_%7Bdx%5Crightarrow0%7D%5Cfrac%7Bf(x%2Bdx)-f(x)%7D%7Bdx%7D$ 可以得出有限差分法的离散化为 $%5Cfrac%7Bf(x%2B%5CDelta%20x)-f(x)%7D%7B%5CDelta%20x%7D%5Capprox%5Cfrac%7Bdf%7D%7Bdx%7D(x)%2C%5C%2C%5CDelta%20x%3E0$ , 但 $%5Clim_%7B%5CDelta%20x%5Crightarrow0%7D$ 暗含着条件 $%5Clim_%7Bdx%5Crightarrow0%5E-%7D%3D%5Clim_%7Bdx%5Crightarrow0%5E%2B%7D$ , 所以有 $%5Cfrac%7Bf(x)-f(x-%5CDelta%20x)%7D%7B%5CDelta%20x%7D%5Capprox%5Cfrac%7Bf(x%2B%5CDelta%20x)-f(x)%7D%7B%5CDelta%20x%7D$ , 于是更准确的离散化应为 $%5Cfrac%7Bdf%7D%7Bdx%7D(x)%5Capprox%5Cfrac%7B1%7D%7B2%7D%5Cleft(%5Cfrac%7Bf(x%2B%5CDelta%20x)-f(x)%7D%7B%5CDelta%20x%7D%2B%5Cfrac%7Bf(x)-f(x-%5CDelta%20x)%7D%7B%5CDelta%20x%7D%5Cright)%3D%5Cfrac%7Bf(x%2B%5CDelta%20x)-f(x-%5CDelta%20x)%7D%7B2%5CDelta%20x%7D$ . 类似地可以得出二阶导数的离散化: $%5Cfrac%7Bd%5E2f%7D%7Bdx%5E2%7D(x)%5Capprox%5Cfrac%7Bf(x-%5CDelta%20x)-2f(x)%2Bf(x%2B%5CDelta%20x)%7D%7B%5CDelta%20x%5E2%7D$ . 虽然可以继续推到更高阶导数的离散化, 但是波动方程只需二阶就足够了.

函数在离散化后可以使用数组储存, 假设 $f(x)$ 在数组表示为 $f_x$ , 其中下标为索引, 那么有 $f(x%2B%5CDelta%20x)$ 在数组中为 $f_%7Bx%2B1%7D$ . 所以导数的离散化可以写为 $%7B%5Cfrac%7Bdf%7D%7Bdx%7D%7D_x%3D%5Cfrac%7Bf_%7Bx%2B1%7D-f_%7Bx-1%7D%7D%7B2%5CDelta%20x%7D$ , 二阶导数为 $%7B%5Cfrac%7Bd%5E2f%7D%7Bdx%5E2%7D%7D_x%3D%5Cfrac%7Bf_%7Bx-1%7D%2Bf_%7Bx%2B1%7D-2f_x%7D%7B%5CDelta%20x%5E2%7D$ .

有限差分法的区间截断等价于微分问题里的边值条件, 通过控制边界的行为去避免读取区间外的数据, 从而使计算范围限制在区间内. 常用的边值条件有 Dirichlet 条件和 Neumann 条件, 假设 $x_1$ 为边界, 那么 Dirichlet 条件表示为 $f(x_1)%3Dk$ , Neumann 条件表示为 $%5Cfrac%7Bdf%7D%7Bdx%7D(x_1)%3Dk$ , 其中 $k$ 为常数.

在使用数值求解时需要注意一个严重的问题: 数值不稳定. 它会造成求解误差较大或者直接发散出现 inf, nan. 由于我自己并没有太深究这个问题, 所以并没有具体分析过数值不稳定, 但从个人经验来说, 数值不稳定与参数选择和模型的具体实现, 甚至初值有关.

实现

下面以一维波动方程为例, 一维波动方程为 $%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20x%5E2%7D$ , 其中函数 $u(x%2Ct)$ 即为所求的解. 方程离散化得 $%5Cfrac%7B1%7D%7B%5CDelta%20t%5E2%7D(u_%7Bx%2Ct-1%7D%2Bu_%7Bx%2Ct%2B1%7D-2u_%7Bx%2Ct%7D)%3D%5Cfrac%7Bc%5E2%7D%7B%5CDelta%20x%5E2%7D(u_%7Bx-1%2Ct%7D%2Bu_%7Bx%2B1%2Ct%7D-2u_%7Bx%2Ct%7D)$ , 移项得 $u_%7Bx%2Ct%2B1%7D%3D%5Cfrac%7Bc%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x%5E2%7D(u_%7Bx-1%2Ct%7D%2Bu_%7Bx%2B1%2Ct%7D)-u_%7Bx%2Ct-1%7D%2B2%5Cleft(1-%5Cfrac%7Bc%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x%5E2%7D%5Cright)u_%7Bx%2Ct%7D$ . 其中两块系数可以提前计算以提升效率.

然后需要设置边值条件, 这里在左边界使用 Dirichlet 条件 $u(%5Ccdot%2Ct)%3D0$ , 右边界使用 Neumann 条件 $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D(%5Ccdot%2Ct)%3D0$ , 假设离散化后的数组索引范围为 $%5B1%2Cl%5D$ (大部分编程语言是 $%5B0%2Cl-1%5D$ , 我都老 julia 人了), 那么根据边值条件有 $u_%7B0%2Ct%7D%3D0$ 和 $u_%7Bl%2B1%2Ct%7D%3Du_%7Bl-1%2Ct%7D$ (对于其他语言是 $u_%7B-1%2Ct%7D%3D0$ 和 $u_%7Bl%2Ct%7D%3Du_%7Bl-2%2Ct%7D$ ).

最后因为波动方程是二阶微分方程, 所以至少需要提供两个初值条件, 这里以 $u(%5Ccdot%2C0)$ 和 $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D(%5Ccdot%2C0)$ 为例, 那么有 $u_%7B%5Ccdot%2C0%7D%3Du(%5Ccdot%2C0)$ 和 $u_%7B%5Ccdot%2C-1%7D%3Du_%7B%5Ccdot%2C0%7D-c%5CDelta%20t%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D(%5Ccdot%2C0)$ .

在这个特定的实现里, 当满足 $k%3D%5Cfrac%7Bc%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x%5E2%7D%3E1$ 时, 会造成空间数值不稳定, 值会以指数速度发散; 当 $k%3C1$ 时为时间数值不稳定, 时间不稳定不会造成值发散, 但会产生类似色散的现象; 只有恰好 $k%3D1$ 时才是数值稳定的. 如下图所示