Next Article in Journal
(Re-)Reading Sklar (1959)—A Personal View on Sklar’s Theorem
Previous Article in Journal
Mitigation of the Collision Risk of a Virtual Impactor Based on the 2011 AG5 Asteroid Using a Kinetic Impactor
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Efficient Quadrature Rule for Highly Oscillatory Integrals with Airy Function

1
School of Mathematics, Nanjing Audit University, Nanjing 211815, China
2
College of Mathematics and Information Science, Zhengzhou University of Light Industry, Zhengzhou 450002, China
3
School of Mathematics and Big Data, Guizhou Education University, Guiyang 550018, China
*
Author to whom correspondence should be addressed.
Mathematics 2024, 12(3), 377; https://doi.org/10.3390/math12030377
Submission received: 1 January 2024 / Revised: 20 January 2024 / Accepted: 23 January 2024 / Published: 24 January 2024
(This article belongs to the Section Computational and Applied Mathematics)

Abstract

:
In this work, our primary focus is on the numerical computation of highly oscillatory integrals involving the Airy function. Specifically, we address integrals of the form 0 b x α f ( x ) Ai ( ω x ) d x over a finite or semi-infinite interval, where the integrand exhibits rapid oscillations when ω 1 . The inherent high oscillation and algebraic singularity of the integrand make traditional quadrature rules impractical. In view of this, we strategically partition the interval into two segments: [ 0 , 1 ] and [ 1 , b ] . For integrals over the interval [ 0 , 1 ] , we introduce a Filon-type method based on a two-point Taylor expansion. In contrast, for integrals over [ 1 , b ] , we transform the Airy function into the first kind of Bessel function. By applying Cauchy’s integration theorem, the integral is then reformulated into several non-oscillatory and exponentially decaying integrals over [ 0 , + ) , which can be accurately approximated by the generalized Gaussian quadrature rule. The proposed methods are accompanied by rigorous error analyses to establish their reliability. Finally, we present a series of numerical examples that not only validate the theoretical results but also showcase the accuracy and efficacy of the proposed method.

1. Introduction

The numerical computation of highly oscillatory integrals arises in various scientific disciplines, ranging from physics to engineering and beyond [1,2,3,4,5,6]. Oscillatory integrals are characterized by rapidly oscillating integrands, and their accurate computation poses a formidable challenge due to the potential for high-frequency oscillations and significant cancellations. This challenge is particularly pronounced in integrals involving Airy functions, which are essential mathematical functions frequently encountered in the study of wave phenomena, quantum mechanics, and the analysis of asymptotic behavior [1,7,8,9,10].
The Airy function Ai ( ω x ) , one of the solutions to the differential equation y + ω 3 x y = 0 , plays a significant role in applied mathematics and has been extensively studied for its oscillatory behavior when ω 1 and x 0 (see, for instance, [11,12,13,14]). Figure 1 illustrates the rapid oscillation of Ai ( ω x ) , with the oscillations becoming more pronounced as ω increases. As x + , the following asymptotic formula is valid [13] [Equation (9.7.9)]:
Ai ( x ) x 1 / 4 π cos ( ζ π / 4 ) k = 0 ( 1 ) k u 2 k ζ 2 k + sin ( ζ π / 4 ) k = 0 ( 1 ) k u 2 k + 1 ζ 2 k + 1 ,
where ζ = 2 3 x 3 / 2 , u 0 = 1 and
u k = ( 2 k + 1 ) ( 2 k + 3 ) ( 6 k 1 ) ( 216 ) k k ! , k = 1 , 2 , .
The numerical computation of integrals involving Airy functions is a difficult task, primarily due to their potential rapid oscillations over the interval. Standard numerical methods, such as Newton–Cotes and Gaussian quadrature rules, often struggle to handle the rapid oscillations in such integrals, necessitating the development of advanced techniques to ensure computational reliability.
This study aims to investigate numerical techniques for efficiently and accurately computing highly oscillatory integrals of the form
I [ f ] = 0 b x α f ( x ) Ai ( ω x ) d x , 0 < b + ,
where α > 1 , and f is assumed to be smooth on [ 0 , b ] and analytic in the complex plane { z | 1 ( z ) b } . In particular, for the case b = + , it is established from (1) that ensuring the existence of (2) requires the presence of constants σ < α 3 / 4 and C such that f ( x ) C x σ as x + . In the past few years, many works have been devoted to the numerical treatment of the integrals in (2). The Levin-type method, developed by Levin [15] and later improved by Olver [12], transforms the solving of highly oscillatory integrals into an ordinary differential equation. Xiang and Wang (2010) introduced the Filon-type method for a more generalized integral 0 1 f ( x ) Ai ( ω g ( x ) ) d x with g ( x ) > 0 , employing a diffeomorphism transformation on the phase function g ( x ) [16]. Xu and Xiang (2014) developed the Clenshaw-Curtis-Filon method specifically for the computation of (2) with b = 1 . Asymptotic expansions of 0 b f ( x ) Ai ( ω x ) d x with b < + have been explored by Olver (2007) and Kang (2018, 2020) [11,17,18]. All these methods share a remarkable advantage in that their accuracy improves rapidly as the frequency ω increases. However, each approach has limitations with effectiveness depending on the specific characteristics of the integral under consideration.
However, these existing numerical methods face limitations when directly applied to the integral (2), primarily due to the singularity at x = 0 and the unbounded nature of the interval when b = + . In order to avoid this defect, we divide the interval into two distinct segments and rewrite (2) as
I [ f ] = 0 a x α f ( x ) Ai ( ω x ) d x + a b x α f ( x ) Ai ( ω x ) d x = : I 1 [ f ] + I 2 [ f ] .
Without loss of generality, we assume b > 1 and set a = 1 in this paper. For b 1 , the interval [ 0 , b ] of (2) can be transformed to [ 0 , 1 ] under a linear transformation. Then, the Filon-type method introduced in Section 2 will be applied to approximate (2). The algebraic singularity at x = 0 is exclusively present in I 1 [ f ] and will be treated by the Filon-type method. For the evaluation of I 2 [ f ] , a complex integration method is introduced by applying Cauchy’s integration theorem.
This work is organized as follows: In Section 2, we introduce the Filon-type method for computing I 1 [ f ] based on a two-point Taylor formula. Section 3 presents the complex integration method for approximating I 2 [ f ] , where the highly oscillatory integrals are transformed into several non-oscillatory and exponentially decaying integrands on [ 0 , ) . The application of a Gauss-Laguerre quadrature rule and its error estimates are discussed. Section 4 provides numerical examples to demonstrate the effectiveness of the proposed methods. Finally, we conclude this paper in Section 5.

