带固定点源阻尼波动方程的有限差分法

Finite Difference Method for Damped Wave Equation with Fixed Point Source

摘要: 本文在周期边界条件下针对一类带有固定点源的阻尼波动方程提出了一种有限差分方法：在空间上采用中心差分格式进行近似，并在点源附近构造辅助函数进行修正，有效处理了解在该点的奇异性；在时间上则采用向后差分格式。数值实验结果验证了格式的有效性和收敛性。

Abstract: In this paper, a finite difference method is proposed for a class of damped wave equations with a fixed point source under periodic boundary conditions. The central difference scheme is used to approximate the space, and an auxiliary function is constructed near the point source to correct the singularity of the solution at this point. In time, the backward difference scheme is used. Numerical experimental results verify the effectiveness and convergence of the scheme.

1. 前言

偏微分方程是描述自然现象的核心数学工具，在数学物理研究中占据着重要地位。在许多实际场景中，系统会受到集中式的能量输入，例如激光点热源等。这类局部化的瞬时激励在数学上通常通过带有狄拉克$\delta$ 函数形式的点源项来建模。本文研究下述带有固定点源的阻尼波动方程。

$a{u}_{tt}+b{u}_t-c^2{u}_{xx}=\delta \left(x-\overline{x}\right)F\left(t\right),x_a<x<x_b,0<t\le T,$ (1)

满足初始条件

$u\left(x,0\right)=u_0\left(x\right),$ (2)

$u_t\left(x,0\right)=v_0\left(x\right),$ (3)

以及周期边界条件，其中$a$ 、$b$ 、$c$ 为给定的常数，$\delta \left(x-\overline{x}\right)$ 表示位于$\overline{x}$ 处的狄拉克$\delta$ 函数，$u\left(x,t\right)$ 为待求解函数，在远离点源处足够光滑，$F\left(t\right)$ 为点源的时间依赖函数。

由于狄拉克$\delta$ 函数所引入的强奇异性，传统数值方法在处理这类问题时面临严峻挑战。针对这一问题，近年来，学者们研究出一些数值求解方法，如有限体积法[1]、有限差分法[2] [3]、有限元法[4]、移动网格方法[5]、间断伽辽金方法[6]等。浸入界面法[7]是修正差分格式的典型代表，它由Leveque和Li于1994年提出，利用跳跃条件进行局部格式的修正，为后续研究开辟了道路。之后，Zhao和Wei等人提出匹配界面和边界方法[8]，他们巧妙地引入虚拟点，并利用最低阶界面跳跃条件迭代确定虚拟值，构造了高精度的有限差分格式，从而实现对复杂几何问题的高精度模拟，在此基础上，学者们进行了大量的改进、优化与研究[9]-[11]。这些方法广泛应用于生物医学工程[12]、电磁学[13]、流体力学[14]等领域。然而，学者们对于带有点源的波动方程的研究较少。

本文剩余部分组织如下：在第2节，为带有固定点源的阻尼波动方程设计数值解法；在第3节，通过数值算例验证方法的有效性；在第4节，给出结论。

2. 数值离散格式

考虑计算区域$\left[x_a,x_b\right]\times \left[0,T\right]$ ，根据点源位置设置空间步长$h=\left(x_b-x_a\right)/N$ ，空间网格点 $x_i=x_a+ih\left(i=0,1,\cdots ,N\right)$ ，使点源位置$\overline{x}$ 恰好落在某个网格点上。设时间步长$\tau =T/M$ ，则时间节点为$t_n=n\tau \left(n=0,1,\cdots ,M\right)$ 。记$u\left(x_i,t_n\right)$ 为$u_i^n$ 。

取任意时刻的点源$\left(\overline{x},\overline{t}\right)$ ，对任意小的正数$\epsilon$ ，在$\left(\overline{x}-\epsilon ,\overline{x}+\epsilon \right)$ 上对(1)左右两边取积分，有

${\displaystyle {\int }_{\overline{x}-\epsilon }^{\overline{x}+\epsilon }\left(a{u}_{tt}+b{u}_t-c^2{u}_{xx}\right)\text{d}x}={\displaystyle {\int }_{\overline{x}-\epsilon }^{\overline{x}+\epsilon }\left(\delta \left(x-\overline{x}\right)F\left(\overline{t}\right)\right)\text{d}x}=F\left(\overline{t}\right),$(4)

即

