数值求解波动方程 [2]

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

周期边界条件

上一篇忘记说另外一个常用的边界条件了, 就是循环边界. 它对应着连续函数里的周期条件: $%5Cforall%20x%3A%5C%3Af(x%2BL)%3Df(x)$ , 其中 $L$ 是周期. 在使用有限差分法截断计算区间时, 边界外部的点重新映射到计算区间内, 即 $f_%7B0%2Ct%7D%3Df_%7Bl%2Ct%7D$ 和 $f_%7Bl%2B1%2Ct%7D%3Df_%7B1%2Ct%7D$ (julia 索引). 一般来说, 循环边界必定成对出现 (即波可以从两边穿过边界到达另一边), 但也可以实现波从一边穿过边界但另一边不能穿过. 下面是循环边界的实现例子:

数值求解的两个问题

在上图可以看到数值求解的一个问题. 上图使用的相速度 $c%3D1$ , 并且初值条件为 $u(x%2C0)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D1%2B%5Ccos(20%5Cpi%20x)%2C%5C%2C%7Cx%7C%3C0.05%5C%5C0%2C%5C%2C%7Cx%7C%5Cgeq0.05%5Cend%7Bmatrix%7D%5Cright.$ 和 $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D(x%2C0)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D20%5Cpi%5Csin(20%5Cpi%20x)%2C%5C%2C%7Cx%7C%3C0.05%5C%5C0%2C%5C%2C%7Cx%7C%5Cgeq0.05%5Cend%7Bmatrix%7D%5Cright.$ , 利用已知的一维波动方程的数学解 (见《数理方程》), 可以得出在无限空间里的特解为 $u(x%2Ct)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D1%2B%5Ccos(20%5Cpi(x-t))%2C%5C%2C%7Cx-t%7C%3C0.05%5C%5C0%2C%5C%2C%7Cx-t%7C%5Cgeq0.05%5Cend%7Bmatrix%7D%5Cright.$ , 即这是单个向右传播的波包.

观察上图, 除了一个向右传播的波包, 还有一个向左传播的比较毛糙的波包. 这是由离散化造成的, 由连续的初值条件离散化得到求解系统的初状态时就已经损失了一部分信息了. 通过增加网格精度 (即减少 $%5CDelta%20x$ ) 可以缓解这个情况, 但这意味着 $%5CDelta%20t$ 也要一并减少以确保不会出现空间数值不稳定, 也就是计算成本会快速上升. 下图是使用两倍网格精度求相同条件的解:

使用数值求解还需要注意另外一样东西: 时间数值不稳定和高频分量. 在部分情况下, 不可以选择参数使得数值稳定, 那么这时候宁愿时间数值不稳定也不要空间数值不稳定.

继续以上面循环边界的情况为例, 比较波包宽度为 0.1 和 0.04 两种情况, 可以看到波包宽度为 0.04 的情况弥散现象更严重. 所以需要选择合适的精度以处理情景中的高频分量.

但就算选择参数恰好是数值稳定的, 也不代表就能很好地处理高频分量. 如下图所示, 当选择 $%5CDelta%20x%3D0.003$ 时, 波包宽度为 0.04 的情况仍然会数值不稳定.

非均匀介质

在研究波的传播时肯定少不了对非均匀介质的研究. 在原波动方程里的相速度 $c$ 为常数. 当在非均匀介质里, 相速度与空间位置有关, 表示为 $c(x)$ , 那么相应的波动方程为 $%5Cfrac%7B%5Cpartial%5E2%20u%7D%7B%5Cpartial%20t%5E2%7D%3Dc(x)%5E2%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20x%5E2%7D$ . 形式上与原波动方程相同, 但因为相速度与位置有关, 所以离散化后得到的系数也与位置有关, 那么离散化后得到 $u_%7Bx%2Ct%2B1%7D%3D%5Cfrac%7Bc_x%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x%5E2%7D(u_%7Bx-1%2Ct%7D%2Bu_%7Bx%2B1%2Ct%7D)%2B2%5Cleft(1-%5Cfrac%7Bc_x%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x%5E2%7D%5Cright)u_%7Bx%2Ct%7D-u_%7Bx%2Ct-1%7D$ .

关于连续波动方程的反射/透射系数可以见专栏底部.