2. The Filon-Type Method for I 1 [ f ]

The Filon method [19], named after the French mathematician L.N.G. Filon, stands out as an effective numerical technique designed specifically for a distinct class of highly oscillatory integrals a b f ( x ) e i ω x d x . In this method, the non-oscillatory component f ( x ) is approximated by a piecewise polynomial. This fundamental concept underwent further refinement by Iserles and Nørsett [5], who made a significant advancement by recognizing that the asymptotic expansion of a b f ( x ) e i ω x d x solely depends on the values and derivatives of f ( x ) evaluated at the endpoints. Recently, this strategic approach has been extended to a broader spectrum of highly oscillatory integrals, including those with Bessel-type kernels [16,20].
In this section, the Filon-type method will be extended to the highly oscillatory integrals I 1 [ f ] with the Airy kernel. Let m 1 be a positive integer. In order to obtain a higher asymptotic order, we seek a Hermite’s interpolation polynomial P 2 m 1 ( x ) of f ( x ) that satisfies
P 2 m 1 ( 0 ) = f ( 0 ) , P 2 m 1 ( 0 ) = f ( 0 ) , , P 2 m 1 ( m 1 ) ( 0 ) = f ( m 1 ) ( 0 ) ; P 2 m 1 ( 1 ) = f ( 1 ) , P 2 m 1 ( 1 ) = f ( 1 ) , , P 2 m 1 ( m 1 ) ( 1 ) = f ( m 1 ) ( 1 ) .
Suppose that f ( x ) is ( m 1 ) -times differentiable on [ 0 , 1 ] , the unique polynomial of degree 2 m 1 that satisfies the condition (4) is the two-point Taylor formula, represented as
P 2 m 1 ( x ) = k = 0 m 1 c k x k + 1 ( x 1 ) k + k = 0 m 1 d k ( x 1 ) k + 1 x k ,
where the coefficients are given by [21]
c k = f ( 1 ) , k = 0 , j = 0 k ( k + j 1 ) ! j ! ( k j ) ! ( 1 ) j k f ( k j ) ( 1 ) ( 1 ) k j f ( k j ) ( 0 ) k ! , k 1 ,
and
d k = f ( 0 ) , k = 0 , j = 0 k ( k + j 1 ) ! j ! ( k j ) ! ( 1 ) j j f ( k j ) ( 1 ) ( 1 ) k k f ( k j ) ( 0 ) k ! , k 1 .
As an approximation to f ( x ) , the two-point Taylor polynomial converges to f ( x ) at the rate o ( x m ) as x 0 , and at the rate o ( ( x 1 ) m ) as x 1 . Here, o ( x m ) is defined in the sense that lim x 0 f ( x ) P 2 m 1 ( x ) x m = 0 .
By substituting P 2 m 1 ( x ) into I 1 [ f ] , we obtain an approximation of I 1 [ f ] given by
I 1 [ f ] Q m [ I 1 ] : = k = 0 m 1 c k M k ( k + 1 + α , ω ) + k = 0 m 1 d k M k + 1 ( k + α , ω ) ,
where M k ( ρ , ω ) are the moments defined by
M k ( ρ , ω ) = 0 1 x ρ ( x 1 ) k Ai ( ω x ) d x .
An explicit formula for the moments is given by Xu [14] using the Meijer G-function.
Lemma 1 
([14]). When ρ > 1 , the moments defined in (9) satisfy
M k ( ρ , ω ) = ( 1 ) k π k ! 3 k + 2 [ 3 3 G 4 , 6 1 , 3 ρ 3 , 1 ρ 3 , 2 ρ 3 , 1 2 0 , 1 3 , 1 2 , ρ k 1 3 , ρ k 3 , ρ k + 1 3 | ω 3 9 + ω 3 3 G 4 , 6 1 , 3 ρ 1 3 , ρ 3 , 1 ρ 3 , 1 2 0 , 1 3 , 1 2 , ρ k 2 3 , ρ k 1 3 , ρ k 3 | ω 3 9 ] ,
where the Meijer G-function involved is defined as (cf. [13])
G p , q m , n a 1 , , a n , a n + 1 , , a p b 1 , , b m , b m + 1 , , b q | z = 1 2 π i L k = 1 m Γ ( b k s ) j = 1 n Γ ( 1 a j + s ) k = m + 1 q Γ ( 1 b k + s ) j = n + 1 p Γ ( a j s ) z s d s .
Remark 1. 
The contour integration technique employed in (11) also offers an efficient and accurate implementation option for Meijer G-functions in software applications, including popular platforms like Matlab [22] and Mathematica.
Now, we turn our attention to the asymptotic properties of the errors in the Filon-type method. For simplicity, let E m ( x ) : = f ( x ) P 2 m 1 ( x ) . It is evident that E m ( k ) ( 0 ) = 0 and E m ( k ) ( 1 ) = 0 for any 0 k m 1 , defining
σ 0 ( x ) = x α E m ( x ) , σ k + 1 = σ k ( x ) x , k 0 ,
the asymptotic property of I 1 [ f ] Q m [ I 1 ] is provided as follows.
Theorem 1. 
Let α > 1 be real but not an integer. As ω , the Filon-type quadrature defined in (8) satisfies
I 1 [ f ] Q m [ I 1 ] = O ( ω m α 1 ) .
Proof. 
Recalling Airy’s equation,
d 2 d x 2 Ai ( x ) x Ai ( x ) = 0 ,
we have
I 1 [ f ] Q m [ I 1 ] = 0 1 σ 0 ( x ) Ai ( ω x ) d x = 1 ω 3 0 1 σ 0 ( x ) x ( Ai ( ω x ) ) d x = 1 ω 3 σ 0 ( x ) x ( Ai ( ω x ) ) | x = 0 1 0 1 σ 0 ( x ) x ( Ai ( ω x ) ) d x = 1 ω 3 σ 0 ( x ) x ( Ai ( ω x ) ) σ 0 ( x ) x Ai ( ω x ) x = 0 1 + 1 ω 3 0 1 σ 1 ( x ) Ai ( ω x ) d x .
Let κ be the largest integer such that κ < min { m + α 3 , m 2 } . After an implementation of κ -times of integration (by parts), we obtain
I 1 [ f ] Q m [ I 1 ] = j = 0 κ 1 1 ω 3 ( j + 1 ) σ j x ( Ai ( ω x ) ) σ j x Ai ( ω x ) x = 0 1 + 1 ω 3 κ 0 1 σ κ ( x ) Ai ( ω x ) d x .
From E m ( j ) ( 0 ) = E m ( j ) ( 1 ) = 0 for j = 0 , 1 , , m 1 , we know that there exist functions h j ( x ) such that
σ j = x m + α 3 j ( x 1 ) m 2 j h j ( x ) , j = 0 , 1 , , κ
where h j ( x ) is smooth with h j ( 0 ) 0 and h j ( 1 ) 0 . Thus, the summation in the right-hand side of (16) vanishes. By substituting (17) into (16), we get
I 1 [ f ] Q m [ I 1 ] = 1 ω 3 κ 0 1 x m + α 3 κ ( x 1 ) m 2 κ h κ ( x ) Ai ( ω x ) d x .
In the following, we consider the asymptotic property of the integral in (18). By denoting
Ψ ω ( x ) = 0 x t m + α 3 κ Ai ( ω t ) d t ,
we know from [14] (Theorem 3.1) that Ψ ω ( x ) = O ( ω 3 κ m α 1 ) holds uniformly for all x [ 0 , 1 ] . As a result, it follows that
I 1 [ f ] Q m [ I 1 ] = 1 ω 3 κ 0 1 ( x 1 ) m 2 κ h κ ( x ) d Ψ ω ( x ) = 1 ω 3 κ ( x 1 ) m 2 κ h κ ( x ) Ψ ω ( x ) | x = 0 1 1 ω 3 κ 0 1 Ψ ω ( x ) d { ( x 1 ) m 2 κ h κ ( x ) } = O ( ω m α 1 ) .
This completes the proof. □
In Figure 2, we demonstrate the asymptotic order of the error of the Filon-type method for I 1 [ f ] , where f ( x ) = e 1 / x 2 C [ 0 , 1 ] is an infinitely differentiable function on [ 0 , 1 ] . The absolute errors of Q m [ I 1 ] for various values of ω are displayed. It is evident that the accuracy of Q m [ I 1 ] improves as ω increases and the decay rate is O ( ω m α 1 ) , which confirms the conclusion in Theorem 1.

