高振荡Bessel变换的新型Hermite卷积积分公式及其应用
Novel Hermite Convolution Quadrature for a Class of Highly Oscillatory Bessel Transforms and Its Application
DOI: 10.12677/AAM.2022.114247, PDF, HTML, XML, 下载: 337  浏览: 531 
作者: 任 浩*, 曾 港, 李恒杰:贵州大学数学与统计学院,贵州 贵阳
关键词: 高振荡积分Hermite插值卷积积分Bessel核Highly Oscillatory Integrals Hermite Interpolation Convolution Integral Bessel Kernel
摘要: 高振荡问题已广泛出现于许多科学工程应用领域,例如电磁、声波散射、量子化学和图像分析等问题。高振荡数值积分是高振荡问题数值解中的重要研究方向,由于其被积函数具有高振荡性,传统数值积分解法面临许多挑战。基于Lubich的卷积求积,本文针对一类高振荡的Bessel变换提出了一种新型的Hermite卷积积分公式,并研究了其在高振荡Volterra积分方程的数值解中的应用。通过理论与数值实验表明,该方法在计算含高振荡Bessel核的卷积积分以及积分方程时,计算精度不受振荡频率的影响,是一种高效的计算方法。
Abstract: Highly oscillatory problems frequently arise in many science and engineering fields, such as electromagnetic, acoustic scattering, quantum chemistry and image analysis problems. Numerical calculation of highly oscillatory integration is an important research interest of numerical studies on highly oscillatory problems. Due to the high oscillation of the integrand, classical numerical integration encounters great challenges. Based on Lubich’s convolution quadrature, this paper proposes novel Hermite convolution quadrature for a class of highly oscillatory Bessel transforms and studies its application in the numerical solutions to highly oscillatory Volterra integral equations. Both theoretical and numerical evidences indicate that the proposed method is particularly efficient in computing highly oscillatory integrals and solving highly oscillatory integral equations, and its accuracy is not affected by the oscillation parameter.
文章引用:任浩, 曾港, 李恒杰. 高振荡Bessel变换的新型Hermite卷积积分公式及其应用[J]. 应用数学进展, 2022, 11(4): 2352-2361. https://doi.org/10.12677/AAM.2022.114247

1. 引言

在很多应用领域,高振荡积分数值计算问题是一个非常重要的研究课题。在量子化学、图像分析、电动力学、正电子发射型断层扫描术、单光子发射型断层扫描术、流体力学等问题 [1] [2] [3] [4] [5] 的处理中,高振荡积分的高效计算往往是其研究的核心。一般地,振荡积分可以写成:

0 b f ( t ) g ( t ) d t

这里 f ( ) 表示高度振荡函数, g ( ) 是缓慢变化的。

高振荡积分中的高振荡性质是导致数值计算极其困难的元凶,因为这使得传统的数值积分方法必须采用非常小的离散步长或是非常多的基底函数才能获得有意义的数值结果,然而这会导致其计算成本高昂到现有计算机系统无法承受的地步。鉴于此,高振荡积分的数值计算至今被认为是一项极具挑战性的数值难题。

过去几十年,高振荡积分计算的研究迅速发展 [6]。基于Filon的思想 [7],Iserles和Nørsett [8] 通过其Hermite插值近似缓慢变化的函数来发展Filon类型的方法。理论和数值结果表明,该方法相对于频率具有高阶收敛速度。为了得到稳定和快速的算法,Domínguez [9] 和Xiang等人 [10] 分别提出了Clenshaw-Curtis-Filon型方法,Majidian等人分析了其稳定性 [11],目前得到了广泛的应用。

虽然Filon的方法衍生出了许多有效的算法,但大多数算法都使矩积分的计算复杂化。解决这一问题的另一种方法是在复平面上改变积分路径,并推导出数值最速下降法 [12]。高振荡积分的另一个重要的求积规则是Levin法 [13],它不计算复杂矩,从而得到广泛的应用。该方法的核心是将积分问题转化为常微分方程,并通过配置法求解该方程。在 [14] 中,Li提出了一种病态系统的SVD求解器。在 [15] 中,Suliman对带有紧支撑径向基函数的高振荡积分构造了一种条件良好且有效的Levin方法。

