物理和工程领域的许多问题,往往可以由偏微分方程来表示,而这些偏微分方程有时可以等价的转换为积分方程^[6-7].特别地,在二维声散射问题中,时谐波的声散射问题可以用如下满足一定条件的Helmholtz方程表示,参见文献[3, 9, 20, 23]

(1.1) $\begin{equation}\label{eq:a1}\Delta u(x) + \omega^2u(x) = 0, \quad x \in \mathbb{R}^2\setminus \bar{\Omega}, \end{equation}$

而这个偏微分方程又可以等价的转换为如下的积分方程

(1.2) $\begin{equation}\label{BIE1}u(x)=\int_{\Gamma}u(y)\frac{\partialG(x, y)}{\partialn(y)}-G(x, y)\frac{\partialu}{\partialn}(y){\rm d}s(y), \quad x\in \mathbb{R}^2\setminus \bar{\Omega}, \end{equation}$

其中

(1.3) $\begin{equation}\label{H01}G(x, y)=\frac{\rm i}{4}H_0^{(1)}(\omega|x-y|), \end{equation}$

$H_0^{(1)}$表示第一类Hankel函数, $n(y)$表示$y$点处的外法向量, $\Gamma$表示障碍物$\Omega$的边界.

对于方程(1.1),经典的理论已经证明了解的存在性和唯一性,参见文献[2, 9].对于散射问题,需要再加上如下的Sommerfeld条件

$\lim\limits_{r \to \infty}r^{1/2}(\frac{\partial u^s}{\partial r}-{\rm i}\omega u^s)=0, \quad r:=|x|, $

其中$u^s:=u-u^i$表示散射场.

在过去几年,提出了许多求解边界积分方程(1.2)的数值方法,其中边界元方法是最有效的方法之一.边界元方法只需要对边界进行离散,这使得离散的网格生成速度更快.但是该方法将导致求解以非对称稠密矩阵为系数矩阵的线性方程组.如果用$N$表示自由度的个数,则需要存储$O(N^2)$元素.如采用直接方法求解这个方程组,将需要$O(N^3)$量级的运算量.为了克服这个困难, Rokhlin和Greengard提出了快速多极方法,参见文献[4-6].该方法一经提出便被广泛应用于求解大规模的工程问题如势问题、弹性力学、声散射问题等.然而,函数(1.3)在零点是弱奇异的并且当$\omega\gg 1$时,该函数还是高振荡的.因此经典的边界元方法或快速多极方法求解这类问题时会遇到困难.

另一方面,当$\omega$很大时,方程(1.2)的解也是高振荡的,因此自由度的个数至少要与$\omega$线性成比例才能保证解的准确性.近几年, Langdon, Chandler-Wilde等将解的高频渐近性质融入到逼近空间以降低计算量.当$\omega\rightarrow \infty$时,在远离边界中的角和阴影部分,该方法用如下函数逼近$\frac{1}{\omega}\frac{\partial u}{\partial n}$