3. Complex Integration Method for I 2 [ f ]

This section is dedicated to the numerical computation of I 2 [ f ] defined in (3), where the function f is assumed to be analytic in the complex plane { z | 1 ( z ) b } . The complex integration method will be applied to transform I 2 [ f ] into some non-oscillatory and exponentially decaying integrals.
By utilizing the identity derived from [13] [Equation (9.6.6)],
Ai ( x ) = x 3 J 1 / 3 2 3 x 3 / 2 + J 1 / 3 2 3 x 3 / 2 ,
where J ν ( x ) denotes the first kind of Bessel function of order ν ; we substitute (21) into I 2 [ f ] , denoting r = 2 3 ω 3 / 2 to obtain
I 2 [ f ] = 1 3 3 2 r 1 / 3 1 b x α + 1 / 2 f ( x ) J 1 / 3 ( r x 3 / 2 ) + J 1 / 3 ( r x 3 / 2 ) d x = 2 9 3 2 r 1 / 3 1 b 3 / 2 x 2 α / 3 f ( x 2 / 3 ) J 1 / 3 ( r x ) + J 1 / 3 ( r x ) d x .
By recalling the integral representation of Bessel function, as given in [23] [Equation (8.411.10)], namely
J ν ( x ) = x ν 2 ν π Γ ( 1 / 2 + ν ) 1 1 ( 1 u 2 ) ν 1 / 2 e i x u d u , ( ν ) > 1 / 2 ,
where i = 1 is the imaginary unit, we apply Cauchy’s integral theorem [24] along with contributions from Chen [25] and Zaman et al. [26]. This allows us to express I 2 [ f ] as
I 2 [ f ] = 6 3 i 9 π Γ ( 5 / 6 ) 1 b 3 / 2 Φ 1 / 3 + ( x ) e i r x d x 1 b 3 / 2 Φ 1 / 3 ( x ) e i r x d x + 2 3 3 i 9 π Γ ( 1 / 6 ) r 2 / 3 1 b 3 / 2 Φ 1 / 3 + ( x ) e i r x d x 1 b 3 / 2 Φ 1 / 3 ( x ) e i r x d x ,
where Φ ν ± ( x ) is defined as
Φ ν ± ( x ) = x 2 α / 3 ν f ( x 2 / 3 ) 0 + y ν 1 / 2 ( y ± 2 i r x ) ν 1 / 2 e y d y .