需要注意应该避免在任意位置上出现空间数值不稳定, 亦即 $%5Cforall%20x%3A%20%5Cfrac%7Bc(x)%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x%5E2%7D%5Cleq1$ . 最后, 这种模型实现不适用于连续变化的介质, 当连续变化介质的相速度会被离散化为多层间断的相速度. 波会在间断的相速度处出现反射波 (上图所示), 但准确的连续变化介质是没有反射波的. 下图是一个连续变化介质的例子, 可以看到反射波随着时间越来越大, 这是极其错误的.

非均匀介质的一小点数学

因为《数理方程》里没有提到非均匀介质, 所以这里稍微小提一下.

假设无限空间里有两种介质在 $x_0$ 处为界, 相速度由 $c(x)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Dc_1%2C%5C%2Cx%3Cx_0%5C%5Cc_2%2C%5C%2Cx%3Ex_0%5Cend%7Bmatrix%7D%5Cright.$ 给出. 其中 $c_1%2C%5C%2Cc_2$ 为常数. 假设有波在无穷远处从左往右传播, 那么有 $u%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Du_i%2Bu_r%2C%5C%2Cx%3Cx_0%5C%5Cu_t%2C%5C%2Cx%3Ex_0%5Cend%7Bmatrix%7D%5Cright.$ , 如下:

原波动方程的特解为 $e%5E%7Bi(kx-%5Comega%20t)%7D$ , 其中 $%5Comega$ 为频率, $k$ 为波数 (wave number), 并且有关系 $%5Comega%3Dc%7Ck%7C$ . 当 $k%3E0$ 时, 解表示为从左往右传播的波, 反之为从右往左, 那么 $u$ 表示为:

$u_i%3De%5E%7Bi(%5Comega%20x%2Fc_1-%5Comega%20t)%7D%3B%5C%3Bu_r%3DAe%5E%7Bi(-%5Comega%20x%2Fc_1-%5Comega%20t)%7D%3B%5C%3Bu_t%3DBe%5E%7Bi(%5Comega%20x%2Fc_2-%5Comega%20t)%7D$

其中 $A%2C%20B$ 为待定系数.

在边界上, 应有 $%5Clim_%7Bx%5Crightarrow%20x_0%5E-%7Du(x%2Ct)%3D%5Clim_%7Bx%5Crightarrow%20x_0%5E%2B%7Du(x%2Ct)$ 和 $%5Clim_%7Bx%5Crightarrow%20x_0%5E-%7D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D(x%2Ct)%3D%5Clim_%7Bx%5Crightarrow%20x_0%5E%2B%7D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D(x%2Ct)$ , 即

$e%5E%7Bi(%5Comega%20x_0%2Fc_1-%5Comega%20t)%7D%2BAe%5E%7Bi(-%5Comega%20x_0%2Fc_1-%5Comega%20t)%7D%3DBe%5E%7Bi(%5Comega%20x_0%2Fc_2-%5Comega%20t)%7D$

$%5Cfrac%7B%5Comega%7D%7Bc_1%7De%5E%7Bi(%5Comega%20x_0%2Fc_1-%5Comega%20t)%7D-A%5Cfrac%7B%5Comega%7D%7Bc_1%7De%5E%7Bi(-%5Comega%20x_0%2Fc_1-%5Comega%20t)%7D%3DB%5Cfrac%7B%5Comega%7D%7Bc_2%7De%5E%7Bi(%5Comega%20x_0%2Fc_2-%5Comega%20t)%7D$

可以解得 $A%3D%5Cfrac%7Bc_2-c_1%7D%7Bc_2%2Bc_1%7D%5Calpha$ 和 $B%3D%5Cfrac%7B2c_2%7D%7Bc_2%2Bc_1%7D%5Cbeta$ , 其中 $%5Calpha%2C%5Cbeta$ 为相位部分, 不太重要. 可以由上面的示例进行验证:

下面是整个程序的代码, 比之前的也没改很多东西

标签：

数值求解波动方程 [2]

数值求解波动方程 [2]的评论 (共条)

你可能也喜欢这些文章

最新发布的文章

数值求解波动方程 [2]

本文作者的其他文章

数值求解波动方程 [2]的评论 (共 条)

你可能也喜欢这些文章

最新发布的文章

数值求解波动方程 [2]的评论 (共条)