计算高振荡积分还有许多其他重要的方法,如同伦摄动法 [16]、广义求积规则 [17]、外推法 [18],然而高振荡积分 [19] [20] [21] 的卷积积分公式被人关注较少,特别是它在高振荡情况下的渐近分析。事实上,卷积积分公式能够有效地计算卷积积分并且在解决振荡和进化问题中发挥着重要作用 [22]。因此,研究高振荡卷积积分的高效卷积求积公式具有重要意义。本文针对高振荡Bessel变换,提出改进的卷积求积公式,新公式在计算高振荡积分时不受振荡频率的影响。

2. 新型Hermite卷积积分公式的构造

考虑含高振荡Bessel核的卷积积分

0 x J m ( ω t ) g ( x t ) d t , w 1 (1)

这里 J m ( ) 指m阶第一类Bessel函数,m为非负整数,并且 J m ( ω s ) 的Laplace变换可以表示为

F J m ( s ) = ( s 2 + ω 2 s ) m ω m s 2 + ω 2

则可以得到

J m ( ω t ) = 1 2 π i Γ F J m ( λ ) e λ t d t (2)

这里的 Γ F J m ( s ) 的解析区域内由 e i ( π φ ) e i ( π φ ) 的一条曲线,将式(2)代入式(1)可得:

0 x J m ( ω t ) g ( x t ) d t = 1 2 π i Γ F J m ( λ ) 0 x e λ t g ( x t ) d t d λ (3)

注意到函数 0 x e λ t g ( x t ) d t 是初值问题