3.1. Complex Integration Formulas

It is noteworthy that the integrands Φ ν ± ( x ) defined in (25) are non-oscillatory and analytic in the complex plane { z | ( z ) 1 } . This characteristic makes the evaluation of the highly oscillatory integrals in (24) a typical problem and can be efficiently solved by the complex integration method [27,28].
Theorem 2. 
Let f ( z ) be an analytic function in the complex plane ( z ) 1 , and assume there exist some constants C and ω 0 ( 0 , ω ) such that | f ( z ) | C e ω 0 ( z ) . Then, we establish
I 2 [ f ] = 6 3 9 π Γ ( 5 / 6 ) [ e i r r 0 + Φ 1 / 3 + 1 i u r e u d u + e i r r 0 + Φ 1 / 3 1 + i u r e u d u e i r b 3 / 2 r 0 + Φ 1 / 3 + b 3 / 2 i u r e u d u e i r b 3 / 2 r 0 + Φ 1 / 3 b 3 / 2 + i u r e u d u ] + 2 3 3 r 2 / 3 9 π Γ ( 1 / 6 ) [ e i r r 0 + Φ 1 / 3 + 1 i u r e u d u + e i r r 0 + Φ 1 / 3 1 + i u r e u d u e i r b 3 / 2 r 0 + Φ 1 / 3 + b 3 / 2 i u r e u d u e i r b 3 / 2 r 0 + Φ 1 / 3 b 3 / 2 + i u r e u d u ] ,
where r = 2 3 ω 3 / 2 and Φ ν ± ( · ) are defined in (25).
Proof. 
According to the assumption based on f, it is evident that Φ ν + ( z ) is analytic in the complex region ( z ) 1 . By using Cauchy’s integration theorem [24], we express the integral as
1 b 3 / 2 Φ ν + ( x ) e i r x d x = Γ 1 + Γ 2 + Γ 3 Φ ν + ( z ) e i r z d z ,
where the integration paths Γ 1 , Γ 2 , and Γ 3 are shown in the left of Figure 3. In the following, we parameterize the integrals over these paths respectively.
For the path Γ 1 , setting z = 1 i u with u [ 0 , R ] yields
Γ 1 Φ ν + ( z ) e i r z d z = i e i r 0 R Φ ν + ( 1 i u ) e r u d u .
In the same way, for Γ 3 , setting z = b 3 / 2 i u yields
Γ 3 Φ ν + ( z ) e i r z d z = i e i r b 3 / 2 0 R Φ ν + ( b 3 / 2 i u ) e r u d u .
For Γ 2 , a change in the variable z = u i R with u [ 1 , b 3 / 2 ] results in
Γ 2 Φ ν + ( z ) e i r z d z = e r R 1 b 3 / 2 Φ ν + ( u i R ) e i r u d u .
Letting R + , it is obvious that
Γ 2 Φ ν + ( z ) e i r z d z e r R 1 b 3 / 2 | Φ ν + ( u i r u ) | d u 0 .
Combined with (27), (28), (31), and (29), we derive
1 b 3 / 2 Φ ν + ( x ) e i r x d x = i e i r r 0 + Φ ν + 1 i u r e u d u + i e i r b 3 / 2 r 0 + Φ ν + b 3 / 2 i u r e u d u .
After a similar argument and employing the integration paths depicted in the right of Figure 3, this leads to
1 b 3 / 2 Φ ν ( x ) e i r x d x = i e i r r 0 + Φ ν 1 + i u r e u d u + i e i r b 3 / 2 r 0 + Φ ν b 3 / 2 + i u r e u d u .
This completes the proof by setting ν = 1 / 3 and ν = 1 / 3 , respectively. □
In particular, when b = + , the complex integration method becomes particularly efficient by altering the integration paths of (24), as illustrated in Figure 4.
Theorem 3. 
Let f ( z ) be an analytic function in the complex plane ( z ) 1 , and assume that there exist some constants C and δ > 0 such that | f ( z ) | C | z | δ for all | z | 1 . For b = + , we establish that
I 2 [ f ] = 6 3 9 π Γ ( 5 / 6 ) e i r r 0 + Φ 1 / 3 + 1 i u r e u d u + e i r r 0 + Φ 1 / 3 1 + i u r e u d u + 2 3 3 r 2 / 3 9 π Γ ( 1 / 6 ) e i r r 0 + Φ 1 / 3 + 1 i u r e u d u + e i r r 0 + Φ 1 / 3 1 + i u r e u d u ,
where r = 2 3 ω 3 / 2 .
Proof. 
Applying Cauchy’s integral theorem [24], we have
1 + Φ ν + ( x ) e i r x d x = Γ 1 + Γ 2 Φ ν + ( z ) e i r z d z , as R + ,
where Γ 1 and Γ 2 are complex paths shown in the left of Figure 4. We parameterize Γ 1 as z = 1 i u , u [ 0 , R ] , then
Γ 1 Φ ν + ( z ) e i r z d z = e i r 0 R Φ ν + ( 1 i u ) e r u d u .
For Γ 2 , setting z = 1 + R e i θ with θ [ 0 , π / 2 ] , we obatin
Γ 2 Φ ν + ( z ) e i r z d z = i e i r R 0 π / 2 e i ( r R cos θ + θ ) e r R sin θ Φ ν + ( 1 + R e i θ ) d θ .
As R + , it follows that
Γ 2 Φ ν + ( z ) e i r z d z R 0 π / 2 e r R sin θ | Φ ν + ( 1 + R e i θ ) | d θ C R 1 δ 0 π / 2 e 2 r R θ / π d θ C π R 1 δ 1 e r R 2 r R 0 ,
where we used sin θ > 2 θ / π . As a result, we obtain that
1 + Φ ν + e i r x d x = e i r 0 + Φ ν + ( 1 i u ) e r u d u .
Similarly, by using the complex paths in the right of Figure 4, we get
0 + Φ ν ( x ) e i r x d x = e i r 0 + Φ ν ( 1 + i u ) e r u d u .
Substituting these into (24) leads to the desired result (34). □