$a{\displaystyle {\int }_{\overline{x}-\epsilon }^{\overline{x}+\epsilon }{u}_{tt}}\text{d}x+b{\displaystyle {\int }_{\overline{x}-\epsilon }^{\overline{x}+\epsilon }{u}_t}\text{d}x-c^2{\displaystyle {\int }_{\overline{x}-\epsilon }^{\overline{x}+\epsilon }{u}_{xx}}\text{d}x=F\left(\overline{t}\right),$ (5)

因为$u$ 连续，并且$u_t$ 和$u_{tt}$ 在$\overline{x}$ 附近有界，所以当$\epsilon \to 0^{+}$ 时

${\displaystyle {\int }_{\overline{x}-\epsilon }^{\overline{x}+\epsilon }{u}_t}\text{d}x\to 0,$(6)

${\displaystyle {\int }_{\overline{x}-\epsilon }^{\overline{x}+\epsilon }{u}_{tt}}\text{d}x\to 0,$(7)

记${\left[u\right]}_{(x,t)}=u\left(x+\epsilon ,t\right)-u\left(x-\epsilon ,t\right)$ ，进而

${\displaystyle {\int }_{\overline{x}-\epsilon }^{\overline{x}+\epsilon }{u}_{xx}}\text{d}x={u_x|}_{\overline{x}-\epsilon }^{\overline{x}+\epsilon }={\left[u_x\right]}_{(\overline{x},\overline{t})}=-\frac{F\left(\overline{t}\right)}{c^2}.$(8)

对(1)在$\left({\overline{x}}^{+},t\right)$ 与$\left({\overline{x}}^{-},t\right)$ 处取值作差，得到关系式

$a{\left[u_{tt}\right]}_{(\overline{x},\overline{t})}+b{\left[u_t\right]}_{(\overline{x},\overline{t})}-c^2{\left[u_{xx}\right]}_{(\overline{x},\overline{t})}=0,$ (9)

因为${\left[u_t\right]}_{(\overline{x},\overline{t})}=0,{\left[u_{tt}\right]}_{(\overline{x},\overline{t})}=0$，所以

${\left[u_{xx}\right]}_{(\overline{x},\overline{t})}=0.$ (10)

至此，我们已经得出了$u$ 的导数在$\overline{x}$ 处的“跳跃”值，接下来，我们将构建离散格式。记点源$\overline{x}$ 落在网格点$x_k$ 处，即$\overline{x}=x_k=x_a+kh$ 。取$\left(x,t\right)=\left(x_i^{-},t_n\right)$ ，有

$a{u}_{tt}\left(x_i^{-},t_n\right)+b{u}_t\left(x_i^{-},t_n\right)-c^2{u}_{xx}\left(x_i^{-},t_n\right)=0,$ (11)

每一项取近似

$u_{tt}\left(x_i^{-},t_n\right)\approx {\left(u_{tt}\right)}_i^n,u_t\left(x_i^{-},t_n\right)\approx {\left(u_t\right)}_i^n,u_{xx}\left(x_i^{-},t_n\right)\approx {\left(u_{xx}\right)}_i^n,$ (12)

可得离散格式

$a{\left(u_{tt}\right)}_i^n+b{\left(u_t\right)}_i^n-c^2{\left(u_{xx}\right)}_i^n=0.$ (13)

对于${\left(u_t\right)}_i^n$ 和${\left(u_{tt}\right)}_i^n$ ：

${\left(u_t\right)}_i^n=\frac{u_i^n-u_i^{n-1}}{\tau },$ (14)