{ d y d x = λ y ( x ) + g ( x ) d 2 y d x 2 = λ d y d x + d g d x y ( 0 ) = 0 y ( 0 ) = g ( 0 ) (4)

的解。

设区间I由网格离散为 I h = { t n : = n h , n = 0 , , N , h 0 , N h = T } ,令 H ( x ) y ( x ) 在区间 [ t n , t n + 1 ] 上的两点三次Hermite插值多项式,即

H ( x ) = 2 x 3 t n + t n + 1 ( t n + 1 t n ) 3 ( x t n + 1 ) 2 y n + 2 x t n + 3 t n + 1 ( t n + 1 t n ) 3 ( x t n ) 2 y n + 1 + ( x t n + 1 ) 2 ( x t n ) ( t n + 1 t n ) 2 y n + ( x t n + 1 ) ( x t n ) 2 ( t n + 1 t n ) 2 y n + 1 , (5)

利用Hermite配置法求解上述初值问题可以近似地得到积分 0 x e λ t g ( x t ) d t 的值,即

( 0 1 β n β n + 1 ) ( y n + 1 y n + 1 ) + ( 0 1 α n α n + 1 ) ( y n y n ) = λ I ( y n + 1 y n + 1 ) + I ( g n + 1 g n + 1 ) (6)

这里 β n = 6 h 2 , β n + 1 = 4 h 1 , α n = 6 h 2 , α n + 1 = 2 h 1 。在式(6)中两边同时乘以 ζ n ( n 0 ) 并做和可得

( 0 1 β n β n + 1 ) Y ( ζ ) + ( 0 1 α n α n + 1 ) ζ Y ( ζ ) = λ I Y ( ζ ) + I G ( ζ ) (7)

其中 Y ( ζ ) = n = 0 ( y n y n ) ζ n G ( ζ ) = n = 0 ( g n g n ) ζ n ,令

δ ( ζ ) = h ( ( 0 1 β n β n + 1 ) + ( 0 1 α n α n + 1 ) ζ ) (8)

可以得到

Y ( ζ ) = ( δ ( ζ ) h λ ) 1 G ( ζ ) (9)

从而由Cauchy积分公式可得

1 2 π i Γ F J m ( λ ) y ( ζ ) d λ = 1 2 π i Γ F J m ( λ ) ( δ ( ζ ) h λ ) g ( ζ ) d λ = F J m ( δ ( ζ ) h ) g ( ζ ) , (10)

这里 ζ n 对应得系数矩阵的第一个元素即为积分(1)在 x = n h 处的一个近似。设函数 F J m ( δ ( ζ ) h ) ζ = 0 处的一个Taylor展开为

F J m ( δ ( ζ ) h ) = j = 0 W j ( h ) ζ j (11)

则计算积分(1)的新型Hermite卷积积分公式定义为

Q h H C Q ( x ) = ( 1 , 0 ) T 0 j h x W j ( h ) G ( x j h ) (12)

由Taylor展开系数的定义,可得

W j ( h ) = 1 2 π i | z | = ρ F J m ( δ ( ζ ) h ) z n 1 d z (13)

其中 ρ 足够小使得圆盘 | z | ρ 落在函数 F J m ( δ ( ζ ) h ) 的解析区域内。在上式中引入变量替换 z = ρ e i θ ,我们可以得到

W j ( h ) = ρ n 2 π 0 2 π F J m ( δ ( ρ e i θ ) h ) e i n θ d θ (14)

利用梯形公式去离散这一积分就有

0 2 π F J m ( δ ( ρ e i θ ) h ) e i n θ d θ = l = 0 L 1 l 2 π L ( l + 1 ) 2 π L F J m ( δ ( ρ e i θ ) h ) e i n θ d θ l = 0 L 1 π L ( F J m ( δ ( ρ e i l 2 π L ) h ) e i n l 2 π L + F J m ( δ ( ρ e i ( l + 1 ) 2 π L ) h ) e i n ( l + 1 ) 2 π L ) = 2 π L l = 0 L 1 F J m ( δ ( ρ e i l 2 π L ) h ) e i n l 2 π L , (15)

这里为了保证 ε 的计算精度,我们选取 L = 10 N ρ n = ε ,而积分权(13)可以通过快速Fourier变换计算,计算量为 O ( N log N )

3. 依振荡频率的收敛性分析

由于卷积积分公式的离散形式与振荡频率 ω 无关,可以预测这种方法在计算高振荡积分的时候不会随着频率的增大而损失计算精度。为了进一步研究这一类公式在计算高振荡积分时的可靠性,我们首先介绍一个引理。

引理1. [23] 假设函数 g ( ) 在区间I上连续,这里I是正半轴上的有限闭区间或者 [ a , ) a 0 且积分 S | g ( t ) | d t , S | g ( t ) | d t 存在。那么对于任意的正数 v 0 ,当 ω 趋向于无穷大时,存在与 ω 无关的常数C使得

| S g ( t ) J m ( ω t ) d t | C ω 1

进一步,如果积分区间S不包含零点,那么有

| S g ( t ) J m ( ω t ) d t | C ω 3 / 2

下面考虑Hermite卷积积分公式计算高振荡积分时的一些性质。

定理1. 在积分(16)中,若函数 g ( ) 在正半轴上二次连续可微,且积分 0 b | g ( t ) | d t 0 b | g ( t ) | d t 0 b | g ( t ) | d t 存在。那么当 ω 趋向于无穷大时,存在常数 C > 0 ,使得Hermite卷积积分公式的误差满足

| Q h H C Q ( b ) I [ g ] | C ω 3 / 2

证明:令 δ j δ ( ζ ) = j = 0 δ j ζ j ζ j 的系数,并且定义

s h G ( x ) : = 1 h 0 j h x δ j G ( x j h ) , 0 < x b (16)

这里 G ( x ) = ( g ( x ) g ( x ) )

b = ( 1 , 0 ) T ,由卷积积分公式的定义,其可以表示成Dunford-Taylor积分的形式

F ( s h ) g ( x ) = 1 2 π i Γ F J m ( λ ) b T ( s h λ ) 1 G ( x ) d λ (17)

又因为

F J m ( s ) = 0 e t s J m ( ω t ) d t (18)

将上式代入等式(17)得

F ( s h ) g ( x ) = 1 2 π i Γ 0 e t λ J m ( ω t ) b T ( s h λ ) 1 G ( x ) d t d λ = 0 J m ( ω t ) b T e t s h G ( x ) d t , (19)

I [ g ] 表示成如下形式

F ( s ) g ( x ) = 0 x J m ( ω t ) g ( x t ) d t = 0 x J m ( ω t ) b T G ( x t ) d t (20)

因此

F ( s h ) g ( x ) F ( s ) g ( x ) = 0 x J m ( ω t ) ( b T e t s h G ( x ) b T G ( x t ) ) d t + x J m ( ω t ) b T e t s h G ( x ) d t . (21)

定义 ϕ 1 ( t ) : = b T e t s h G ( x ) b T G ( x t ) ϕ 2 ( t ) : = b T e t s h G ( x ) ,由Bessel函数的求导性质可得

F ( s h ) g ( x ) F ( s ) g ( x ) = 0 x J m ( ω t ) ϕ 1 ( t ) d t + x J m ( ω t ) ϕ 2 ( t ) d t = ϕ 1 ( 0 ) 0 x J m ( ω t ) d t + 1 ω ( ϕ 1 ( t ) ϕ 1 ( 0 ) ) J m + 1 ( ω t ) | t = 0 x 1 ω 0 x ( ϕ 1 ( t ) m + 1 t ( ϕ 1 ( t ) ϕ 1 ( 0 ) ) ) J m + 1 ( ω t ) d t + 1 ω ϕ 2 ( t ) J m + 1 ( ω t ) | t = x 1 ω x ( ϕ 2 ( t ) m + 1 t ϕ 2 ( t ) ) J m + 1 ( ω t ) d t , (22)

ϕ 1 ( 0 ) = 0 和引理1可得,存在常数 C > 0 使得

| F ( s h ) f ( x ) F ( s ) f ( x ) | C ω 3 / 2

证毕。

4. 数值实验

在本节中,我们通过几个数值算例来诠释上文提出的Hermite卷积积分公式在计算高振荡积分时的高效性。

例1. 我们考虑以下两个高振荡积分

I 1 = 0 2 J 0 ( ω ( 2 t ) ) 1 1 + 25 t 2 d t I 2 = 0 1 J 1 ( ω ( 1 t ) ) cos ( t ) e t d t (23)

表1中,我们可以看到这种方法计算上述积分得到的关于不同频率 ω 的绝对误差。在图1图2中,我们展示了Hermite卷积求积公式计算此类高振荡积分的收敛速度。数值结果表明,Hermite卷积积分公式在计算此类高振荡积分方面具有较好效果,其关于频率 ω 渐近阶能够达到 O ( ω 3 / 2 ) 。可以看出,这些计算结果与之前的定理中给出的渐近阶的估计是吻合的。

Table 1. The absolute errors of the numerical solution obtained by using the Hermite Convolution Quadtrature to solve integral (23)

表1. 利用Hermite卷积积分公式求(23)所得数值解的绝对误差

5. 在高振荡Volterra积分方程数值解中的应用

5.1. Hermite卷积积分公式求解高振荡Volterra积分方程

在物理学、工程等领域的诸多问题最终都可以转化为第一类卷积型Volterra积分方程的求解。例如,对于二维声学散射问题中的单层势能积分方程

1 2 π R 2 u ( x , t | x x | ) | x x | d x = a ( x , t ) , ( x , t ) R 2 × ( 0 , T )

Figure 1. Absolute errors of I 1 (left) and its transformation law with respect to ω (right)

图1. I 1 的绝对误差(左)与绝对误差关于 ω 的变换规律(右)

Figure 2. Absolute errors of I 2 (left) and its transformation law with respect to ω (right)

图2. I 2 的绝对误差(左)与绝对误差关于 ω 的变换规律(右)

t 0 ,满足 u 0 a 0 时,上式可以利用连续Fourier变换转化为带Bessel核的第一类卷积型Volterra积分方程:

0 x J 0 ( ω ( x t ) ) u ( t ) d t = g ( x ) , x [ 0 , T ] (24)

这里 g ( ) 足够光滑且 g ( 0 ) = 0

此时 I h = { t n : = n h , n = 0 , , N , h 0 , N h = T } [ 0 , T ] 上的等距网格,且 h = T N 。对此方程应用Hermite卷积积分公式,可以得到

( 1 , 0 ) T j = 0 k W k j ( u j u j ) = g k (25)

对(24)求导可得

0 x ω J 0 ( ω ( x t ) ) u ( t ) d t + u ( x ) = g ( x ) (26)

对此方程应用Hermite卷积积分公式,可以得到

( 1 , 0 ) T j = 0 k W k j ( u j u j ) + u k = g k (27)

这里的W和 W 是Hermite卷积积分对应的权矩阵, u j u j 分别表示 u ( x ) u ( x ) t j 处的数值解, g j g j 分别表示 g ( x ) g ( x ) t j 处的值。联立(25)和(27)得到一个线性系统,对此线性系统求解即可得到均匀网格 I h 处的数值解。

5.2. 数值实验

卷积求积公式是求解具有卷积核的Volterra方程的重要工具。许多数值实验表明,它们对于求解一些高振荡的Volterra积分方程是有效的。

例2. 考虑以下高振荡Volterra积分方程

0 x J 0 ( ω ( x t ) ) u ( t ) d t = x e x (28)

在这个例子中,我们考虑了Hermite卷积积分公式在求解高振荡Volterra积分方程中的应用。首先,通过让 T = 2 N = 4 ,我们在表2图3中给出了此方法关于不同频率 ω 的绝对误差。这些表和图中列出的绝对误差表明,此方法具有振荡频率越高,近似值就越好的性质。

Table 2. The absolute errors of the numerical solution obtained by using the Hermite convolution quadtrature to solve Equation (28)

表2. 利用Hermite卷积积分公式求解方程(28)所得数值解的绝对误差

Figure 3. The absolute errors of the numerical solution obtained by using the Hermite convolution quadtrature to solve equation (28) ( ω = 100 , 1000 )

图3. Hermite卷积积分公式求解(28)所得数值解的绝对误差( ω = 100 , 1000 )

6. 总结

本文提出了高振荡Bessel变换的新型Hermite卷积积分公式,研究了此公式在处理高振荡积分数值计算问题时的一些性质。第2节的理论结果揭示了该公式相对于频率 ω 的收敛速度,当 ω 趋于无穷大时,Hermite卷积积分公式以频率 ω 的负幂收敛,即 O ( ω 3 / 2 ) 。在第3节中,我们通过两个数值算例验证了第2节中理论的准确性。在第4节我们应用Hermite卷积积分公式求解带高振荡Bessel核的Volterra卷积积分方程,通过数值实验表明,此类方法是有效的。此外,本文仅仅为高振荡积分数值计算问题的新型Hermite求积公式收敛理论打开了一扇窗,未来需要对此方法进行改进,从而得到更高的收敛阶。

参考文献

NOTES

*通讯作者。

参考文献

[1] Asheim, A. and Huybrechs, D. (2009) Local Solutions to High Frequency 2D Scattering Problems. Journal of Computational Physics, 229, 5357-5372.
https://doi.org/10.1016/j.jcp.2010.03.034
[2] Cakoni, F. and Colton, D. (2006) Qualitative Methods in Inverse Scattering Theory. Springer, Berlin.
https://doi.org/10.1515/jiip.2007.027
[3] Colton, D. and Kress, R. (1983) Integral Equation Methods in Scattering Theory. Wiley, Hoboken.
[4] Langdon, S. and Chandler-Wilde, S.N. (2006) A Wavenumber Independent Boundary Element Method for an Acoustic Scattering Problem. SIAM Journal on Numerical Analysis, 43, 2450-2477.
https://doi.org/10.1137/S0036142903431936
[5] Nédélec, J. (2001) Acoustic and Electromagnetic Equations. Springer, Berlin, 144.
https://doi.org/10.1007/978-1-4757-4393-7
[6] 向淑晃. 一些高振荡积分、高振荡积分方程的高性能计算[J]. 中国科学: 数学, 2012, 42(7): 651-670.
[7] Filon, L. (1930) On a Quadrature Formula for Trigonometric Integrals. Proceedings of the Royal Society of Edinburgh, 49, 38-47.
https://doi.org/10.1017/S0370164600026262
[8] Nørsett, S. and Iserles, A. (2005) Efficient Quadrature of Highly Oscillatory Integrals Using Derivatives. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 461, 1383-1399.
https://doi.org/10.1098/rspa.2004.1401
[9] Domínguez, V., Graham, I. and Smyshlyaev, V. (2011) Stability and Error Estimates for Filon-Clenshaw-Curtis Rules for Highly-Oscillatory Integrals. IMA Journal of Numerical Analysis, 31, 1253-1280.
https://doi.org/10.1093/imanum/drq036
[10] Xiang, S., Cho, Y., Wang, H., et al. (2011) Clenshaw-Curtis-Filon-Type Methods for Highly Oscillatory Bessel Transforms and Applications. IMA Journal of Numerical Analysis, 31, 1281-1314.
https://doi.org/10.1093/imanum/drq035
[11] Majidian, H., Firouzi, M. and Alipanah, A. (2022) On the Stability of Filon-Clenshaw-Curtis Rules. Bulletin of the Iranian Mathematical Society, 1-22.
https://doi.org/10.1007/s41980-022-00681-4
[12] Wu, Q. and Sun, M. (2021) Numerical Steepest Descent Method for Hankel Type of Hypersingular Oscillatory Integrals in Electromagnetic Scattering Problems. Advances in Mathematical Physics, 2021, Article ID: 8021050.
https://doi.org/10.1155/2021/8021050
[13] Levin, D. (1982) Procedures for Computing One- and Two-Dimensional Integrals of Functions with Rapid Irregular Oscillations. Mathematics of Computation, 38, 531-538.
https://doi.org/10.1090/S0025-5718-1982-0645668-7
[14] Li, J., Wang, X., Wang, T., et al. (2010) An Improved Levin Quadrature Method for Highly Oscillatory Integrals. Applied Numerical Mathematics, 60, 833-842.
https://doi.org/10.1016/j.apnum.2010.04.009
[15] Khan, S., Zaman, S., Arshad, M., et al. (2021) A Well-Conditioned and Efficient Levin Method for Highly Oscillatory Integrals with Compactly Supported Radial Basis Functions. Engineering Analysis with Boundary Elements, 131, 51-63.
https://doi.org/10.1016/j.enganabound.2021.06.012
[16] He, J. and El-Dib, Y. (2021) Homotopy Perturbation Method with Three Expansions for Helmholtz-Fangzhu Oscillator. International Journal of Modern Physics B, 35, Article ID: 2150244.
https://doi.org/10.1142/S0217979221502441
[17] Kang, H., Xiang, C., Xu, Z., et al. (2021) Efficient Quadrature Rules for the Singularly Oscillatory Bessel Transforms and Their Error Analysis. Numerical Algorithms, 88, 1493-1521.
https://doi.org/10.1007/s11075-021-01083-z
[18] Ke, X., Ying, S. and Wang, Y. (2021) One-Step Extrapolation Method Based on the Multiple-Angle Formula. Chinese Journal of Geophysics, 64, 2480-2493.
[19] Lubich, C. (1988) Convolution Quadrature and Discretized Operational Calculus. I. Numerische Mathematik, 52, 129-145.
https://doi.org/10.1007/BF01398686
[20] Lubich, C. (1988) Convolution Quadrature and Discretized Operational Calculus. II. Numerische Mathematik, 52, 413-425.
https://doi.org/10.1007/BF01462237
[21] Lubich, C. and Ostermann, A. (1993) Runge-Kutta Methods for Parabolic Equations and Convolution Quadrature. Mathematics of Computation, 60, 105-131.
https://doi.org/10.1090/S0025-5718-1993-1153166-7
[22] Xiang, S. and Brunner, H. (2013) Efficient Methods for Volterra Integral Equations with Highly Oscillatory Bessel Kernels. BIT Numerical Mathematics, 53, 241-263.
https://doi.org/10.1007/s10543-012-0399-8
[23] Ma, J. and Liu, H. (2018) On the Convolution Quadrature Rule for Integral Transforms with Oscillatory Bessel Kernels. Symmetry, 10, 239.
https://doi.org/10.3390/sym10070239