3.2. Numerical Algorithms and Error Estimates

Apparently, the integrands in (26) are smooth and accompanied by the Laguerre weight function. This characteristic makes them amenable to be accurately approximated through the Gaussian quadrature rule. Consequently, employing the n-point Gauss-Laguerre quadrature rule provides a robust numerical method for evaluating I 2 [ f ] :
Q n [ I 2 ] : = 6 3 9 π Γ ( 5 / 6 ) k = 1 n w k ( 0 ) [ e i r r Φ ^ 1 / 3 + 1 i r u k ( 0 ) + e i r r Φ ^ 1 / 3 1 + i r u k ( 0 ) e i r b 3 / 2 r Φ ^ 1 / 3 + b 3 / 2 i r u k ( 0 ) e i r b 3 / 2 r Φ ^ 1 / 3 b 3 / 2 + i r u k ( 0 ) ] + 2 3 3 r 2 / 3 9 π Γ ( 1 / 6 ) k = 1 n w k ( 0 ) [ e i r r Φ ^ 1 / 3 + 1 i r u k ( 0 ) + e i r r Φ ^ 1 / 3 1 + i r u k ( 0 ) e i r b 3 / 2 r Φ ^ 1 / 3 + b 3 / 2 i r u k ( 0 ) e i r b 3 / 2 r Φ ^ 1 / 3 b 3 / 2 + i r u k ( 0 ) ] ,
where Φ ^ ν ± ( x ) is defined as a numerical approximation of Φ ν ± ( x ) by the Gaussian quadrature rule such that
Φ ν ± ( x ) Φ ^ ν ± ( x ) : = x 2 α / 3 ν f ( x 2 / 3 ) k = 1 n w k ( ν 1 / 2 ) u k ( ν 1 / 2 ) ± 2 i r x ν 1 / 2 .
Here, by using { u k ( γ ) , w k ( γ ) } , we denote the quadrature nodes and weights of the generalized Laguerre quadrature rule corresponding to the weight function u γ e u with γ > 1 .
Remark 2. 
Algorithms for computing Gaussian quadrature nodes and weights have been extensively studied over the past several decades. The classic approach is the Golub-Welsch algorithm, which solves an eigenvalue problem and requires O ( n 2 ) operations. However, when n is large, several fast algorithms have been presented by Glaser, Liu, and Rokhlin [29], as well as Hale and Trefethen [30]. These algorithms require only O ( n ) operations by leveraging the asymptotic properties of orthogonal polynomials.
Theorem 4. 
Suppose that f ( x ) satisfies the condition in Theorem 2. Then, as ω + , the absolute errors of (41) satisfy that
| I 2 [ f ] Q n [ I 2 ] | = O ( ω 3 n 7 / 4 ) .
Proof. 
As is known, the error of the n-point generalized Gaussian Laguerre quadrature rule to the integral [31] 0 + x γ e x φ ( x ) d x is
E n [ φ ] = n ! Γ ( n + γ + 1 ) ( 2 n ) ! φ ( 2 n ) ( ξ ) , ξ ( 0 , + ) .
By substituting this into (26) and (41), it is not difficult to derive that
| I 2 [ f ] Q n [ I 2 ] | = O ( r 2 n 1 1 / 6 ) = O ( ω 3 n 7 / 4 ) ,
where we used r = 2 3 ω 3 / 2 . □
For the case b = + , after applying the Gaussian quadrature rule to (34), we get the numerical algorithm for I 2 [ f ] such that
Q n [ I 2 ] : = 6 3 9 π Γ ( 5 / 6 ) k = 1 n w k ( 0 ) e i r r Φ ^ 1 / 3 + 1 i r u k ( 0 ) + e i r r Φ ^ 1 / 3 1 + i r u k ( 0 ) + 2 3 3 r 2 / 3 9 π Γ ( 1 / 6 ) k = 1 n w k ( 0 ) e i r r Φ ^ 1 / 3 + 1 i r u k ( 0 ) + e i r r Φ ^ 1 / 3 1 + i r u k ( 0 ) ,
where Φ ν ± ( x ) is defined in (42). As ω + , the error of the quadrature rule Q n [ I 2 ] satisfies the estimate
| I 2 [ f ] Q n [ I 2 ] | = O ( ω 3 n 7 / 4 ) .
In the following, we provide a numerical test of the error estimates of the complex integration method for I 2 [ f ] = 1 b x α f ( x ) Ai ( ω x ) d x . We select f ( x ) = 1 / ( 1 + x 2 ) as the Runge function, which is analytic in the complex plane ( z ) 1 but has two poles at z = ± i . In Figure 5 and Figure 6, the absolute errors of Q n [ I 2 ] for various values of ω are depicted in the left panel, while the absolute errors scaled by the factor ω 3 n + 7 / 4 are presented in the right panel. It is evident that the accuracy of Q n [ I 2 ] improves as ω increases, and the convergence rate is observed to be ω 3 n 7 / 4 . This observation aligns with the theoretical results stated in Theorem 4.

