一种用于高频和半无限贝塞尔积分的新方法

《Applied Mathematics and Computation》：A novel method for high-frequency and semi-infinite Bessel integrals

　　陈瑞云|洪都，中国广东省湛江市广东海洋大学数学与计算机学院，邮编524088

摘要
本文提出了一种用于高频半无限贝塞尔积分的新方法。我们将积分区间分解为有限部分和尾部部分，然后分别利用数值最速下降法高效地近似这两部分。通过理论分析及数值算例展示了该方法的效率。

引言
本文介绍了一种用于近似包含高频贝塞尔函数的半无限积分的新的数值方法。这类积分定义为：
$$I[f;ω]=∫_0^∞f(x)Ju(ωx)dx$$
其中$u$为非负整数，$f(x)$在$[0, ∞)$上具有足够的平滑性，$ω?1$表示较大的频率参数，$Ju(·)$表示阶数为$u$的第一类贝塞尔函数。这类积分在数学物理和工程学的诸多领域中起着重要作用[1]、[2]、[3]、[4]。一个典型的应用场景是固体圆柱体的热传导问题，其柱坐标下的热方程为：
$$?T/?t=α(?^2T/?r^2+1/r?T/?r)，0≤r0$$
其中$a$为固体圆柱体的半径，$T(r, t)$为温度分布，$α=kρcp$为热扩散率，$k$为热导率，$\rho$为密度，$cp$为恒压比热容。该方程满足边界条件和初始条件：
$$T(a,t)=0，t>0；|T(0,t)|<∞；T(r,0)=f(r)，0≤r≤a$$
对于函数$T(r, t)$，其有限汉克尔变换及其逆变换定义为：
$$T?(λn,t)=∫_0^a rT(r,t)J_0(λr)dr；T(r,t)=2a^2∑_{n=1}^∞J_0(λnr)[J_1(λna)]^2T?(λn,t)$$
其中$λn=αna$，$αn$是$J_0(x)=0$的正根。根据[5]，方程(1.2)的解为：
$$T(r,t)=2a^2∑_{n=1}^∞J_0(λnr)[J_1(λna)]^2e^{-αλn^2t}∫_0^a f(r)J_0(λnr)dr$$
尽管包括辛普森法则和高斯求积法在内的经典求积规则能够可靠地计算形如$∫_a^b f(x)dx$的积分，但对于形如$∫_a^b f(x)e^{iωg(x)}dx$的积分，由于存在剧烈振荡，这些方法的计算效果较差。为解决这一难题，近年来人们对这类振荡积分的数值近似方法进行了大量研究，由此发展出了多种高效方法，包括菲隆型方法[6]、[7]、莱文型方法[8]、[9]以及数值最速下降法[10]、[11]、[12]。
菲隆型方法用多项式插值函数替代振幅$f(x)$（如在端点及驻点处），并精确计算振荡部分，当有$m$个插值点时，其误差可达到$O(ω^{-m-1})$阶。莱文型方法将积分转化为常微分方程$(F(x)e^{iωg(x)})′=f(x)e^{iωg(x)}$，再通过配置法求解$F(x)$，从而避免显式计算矩值，并能有效处理任意相位函数。数值最速下降法则将积分路径变形到复平面上，沿着最速下降路径行进，消除振荡，还将积分转化为一系列快速衰减的积分，便于使用标准高斯求积法进行计算。
在有限振荡积分的研究中，也发展出了专门的技术，包括复变积分法[10]和数值最速下降法[11]。复变积分法适用于形如$∫_{-1}^1 f(x)e^{iωx}dx$、$∫_0^∞ f(x)e^{iωx}dx$的积分，其中$f(x)$为解析实值函数。该方法将原始积分路径从实轴变形到复平面上，通过选择一条避开奇点且振荡核呈指数衰减的合适路径，将原本具有剧烈振荡或奇点的积分转化为一系列无振荡积分，这类积分非常适合使用标准求积规则，如高斯求积法。数值最速下降法则是针对形如$∫_a^b f(x)e^{iωg(x)}dx$的剧烈振荡傅里叶型积分设计的，其中$f(x)$和$g(x)$均为解析函数，且$g(x)$在$[a, b]$上为实值函数。该方法是对渐近分析中经典最速下降原理的数值化改进，其核心思路是利用柯西定理将积分路径变形到复平面上，沿着$g(x)$的虚部保持恒定的路径行进。在这些路径上，振荡因子$e^{iωg(x)}$会呈指数衰减，从而使原始积分转化为拉普拉斯型积分，此类积分便可使用高斯-拉盖尔求积法准确计算。这两种方法都采用了在复平面上变形积分路径以消除振荡的策略，但适用范围有所不同：米洛瓦诺维奇的方法适用于包含振荡核和奇点核的更广泛类别的积分，而许布雷赫斯和范德瓦勒则专注于具有解析相位函数的傅里叶型积分，并为最速下降原理提供了系统的数值实现方法。
在[13]中，我们得到了一个重要恒等式：
$$Ju(ωx)=iu+\frac{1}{π}K_u(iωx)+(-i)\frac{1}{π}K_u(-iωx)$$
其中$K_u(x)$表示第二类修正贝塞尔函数。在米洛瓦诺维奇和许布雷赫斯[10]、[11]的工作基础上，我们还将数值最速下降法扩展到了处理有限区间上的积分[12]、[13]、[14]、[15]，即$∫_a^b f(x)J_v(ωx)dx$。关于将该方法进一步应用于类型为(1.9)的积分的相关研究，可参考[16]、[17]、[18]、[19]、[20]。
本文将数值最速下降法应用于高频半无限贝塞尔积分。计算积分(1.1)时面临两个主要挑战：$x=0$处的奇点以及贝塞尔函数的高振荡性。为解决奇点问题，我们通过减去$f(x)$在$x=0$处的$2n+u$阶泰勒展开式来处理。针对振荡问题，我们利用恒等式(1.8)将贝塞尔积分转化为傅里叶积分，进而引入数值最速下降法。
本文的结构如下：第2节详细介绍了如何通过超几何函数计算矩值$\rho_k(ω)$，并阐述了所提方法的构建过程。第3节对所提方法的误差进行了分析。第4节通过数值算例验证了该方法的有效性。