${\left(u_{tt}\right)}_i^n=\{\begin{array}{l}\frac{2\left(u_i^1-u_i^0-\tau v_0\left(x_i\right)\right)}{{\tau }^2},n=1,\hfill \\ \frac{u_i^n-2u_i^{n-1}+u_i^{n-2}}{{\tau }^2},n>1.\hfill \end{array}$ (15)

受点源影响，解$u$ 的导函数$u_x$ 在点源$x=\overline{x}$ 处存在间断，在远离点源处仍保持其原有的光滑性。因此，${(u_{xx})}_i^n$ 的离散格式有以下几种情况：

当$i\ne k$ 时，解$u$ 在区间$\left[x_{i-1},x_{i+1}\right]$ 上足够光滑，有

${\left(u_{xx}\right)}_i^n=\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{h^2},$ (16)

当$i=k$ 时，构造函数

$\tilde{u}\left(x,t_n\right)=\{\begin{array}{l}u\left(x,t_n\right)+\left(x-x_i\right){\left[u_x\right]}_{(\overline{x},t_n)}+\frac{1}{2}{\left(x-x_i\right)}^2{\left[u_{xx}\right]}_{(\overline{x},t_n)},x\in \left[x_{i-1},x_i\right),\hfill \\ u\left(x,t_n\right),x\in \left[x_i,x_{i+1}\right].\hfill \end{array}$(17)

可知$\tilde{u}\left(x,\overline{t}\right)\in C^2\left[x_{i-1},x_{i+1}\right]$ ，并且$\tilde{u}\left(x,\overline{t}\right)$ 对于$x$ 的三阶导数有界。因此

$\begin{array}{l}{\left(u_{xx}\right)}_i^n=u_{xx}\left(x_i,t_n\right)={\tilde{u}}_{xx}\left(x_i,t_n\right)=\frac{\tilde{u}\left(x_{i+1},t_n\right)-2\tilde{u}\left(x_i,t_n\right)+\tilde{u}\left(x_{i-1},t_n\right)}{h^2}\\ =\frac{u\left(x_{i+1},t_n\right)-2u\left(x_i,t_n\right)+u\left(x_{i-1},t_n\right)}{h^2}+\frac{-2{\left[u_x\right]}_{(\overline{x},t_n)}+h{\left[u_{xx}\right]}_{(\overline{x},t_n)}}{2h}.\end{array}$(18)

值得说明的是，本文采用的数值格式为隐式格式，格式是无条件稳定的，其数值稳定性不依赖于时间步长和空间步长的比值，即不受经典显式格式中CFL条件(Courant-Friedrichs-Lewy condition)的限制。但在实际计算时，为保证数值精度和物理上的合理性，仍需选择适当的步长。

3. 数值实验

考虑带有固定点源的阻尼波动方程

$\{\begin{array}{l}a{u}_{tt}+b{u}_t-c^2{u}_{xx}=\delta \left(x-\overline{x}\right)F\left(t\right),-\pi <x<\pi ,t>0,\hfill \\ u\left(x,0\right)=u_0\left(x\right)=\sin x,\hfill \\ u_t\left(x,0\right)=0,\hfill \end{array}$ (19)

具有周期边界条件，$F\left(t\right)=10t+10$ 是随时间线性增长的点源(其中常数项表示某种恒定的激励，线性项表示激励随时间线性增强的大小，模拟物理系统中外力或能量不断输入的过程)，固定点源$\overline{x}=\pi /4$ 。取$a=1,b=1,c=4$ ，运用上述数值方法，可以得到解的图像如图1，其中固定点源影响情况见图2。

Figure 1. Numerical solution for a = 1, b = 1, c = 4

图1. a = 1，b = 1，c = 4时的数值解

Figure 2. The influence of fixed point source when a = 1, b = 1, c = 4

图2. a = 1，b = 1，c = 4时固定点源影响情况

$t=1$ 时刻整个空间区域上解的分布见图3，其中红色竖线表示点源位置。

Figure 3. The distribution of solutions at time t = 1

图3. t = 1时刻解的分布

时间份数固定为3200时误差及收敛阶随空间份数变化情况见表1。

Table 1. The error and convergence order when the time fraction is fixed to 3200

表1. 时间份数固定为3200时的误差及收敛阶

空间份数	误差	收敛阶
16	1.439544 × 10^−^2	−
32	4.790466 × 10^−^3	1.5874
64	1.509366 × 10^−^3	1.6662
128	4.104785 × 10^−^4	1.8786
256	1.007884 × 10^−^4	2.0260
512	2.483341 × 10^−^5	2.0210
1024	6.181888 × 10^−^6	2.0062

空间份数固定为100时误差及收敛阶随时间份数变化情况见表2。

Table 2. The error and convergence order when the space fraction is fixed to 100

表2. 空间份数固定为100时的误差及收敛阶

时间份数	误差	收敛阶
100	2.121945 × 10^−^2	−
200	1.081448 × 10^−^2	0.9724
400	5.459768 × 10^−^3	0.9861
800	2.743270 × 10^−^3	0.9929
1600	1.375036 × 10^−^3	0.9964
3200	6.883807 × 10^−^4	0.9982
6400	3.444088 × 10^−^4	0.9991

4. 结论

本文在均匀网格上对带有固定点源的阻尼波动方程进行了离散，通过构造辅助函数有效处理了狄拉克$\delta$ 函数所引入的奇异性，数值实验则验证了格式的有效性，结果表明，我们的方法在空间上具有二阶收敛精度，在时间上具有一阶收敛精度。然而，当点源变为移动点源时，求解将更加困难，下一步，可以在均匀网格上探索带有移动点源的阻尼波动方程的数值解法。