4. Numerical Examples

This section is dedicated to presenting the efficiency and accuracy of the proposed method by several typical numerical examples. The exact values used were obtained from Mathematica 12.1 with high-precision arithmetic, whereas Matlab 2016 was used to compute the numerical results of the proposed method.
Example 1. 
We consider the highly oscillatory integral over a finite interval given by
I [ f ] = 0 5 x 1 / 2 sin x Ai ( ω x ) d x .
After dividing the interval, we compute the integral over [ 0 , 1 ] using the Filon-type method Q m [ I 1 ] described in Section 2, and the integral over [ 1 , 5 ] using the complex integration method Q n [ I 2 ] described in Section 3. The exact values of I [ f ] in the bottom row of Table 1 are computed by  Mathematica 12.1 with 16-digit precision. The absolute errors are displayed in Table 1, which shows the efficiency of the proposed method.
Example 2. 
Now, let us consider the highly oscillatory integral over the semi-infinity interval
I [ f ] = 0 + x 1 / 2 1 100 + x 2 Ai ( ω x ) d x .
Similarly, we compute the integral over [ 0 , 1 ] by using the Filon-type method described in Section 2 and the integral over [ 1 , + ) by using the complex integration method described in Section 3. The absolute errors are displayed in Table 2. The bottom line of Table 2 gives the exact values of I [ f ] obtained from the software  Mathematica  12.1 with 16-digit precision.
From Table 1 and Table 2, it can be readily concluded that, for any fixed ω , the accuracy of the proposed method improves significantly with an increase in n. The method requires only a small number of function evaluations to yield highly accurate approximations of I [ f ] . Furthermore, for a fixed value of n, the accuracy of the method improves with increasing values of ω . The numerical examples presented in each section underscore the efficiency of our approach. For comparison, we provide the absolute errors of the proposed method for the non-oscillatory integrals 0 b x α f ( x ) Ai ( x ) d x (see the second column of Table 1 and Table 2). The results illustrate the efficiency and accuracy of the proposed method.

5. Conclusions

In conclusion, this study addresses the numerical computation of highly oscillatory integrals involving the Airy function over both finite and semi-infinite intervals. The proposed method, which employs a partitioning strategy and complex integration techniques, demonstrates efficiency and accuracy in handling the challenges posed by the oscillatory nature of the integrands. This approach is easy to implement, and rigorous error analyses support its reliability. Numerical examples are provided to validate the theoretical results and highlight the accuracy of the methods. Furthermore, it is worth noting that the proposed method in this work achieves higher asymptotic rates concerning the frequency ω compared to those discussed in the previously mentioned references.
In the implementation of our methods, the function f is required to be analytic in the complex plane { z | 1 ( z ) b } . However, this requirement may be overly restrictive. In situations where f ( x ) is only differentiable but not analytic, it could be worthwhile to explore alternative approaches, such as employing suitably designed asymptotic methods or Filon-type methods. These considerations will be the focus of our future work.

Author Contributions