(1.4) $\begin{eqnarray}\label{Asy1}\Psi :=\left\{\begin{array}{ll} \frac{2}{\omega}\frac{\partial u^i}{\partial n}, \; & \mbox{in the illuminated region }(n\cdot d<0), \\ 0, & \mbox{in the shadow region }(n\cdot d>0).\end{array}\right. \end{eqnarray}$

对于有尖角的边界如凸多边形, $\frac{1}{\omega}\frac{\partial u}{\partial n}$采用如下函数表示

(1.5) $\begin{eqnarray}\label{bijjie} \frac{1}{\omega}\frac{\partial u}{\partial n}(x)=\Psi(x) + \omega \sum\limits_{m=1}^{n} [{\rm e}^{{\rm i}\omega x\cdot d_m}v_m^+(x, \omega)+{\rm e}^{-{\rm i}\omega x\cdot d_m}v_m^-(x, \omega)], x\in\Gamma, \end{eqnarray}$

其中$n$表示多边形的边数, $d_m$是与第$m$条边$\Gamma_m$平行的单位向量.理论分析和数值结果显示该方法保证精度所需要的自由度只需要$O(\omega^{1/3})$量级.再后来通过网格改进,保证精度所需的自由度的量级只需要$O(\omega^{1/9})$量级,到最后甚至与$\omega$无关,参见文献[21].这些研究结果对高频散射问题中高振荡积方程求解的数值方法的设计非常重要.

本文,将考虑如下的高振荡积分方程

(1.6) $\begin{equation}\label{BIE2}{\displaystyle\frac{1}{2}q(x)+\int_{\Gamma}\left(\frac{\partialG(x, y)}{\partialn(x)}+{\rm i}\eta G(x, y)\right) q(y){\rm d}s(y)=f(x), \quad x\in \Gamma, }\end{equation}$

其中$q(x):=\frac{\partial u}{\partial n}(x)$是未知函数, $f(x):={\rm i}\eta u^i+\partial u^i/\partial n$是已知函数.事实上,该方程是散射问题中当边界条件为$u(x)=0, x \in \Gamma$时的Burton-Miller形式.我们采用常数元离散积分方程,利用快速多极方法对算法进行加速,所涉及的高振荡积分通过CCF方法^[17-19]进行计算.

本节将介绍快速多极方法的一些基本公式.为了数值求解边界积分方程,需要对方程进行离散化.假设边界$\Gamma$可以分为$\Gamma_m, m=1, \cdots, N$,将未知函数假设成(1.5)式的形式,则边界积分方程(1.6)可写为

(2.1) $\begin{eqnarray}\label{BIE2j} \nonumber &&\sum\limits_{m}\int_{\Gamma_m} \left(2\frac{\partial G(x, y)}{\partial n(x)}-2{\rm i}\omega G(x, y)\right)(\exp({\rm i}\omega x\cdot d_m)v_{m}^+ +\exp(-{\rm i}\omega x\cdot d_m)v^{-}_m){\rm d}s(y)\\ &&+\sum\limits_{m}\exp({\rm i}\omega x\cdot d_m)v_m^++\exp(-{\rm i}\omega x\cdot d_m)v_m^- \nonumber\\ &=&2\frac{\partial u^i}{\partial n}\frac{1}{\omega}-2{\rm i}u^i-\int_{\Gamma}\left(2\frac{G(x, y)}{n(x)}-2{\rm i}\omega G(x, y)\right)\Psi {\rm d}s(y)-\Psi.\end{eqnarray}$

然后再利用常数元逼近$v_m^+, v_m^-$.

对于多边形散射体,利用类似文献[10, 20]中的复化网格,该网格在$[0, \lambda]$和$[\lambda, L]$采用不同的网格,其中$\lambda=2\pi/\omega$.

定义 2.1   对于$L> \lambda> 0, q>0, N=2, 3, \cdots$,网格$I_{N, \hat{N}}:=\{y_0, \cdots, y_{N+\hat{N}}\}$由如下点构成

$ y_i=\lambda(\frac{\rm i}{N})^q, \quad i=0, \cdots, N, $

$y_{N+j}:=\lambda(\frac{L}{\lambda})^{\lambda/\hat{N}}, \quad j=1, \cdots, \hat{N}, $

其中$\hat{N}$是大于或等于$N^*$的最小整数,并且$N^*=\frac{-\log(L/\lambda)}{q\log(1-1/N)}.$

这个复化网格在$[0, \lambda]$中有$N$个点,在$[\lambda, L]$中有$\hat{N}$个点.利用该网格的好处可以参看文献[18].

类似于文献[18],定义$q_j:=(2v+3)/(2\pi/\Omega_j-1), j=1, \cdots, n$,其中$\Omega_j\in (\pi, 2\pi)$表示顶点$P_j$处的外角.在多边形的每一边,取每个区间的中点作为配置点.

将配置点代入方程(2.1),可得一个线性方程组${\rm{A}}{\bf{q}} = {\bf{b}}$.若用高斯消去法求解该方程组需要$O(N^3)$量级的运算量,若用迭代方法求解由于需要做矩阵向量的乘积,也需要$O(N^2)$量级的运算量.通过对核函数的多极展开,并利用源点和场点之间的关系,快速多极方法可以加速矩阵向量的乘积运算.

接下来,介绍快速多极对角展开的一些基本公式^[15].为了得到核函数$G(x, y)$的多极展开,我们先介绍如下的引理和公式^[1].

引理 2.1   假设极坐标下$x=(r_x, \theta_x)$, $y=(r_y, \theta_y)$其中$r_x>r_y$,则

(2.2) $\begin{equation}\label{L1} H_v(\omega \|x-y\|){\rm e}^{\pm v \theta_{x-y}}=\sum\limits_{m=-\infty}^{+\infty}H_{v+m}(\omega r_x){\rm e}^{\pm (v+m) \theta_x}J_{m}(\omega r_y){\rm e}^{\mp m\theta_y}, \end{equation} $

(2.3) $\begin{equation} J_v(\omega \|x-y\|){\rm e}^{\pm v \theta_{x-y}}=\sum\limits_{m=-\infty}^{+\infty}J_{v+m}(\omega r_x){\rm e}^{\pm (v+m) \theta_x}J_{m}(\omega r_y){\rm e}^{\mp m\theta_y}, \end{equation}$

其中$m$为整数, $H_v(x), J_{m}(x)$分别为第一类汉克函数和贝塞尔函数, $\theta_{x-y}$是向量$x-y$与$x$轴的夹角.

对于贝塞尔函数$J_m( x)$,有如下的积分等式,参见文献[1].

(2.4) $\begin{equation}\label{J1} J_m( x)=\frac{1}{2\pi {\rm i}^m}\int_0^{2\pi}{\rm e}^{{\rm i}(x\cos(t)-mt)}{\rm d}t, \end{equation}$

其中i$:=\sqrt{-1}.$

由(2.2)和(2.4)式,当$\mid\parallel x-x_l\parallel-\parallel x_c-y\parallel\mid<\parallel x_l-x_c\parallel$时, $G(x, y)$可写成如下形式^[15-16]

(2.5) $\begin{equation} G(x, y)\approx \frac{\rm i}{4L}\textbf{h}_{\omega}(x, x_l)T_{\omega}(x_l, x_c)\textbf{g}_{\omega}(x_c, y), \end{equation}$

其中$\textbf{h}_{\omega}(x, x_l)=[h_{\omega}^0, \cdots, h_{\omega}^{L-1}]$, $T_{\omega}={\rm diag}(t_{\omega}^0, \cdots, t_{\omega}^{L-1})$, $\textbf{g}_{\omega}=[g_{\omega}^0, \cdots, g_{\omega}^{L-1}]^T$,其元素定义如下

(2.6) $ \begin{equation} h_{\omega}^k(x, x_l)={\rm e}^{{\rm i}\omega (x-x_l)\cdot S(\beta_k)}\label{h}, \end{equation} $

(2.7) $\begin{equation} g_{\omega}^k(x_c, y)={\rm e}^{-{\rm i}\omega (x_c-y)\cdot S(\beta_k)}\label{g}, \end{equation} $

(2.8) $\begin{equation} t_{\omega}^k(x_l, x_c)=\sum\limits_{j=-M}^M {\rm i}^{-j}H_{j}(\omega\|x_c-x_l\|){\rm e}^{{\rm i}j(\theta_{x_c-x_l}-\beta_k)}\label{T}, \end{equation}$

其中$S(\beta_k):=(\cos(\beta_k), \sin(\beta_k))^T, \quad \beta_k := 2\pi k/L, \quad k=0, \cdots, L-1$.

公式(2.6)和(2.8)中整数$L$表示截断项(2.2)式, $M$是用梯形公式离散积分(2.4)时点的个数.根据经验一般取$M=\omega R+5\log(\omega R+\pi)$,其中$R$是当前层的最大半径. $L$一般取为$2M+1$,参见文献[14-15].方程(1.6)中的积分可以写成

(2.9) $\begin{eqnarray}\nonumber \int_{\Gamma_j} G(x_k, y)q(y){\rm d}s(y)&\approx &\int_{\Gamma_j}\frac{\rm i}{4L}\textbf{h}_{\omega}(x_k, x_l)T_{\omega}(x_l, x_c)\textbf{g}_{\omega}(x_c, y)q(y){\rm d}s(y)\\&=&\frac{\rm i}{4L}\textbf{h}_{\omega}(x_k, x_l)T_{\omega}(x_l, x_c)M_{\omega}(x_c), \label{M1}\end{eqnarray}$

其中

$\begin{eqnarray*}M_{\omega} (x_c)=[M_{\omega}^0 (x_c), \cdots, M_{\omega}^{L-1} (x_c)]^T , \quad M_{\omega}^k (x_c)=\int_{\Gamma_j}g_{\omega}^k(x_c, y)q(y){\rm d}s(y), \end{eqnarray*}$

并称之为$x_c$点处的矩,且该矩是与$x_k$点无关的.如果展开点由$x_c$移动到$x_{c'}$,有如下的转移公式

$g_{\omega}^k(x_{c'}, y)={\rm e}^{-{\rm i}\omega (x_{c'}-y)\cdot S(\beta_k)} = {\rm e}^{-{\rm i}\omega (x_{c'}-x_c)}g_{\omega}^k(x_{c}, y), $

并且可得

(2.10) $\begin{equation} \label{M2M}M_{\omega}^k (x_{c'})={\rm e}^{-{\rm i}\omega (x_{c'}-x_c)}M_{\omega}^k (x_{c}).\end{equation} $

该公式称为矩-矩的转移公式,记为(M2M).

另一方面, (2.9)式可写成

$\int_{\Gamma_j} G(x_k, y)q(y){\rm d}s(y)\approx\frac{\rm i}{4L}\textbf{h}_{\omega}(x_k, x_l)L_{\omega}(x_l), \label{L2}$

其中

(2.11) $\begin{equation} L_{\omega}(x_l)=T_{\omega}(x_l, x_c)M_{\omega}(x_c), \end{equation}$

该公式称为矩-局部的转移公式,记为(M2L).

如果局部点由$x_l$移动到$x_{l'}$,可得如下的转移公式

$h_{\omega}^k(x, x_{l'})={\rm e}^{{\rm i}\omega (x-x_{l'})\cdot S(\beta_k)} = {\rm e}^{{\rm i}\omega (x_{l}-x_{l'})}h_{\omega}^k(x, x_l), $

由此可得

(2.12) $\begin{equation} \label{L2L}L_{\omega} (x_{l'})={\rm e}^{{\rm i}\omega (x_{l}-x_{l'})}L_{\omega} (x_{l}), \end{equation} $

该公式称为局部-局部的转移公式(L2L).

对于积分

$\begin{eqnarray*}\int_{\Gamma_j}\frac{\partialG(x_k, y)}{\partialn(x)} q(y){\rm d}s(y)&\approx &\int_{\Gamma_j}\frac{\rm i}{4L}\frac{\partial \textbf{h}_{\omega}(x_k, x_l)}{\partial n(x)}T_{\omega}(x_l, x_c)\textbf{g}_{\omega}(x_c, y){\rm d}s(y)\\&=&\frac{\rm i}{4L}\frac{\partial \textbf{h}_{\omega}(x_k, x_l)}{\partial n(x)}T_{\omega}(x_l, x_c)M_{\omega}(x_c), \end{eqnarray*}$

矩, M2M, M2L和L2L转移公式与$G(x_k, y)$的相同且$\textbf{h}_{\omega}(x, x_l)$的导数具有如下的形式

$\frac{\partial h_{\omega}^k(x, x_l)}{\partial n(x)}={\rm i}\omega n(x)\cdot S(\beta_k) h_{\omega}^k(x, x_l).$

本节介绍计算高振荡积分的Clenshaw-Curtis Filon方法,假设边界$\Gamma$有如下的参数方程形式

$ x(s)=(x_1(s), x_2(s)), y(t)=(y_1(t), y_2(t)) \quad s, t\in [a, b], $

${\rm d}s(y)=\sqrt{y_1'(t)^2+y_2'(t)^2}.$

则多极展开公式(2.9)中的

$M_{\omega}^k (x_c)=\int_{\Gamma_j} {\rm e}^{-{\rm i}\omega (x_c-y)\cdot S(\beta_k)}q(y){\rm d}s(y)$

可以表示为

$ M_{\omega}^k (x_c)=\int_{\Gamma_j} {\rm e}^{-{\rm i}\omega ((c_1-y_1(t))\cos\beta_k+(c_2-y_2(t))\sin\beta_k)}\hat{q}(t)\sqrt{y_1'(t)^2+y_2'(t)^2}{\rm d}t, $

其中$x_c=(c_1, c_2), \hat{q}(t)=q(y(t))$.如果通过(1.5)式来近似$q(t)$,可得如下一般形式的积分

(3.1) $ \begin{equation} \label{Int2} I[f]=\int_a^b f(x){\rm e}^{{\rm i}\omega g(x)} {\rm d}x. \end{equation} $

当$\omega\gg 1$时,被积函数是高振荡的,经典的数值积分法计算这类积分是很困难的,参看文献[12-13, 22].

对于一般的振荡因子$g(\tau)$,积分没有解析表达式,必须用数值方法进行计算. CCF方法是计算这类积分的重要方法.

计算积分$I[f]=\int_{-1}^1 f(x){\rm e}^{{\rm i}\omega x} {\rm d}x$的CCF方法:设$Q_N$在Clenshaw-Curtis点处的插值多项式,则$I[f]$通过下式近似

(3.2) $\begin{eqnarray}I_N[f]:=\int_{-1}^1 Q_N(x){\rm e}^{{\rm i}\omega x}{\rm d}x=\sum\limits_{n=0}^{N}{''}a_{n, N}(f)W_n, \end{eqnarray}$

其中对于$n\geq 0$, $W_n:=\int_{-1}^1T_n(x){\rm e}^{{\rm i}\omega x}{\rm d}x$, $T_n(x)=\cos(n\arccos(x))$,并且

$a_{n, N}(f)=\frac{2}{N}\sum\limits_{j=0}^{N}{''}\cos(\frac{jn\pi}{N})f(t_{j, N}), \quad n=0, \cdots, N, $

可以由FFT快速计算, $t_{j, N}$表示Clenshaw-Curtis点(定义为$t_{j, N}:=\cos(j\pi/N)$),符号$\sum{''}$表示第一项和最后一项乘以$\frac{1}{2}$.

利用等式$2T_n(x)=U_n(x)-U_{n-2}(x)$,其中$n\geq2$, $U_n(x)=\frac{1}{n+1}T'_{n+1}(x)$是第二类Chebyshev多项式,可得$W_n$的递推公式,即

(3.3) $\begin{equation}\label{Chebyrec} 2W_n=\rho_{n+1}-\rho_{n-1}, \quad n\geq 2, \end{equation}$

其中

$\rho_n:=\int_{-1}^1U_{n-1}(x){\rm e}^{{\rm i}\omega x}{\rm d}x.$

另一方面,对$W_n:=\int_{-1}^1T_n(x){\rm e}^{{\rm i}\omega x}{\rm d}x$进行分部积分,可得

(3.4) $\begin{equation}\label{ditui} W_0=\sigma_0 , \quad W_n=\sigma_n-\frac{n}{{\rm i}\omega}\rho_{n}, \quad n\geq 1, \end{equation}$

其中

${\sigma _n} = \frac{1}{{{\rm{i}}\omega }}{T_n}(x){{\rm{e}}^{{\rm{i}}\omega x}}\left| {\begin{array}{*{20}{c}}{x = 1}\\{x = - 1}\end{array}} \right. = \left\{ {\begin{array}{*{20}{c}}{\frac{{2\sin \omega }}{\omega },{\rm{for}}\;{\rm{even}}\quad n,}\\{\frac{{2\cos \omega }}{{{\rm{i}}\omega }},{\rm{for}}\;{\rm{odd}}\quad n,}\end{array}} \right.$

结合(3.3)和(3.4)式,可得如下递推公式(参看文献[11])

(3.5) $\begin{equation}\label{gongshi} 2\sigma_n-\frac{2n}{{\rm i}\omega}\rho_n=\rho_{n+1}-\rho_{n-1}, \quad n\geq 2.\end{equation}$

更进一步有

$U_0(x) = 1 , \; U_1(x) = 2x , \; \rho_1 = \sigma_0 , \; \rho_2 = 2\sigma_1-\frac{2}{{\rm i}\omega}\sigma_0.$

对于一般的区间$[a, b]$,首先利用如下的变量替换

$x=p+qt, \quad t\in[-1, 1], \quad p:=\frac{b+a}{2}, \quad q:=\frac{b-a}{2}, $

将积分化为$I[f]:=\int_a^bf(x){\rm e}^{{\rm i}\omega x}{\rm d}x=q\int_{-1}^1f(p+qt){\rm e}^{{\rm i}\omega(p+qt)}{\rm d}t$.然后再利用CCF方法.

对于具有一般振荡因子$g(x)$的积分(3.1),为了简单起见,假设$g(x)\in C^{\infty}[a, b]$:

$\bullet$如果$g{'}(x)\neq 0$, $x\in[a, b]$,利用变量替换$\tau=g(x) , $则$g^{-1}(\tau)\in C^{\infty}$,且

(3.6) $ \begin{eqnarray} I^{[a, b]}[f]=I^{[g(a), g(b)]}[F], \; \mbox{其中}\ F=(f\circ g^{-1})|(g^{-1})'|, \end{eqnarray} $

则(3.1)式可以用CCF方法计算.

$\bullet$如果存在$\xi\in [a, b]$使得$g'(\xi)=g''(\xi)= \cdots = g^{(r)}(\xi)=0$, $g^{(r+1)}(\xi) \neq 0$且$g{'}(x)\neq 0$, $x\neq \xi$.不失一般性,假设$g^{(r+1)}(\xi)>0$.在每个区间上利用变量替换可得

$I^{[a, b]}[f]=\left( \int_{g(a)}^{g(\xi)}+\int_{g(\xi)}^{g(b)}\right) F(\tau){\rm e}^{{\rm i}\omega \tau}{\rm d}\tau, $

则可以用CCF方法计算.

对于CCF方法,其详细的研究可参看文献[17-19, 24-25],接下来我们给出几个算例来说明CCF方法的有效性.

算例1   考虑积分$I_1[f]=\int_{0.5}^1 \cos(x) {\rm e}^{{\rm i}\omega (x^2-x)}{\rm d}x$,其中$\xi=0.5$是$1$阶驻点. $m+1$为插值点的个数.

注 3.1   由于存在驻点,函数$F=(f\circ g^{-1})| (g^{-1})'|$在$\xi$是弱奇异的,因此将需要在一个复化网格上运用CCF方法,复化的网格定义为$\{(\frac{j}{M})^q, \ j=0, \cdots, M\}$.我们采用与文献[18]相同的策略计算两个小区间上的积分.

算例 2   考虑积分$I_2[f]=\int_{0}^1{\rm e}^x {\rm e}^{{\rm i}\omega cos(x)}{\rm d}x$,其中$\xi=0$是$1$阶驻点.

由表 1和表 2可知CCF方法是很有效的.

本节通过数值例子来验证方法的精确性和有效性,所有的计算都是在Matlab2014a上进行.考虑一个障碍物为边长为1的正方形的散射问题,入射角为$\pi/4$.将本文的方法与文献[10]中用到的经典配置方法进行比较.

记$Cu_N^{\omega}$经典配置方法得到的解, $u_N^{\omega}$为本文方法计算得到的解,计算结果见表 3.

由表 3可知当$\omega=100, 160, 200$时,本文的快速多极方法的精度可以与经典配置方法相匹敌.当$N=8, 16$时,本文方法的计算时间比经典方法要少很多.由于对角形式的快速多极方法对于低频问题不适用,所以当$\omega=10$时,方法的精度有所下降.如果不参考散射体的尺寸,单独因为$\omega$很小,得出问题是是低频是没有意义的.因此考虑用$\omega\times L$来刻画问题的振荡性,其中$L$是散射体的尺寸.由于采用复化网格,所以当$N$增加时,配置点将集中在正方形的四个角附近,尽管$\omega$很大,但是$\omega \times l$可能很小,其中$l$表示配置点所在小区间的长度.这就导致了一个低频问题,使得对角快速多极方法的精度降低,详细的误差分析可以参看文献[8].

对于具有高振荡核的积分方程,经典的数值积分法不适用,并且当频率增加时,计算代价增长很快.本文结合CCF方法、快速多极方法、配置方法,研究一类带有Hankel核函数的边界积分法的数值方法.数值算例显示该方法的精度与传统方法相当,但是计算时间要少很多.