数值方法的发展
我们首先将原始积分(1.1)分解为三个部分：
$$I[f;ω]=I_1[H;ω]+I_2[f;ω]+∑_{k=0}^{2n+u}uf(k)(0)k!ρ_k(ω)$$
其中前两个部分分别是有限积分和半无限积分，定义如下：
$$I_1[f;ω]=∫_0^1H(x)Ju(ωx)dx；I_2[f;ω]=∫_1^∞f(x)Ju(ωx)dx$$
此处函数$H(x)$是通过减去$f(x)$在$x=0$处的$2n+u$阶泰勒展开式来定义的：
$$H(x)=f(x)?∑_{k=0}^{2n+u}uf(k)(0)k!x^k$$
这种减法操作消除了被积函数在$x=0$处的奇点。辅助矩值$\rho_k(ω)$则...

误差分析
在介绍了数值方法之后，我们现在对其误差行为进行分析。
定理3.1
在区域${0≤Re(z)≤1}$内，$f(x)$为解析函数，且满足$|f(p±iR)|≤MR^α$，当$R→∞$时，其中$\alpha<3/2$，$M＞0$为常数。那么，当$ω?1$时，数值方法(2.22)用于近似$I_1[f; ω]$的误差满足：
$$|I_1[f;ω]?Q_{1,n}NSD[f;ω]|=O((n!)^2(2n)!ω^{-2n-1})$$
由于该证明与[13]中的定理3类似，此处不再详述。
定理3.2
在与前一个定理相同的解析性和增长条件下，...

数值算例
我们通过数值实验来验证理论分析的结果，同时展示公式(2.37)在计算高频半无限贝塞尔积分时的高效性。性能评估指标包括计算精度和收敛阶数。
在以下结果中，“AE”表示绝对误差，“RE”表示相对误差。所有实验均在Unis?14B?0023环境中实现。
示例4.1
考虑积分：
$$I[f;ω]=∫_0^∞\frac{1}{1+x^2}J_2(ωx)dx$$
示例4.2
考虑积分：
$$I[f;ω]=∫_1^∞\frac{x}{1+x^2}J_2(ωx)dx$$

结论
本文针对一类包含高频贝塞尔函数的积分提出了新方法，推导出了绝对误差与ω的依赖关系。未来的工作将致力于解决包含其他高频贝塞尔型函数的更一般积分的计算难题。