Methodology, G.L.; software, G.L. and B.L.; formal analysis, G.L., Z.X. and B.L.; writing—original draft preparation, G.L.; writing—review and editing, Z.X. and B.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Natural Science Foundation of Henan (grand number 232300420118), Guizhou Province 2022 Philosophy and Social Science Planning Major Project (grant number: 22GZZB09), and Guizhou Provincial Science and Technology Plan Project (grant number: QKHJC-ZK[2021]YB018).

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Cai, W. Computational Methods for Electromagnetic Phenomena: Electrostatics in Solvation, Scattering, and Electron Transport; Cambridge University Press: Cambridge, UK, 2013. [Google Scholar]
  2. Colton, D.; Kress, R. Integral Equation Methods in Scattering Theory; SIAM: Philadelphia, PA, USA, 2013. [Google Scholar]
  3. Iserles, A. On the numerical quadrature of highly-oscillating integrals I: Fourier transforms. IMA J. Numer. Anal. 2004, 24, 365–391. [Google Scholar] [CrossRef]
  4. Iserles, A. On the numerical quadrature of highly-oscillating integrals II: Irregular oscillators. IMA J. Numer. Anal. 2005, 25, 25–44. [Google Scholar] [CrossRef]
  5. Iserles, A.; Nørsett, S. Efficient quadrature of highly oscillatory integrals using derivatives. Proc. R. Soc. A 2005, 461, 1383–1399. [Google Scholar] [CrossRef]
  6. Patidar, K.; Shikongo, A. The application of asymptotic analysis for developing reliable numerical method for a model singular perturbation problem. J. Numer. Anal. Ind. Appl. Math. 2007, 2, 193–207. [Google Scholar]
  7. Chung, K.; Evans, G.; Webster, J. A method to generate generalized quadrature rules for oscillatory integrals. Appl. Numer. Math. 2000, 34, 85–93. [Google Scholar] [CrossRef]
  8. Debnath, P.; Srivastava, H.; Chakraborty, K.; Kumam, P. Advances in Number Theory and Applied Analysis; World Scientific: Singapore, 2023. [Google Scholar]
  9. Engquist, B.; Fokas, A.; Hairer, E.; Iserles, A. Highly Oscillatory Problems; Cambridge University Press: Cambridge, UK, 2009. [Google Scholar]
  10. Hazarika, B.; Acharjee, S.; Srivastava, H. Advances in Mathematical Analysis and Its Applications; CRC Press: New York, NY, USA, 2022. [Google Scholar]
  11. Kang, H. Numerical integration of oscillatory Airy integrals with singularities on an infinite interval. J. Comput. Appl. Math. 2018, 333, 314–326. [Google Scholar] [CrossRef]
  12. Olver, S. Numerical approximation of vector-valued highly oscillatory integrals. BIT Numer. Math. 2007, 47, 637–655. [Google Scholar] [CrossRef]
  13. Olver, F.; Lozier, D.; Boisvert, R.; Clark, C. NIST Handbook of Mathematical Functions; Cambridge University Press: Cambridge, UK, 2010. [Google Scholar]
  14. Xu, Z.; Xiang, S. Numerical evaluation of a class of highly oscillatory integrals involving Airy functions. Appl. Math. Comput. 2014, 246, 54–63. [Google Scholar] [CrossRef]
  15. Levin, D. Fast integration of rapidly oscillatory functions. J. Comput. Appl. Math. 1996, 67, 95–101. [Google Scholar] [CrossRef]
  16. Xiang, S.; Wang, H. Fast integration of highly oscillatory integrals with exotic oscillators. Math. Comp. 2010, 79, 829–844. [Google Scholar] [CrossRef]
  17. Kang, H.; Wang, H. Asymptotic analysis and numerical methods for oscillatory infinite generalized Bessel transforms with an irregular oscillator. J. Sci. Comput. 2020, 82, 29. [Google Scholar] [CrossRef]
  18. Olver, S. Numerical Approximation of Highly Oscillatory Integrals. Ph.D. Thesis, University of Cambridge, Cambridge, UK, 2007. [Google Scholar]
  19. Filon, L. On the quadrature formula for trigonometric integrals. Proc. R. Soc. Edinb. 1928, 49, 38–47. [Google Scholar] [CrossRef]
  20. Wang, H. A unified framework for asymptotic analysis and computation offinite Hankel transform. J. Math. Anal. Appl. 2020, 483, 123640. [Google Scholar] [CrossRef]
  21. López, J.; Temme, N. Two-point Taylor expansions of analytic functions. Stud. Appl. Math. 2002, 109, 297–311. [Google Scholar] [CrossRef]
  22. Oreshkin, B. MeijerG, MATLAB Central File Exchange. 2023. Available online: https://www.mathworks.com/matlabcentral/fileexchange/31490-meijerg (accessed on 19 December 2023).
  23. Gradshteyn, I.; Ryzhik, I. Table of Integrals, Series, and Products, 8th ed.; Academic Press: Cambridge, MA, USA, 2015. [Google Scholar]
  24. Lang, S. Complex Analysis, 4th ed.; Springer: New York, NY, USA, 1999. [Google Scholar]
  25. Chen, R. Numerical approximations to integrals with a highly oscillatory Bessel kernel. Appl. Numer. Math. 2012, 62, 636–648. [Google Scholar] [CrossRef]
  26. Zaman, S.; Siraj-ul-Islam; Khan, M.; Ahmad, I. New algorithms for approximation of Bessel transforms with high frequency parameter. J. Comput. Appl. Math. 2022, 299, 113705. [Google Scholar] [CrossRef]
  27. Huybrechs, D.; Vandewalle, S. On the evaluation of highly oscillatory integrals by analytic continuation. SIAM J. Numer. Anal. 2006, 44, 1026–1048. [Google Scholar] [CrossRef]
  28. Milovanović, G. Numerical calculation of integrals involving oscillatory and singular kernels and some applications of quadratures. Comput. Math. Appl. 1998, 36, 19–39. [Google Scholar] [CrossRef]
  29. Glaser, A.; Liu, X.; Rokhlin, V. A fast algorithm for the calculation of the roots of special functions. SIAM J. Sci. Comput. 2007, 29, 1420–1438. [Google Scholar] [CrossRef]
  30. Hale, N.; Trefethen, L. Chebfun and numerical quadrature. Sci. China Math. 2012, 55, 1749–1760. [Google Scholar] [CrossRef]
  31. Evans, G.; Chung, K. Some theoretical aspects of generalised quadrature methods. J. Complex. 2003, 19, 272–285. [Google Scholar] [CrossRef]
Figure 1. Plot of the Airy function Ai ( ω x ) for ω = 20 (left) and ω = 50 (right).
Figure 1. Plot of the Airy function Ai ( ω x ) for ω = 20 (left) and ω = 50 (right).
Mathematics 12 00377 g001
Figure 2. The absolute errors of Q m [ I 1 ] with f ( x ) = e 1 / x 2 and varying ω = 10 2 : 20 : 10 3 . The dashed lines depict the decay rates according to O ( ω m α 1 ) .
Figure 2. The absolute errors of Q m [ I 1 ] with f ( x ) = e 1 / x 2 and varying ω = 10 2 : 20 : 10 3 . The dashed lines depict the decay rates according to O ( ω m α 1 ) .
Mathematics 12 00377 g002
Figure 3. The complex integration paths for the integrals in (24). The left one is for integrals with Φ ν + ( x ) , while the right one is for integrals with Φ ν ( x ) .
Figure 3. The complex integration paths for the integrals in (24). The left one is for integrals with Φ ν + ( x ) , while the right one is for integrals with Φ ν ( x ) .
Mathematics 12 00377 g003
Figure 4. The complex integration paths for the integrals in (24) when b = + . The left one is for integrals with Φ ν + ( x ) , while the right one is for integrals with Φ ν ( x ) .
Figure 4. The complex integration paths for the integrals in (24) when b = + . The left one is for integrals with Φ ν + ( x ) , while the right one is for integrals with Φ ν ( x ) .
Mathematics 12 00377 g004
Figure 5. The absolute errors of Q n [ I 2 ] (left) and the absolute errors scaled by ω 3 n + 7 / 4 (right) for the integral I 2 [ f ] = 1 2 x α f ( x ) Ai ( ω x ) d x with f ( x ) = 1 / ( 1 + x 2 ) , where ω runs over [ 10 , 100 ] . The exact value of I 2 [ f ] is approximated by the software Mathematica 12.1 with high-precision arithmetic.
Figure 5. The absolute errors of Q n [ I 2 ] (left) and the absolute errors scaled by ω 3 n + 7 / 4 (right) for the integral I 2 [ f ] = 1 2 x α f ( x ) Ai ( ω x ) d x with f ( x ) = 1 / ( 1 + x 2 ) , where ω runs over [ 10 , 100 ] . The exact value of I 2 [ f ] is approximated by the software Mathematica 12.1 with high-precision arithmetic.
Mathematics 12 00377 g005
Figure 6. The absolute errors of Q n [ I 2 ] (left) and the absolute errors scaled by ω 3 n + 7 / 4 (right) for the integral I 2 [ f ] = 1 + x α f ( x ) Ai ( ω x ) d x with f ( x ) = 1 / ( 1 + x 2 ) , where ω runs over [ 10 , 100 ] . The exact value of I 2 [ f ] is approximated by the software Mathematica 12.1 with high-precision arithmetic.
Figure 6. The absolute errors of Q n [ I 2 ] (left) and the absolute errors scaled by ω 3 n + 7 / 4 (right) for the integral I 2 [ f ] = 1 + x α f ( x ) Ai ( ω x ) d x with f ( x ) = 1 / ( 1 + x 2 ) , where ω runs over [ 10 , 100 ] . The exact value of I 2 [ f ] is approximated by the software Mathematica 12.1 with high-precision arithmetic.
Mathematics 12 00377 g006
Table 1. Errors of the proposed numerical method for (48).
Table 1. Errors of the proposed numerical method for (48).
n ( = m ) ω = 1 ω = 10 ω = 20 ω = 40 ω = 80 ω = 160
1 3.2 × 10 2 2.3 × 10 3 8.2 × 10 4 2.9 × 10 4 1.0 × 10 4 3.6 × 10 5
2 5.9 × 10 3 5.7 × 10 6 5.4 × 10 7 4.9 × 10 8 4.3 × 10 9 3.8 × 10 10
3 1.5 × 10 3 6.4 × 10 8 9.3 × 10 9 9.8 × 10 10 9.4 × 10 11 8.6 × 10 12
4 3.1 × 10 4 3.9 × 10 10 2.3 × 10 11 1.1 × 10 12 4.9 × 10 14 2.2 × 10 15
5 1.3 × 10 5 2.9 × 10 18 1.1 × 10 12 1.8 × 10 14 2.6 × 10 16 1.2 × 10 18
Exact values0.57864500.01372620.00519140.00182970.00067400.0002367
(Mathematica)008553076197215839052109309685469185412476298767
Table 2. Errors of the proposed numerical method for (49).
Table 2. Errors of the proposed numerical method for (49).
n ( = m ) ω = 1 ω = 10 ω = 20 ω = 40 ω = 80 ω = 160
1 4.3 × 10 4 1.6 × 10 6 5.2 × 10 7 1.8 × 10 7 6.4 × 10 8 2.2 × 10 8
2 4.0 × 10 4 5.5 × 10 10 3.6 × 10 11 3.9 × 10 12 2.5 × 10 13 2.8 × 10 13
3 1.9 × 10 4 6.1 × 10 12 1.2 × 10 13 1.4 × 10 14 1.7 × 10 13 5.9 × 10 14
4 6.0 × 10 5 2.2 × 10 13 1.9 × 10 14 2.1 × 10 15 2.7 × 10 14 5.8 × 10 14
5 4.3 × 10 6 2.8 × 10 14 1.9 × 10 14 2.1 × 10 15 2.7 × 10 14 5.8 × 10 14
Exact values0.01088550.00344280.00243440.00172140.00121720.0008607
(Mathematica)105680928487379187277255424640434363944521232015
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Liu, G.; Xu, Z.; Li, B. An Efficient Quadrature Rule for Highly Oscillatory Integrals with Airy Function. Mathematics 2024, 12, 377. https://doi.org/10.3390/math12030377

AMA Style

Liu G, Xu Z, Li B. An Efficient Quadrature Rule for Highly Oscillatory Integrals with Airy Function. Mathematics. 2024; 12(3):377. https://doi.org/10.3390/math12030377

Chicago/Turabian Style

Liu, Guidong, Zhenhua Xu, and Bin Li. 2024. "An Efficient Quadrature Rule for Highly Oscillatory Integrals with Airy Function" Mathematics 12, no. 3: 377. https://doi.org/10.3390/math12030377

APA Style

Liu, G., Xu, Z., & Li, B. (2024). An Efficient Quadrature Rule for Highly Oscillatory Integrals with Airy Function. Mathematics, 12(3), 377. https://doi.org/10.3390/math12030377

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop