On the Convolution Quadrature Rule for Integral Transforms with Oscillatory Bessel Kernels

Abstract: Lubich’s convolution quadrature rule provides efficient approximations to integrals with special kernels. Particularly, when it is applied to computing highly oscillatory integrals, numerical tests show it does not suffer from fast oscillation. This paper is devoted to studying the convergence property of the convolution quadrature rule for highly oscillatory problems. With the help of operational calculus, the convergence rate of the convolution quadrature rule with respect to the frequency is derived. Furthermore, its application to highly oscillatory integral equations is also investigated. Numerical results are presented to verify the effectiveness of the convolution quadrature rule in solving highly oscillatory problems. It is found from theoretical and numerical results that the convolution quadrature rule for solving highly oscillatory problems is efficient and high-potential.


Introduction
Highly oscillatory integrals (HOI) arise frequently in antenna problems involving Sommerfeld integrals (see [1,2]), computation of mutual impedance between conductors (see [3,4]) and many other oscillatory problems (HOP).Generally, an oscillatory integral can be written as Here W(•) denotes a highly oscillatory function and f (•) is slowly varied.Due to the high oscillation, classical quadrature rules (e.g., Newton-Cotes and Gauss rules ) are often ineffective, and calculation of this class of integrals is deemed to be a challenging problem ([5]).
Past decades witness a rapid development of researches on calculation of HOIs.Based on Filon's idea ( [6]), Iserles and Nørsett developed the Filon-type method by approximating the slowly varied function by its Hermite interpolant.Both of theoretical and numerical results manifested that this method enjoyed high-order convergence rates with respect to the frequency ( [7]).To get stable and fast algorithms, Domínguez,et al. ([8]), and Xiang, et al. ( [9]), proposed the Clenshaw-Curtis-Filon-type method, respectively, which enjoyed extensive applications at present.Although Filon's methodology leads to many efficient algorithms, most of them suffer to complicate computation of moment integrals.An alternative way to addressing this problem is transforming the integral interval into the complex plane, which derives the numerical steepest descent method.In [10], the complex integration method in the standard case was discussed extensively.In [11], the general form and corresponding error analysis were studied.This method was extended to the case of HOIs on semi-finite intervals in [12].It is notable that the numerical steepest descent method with Gauss-Laguerre quadrature may not provide satisfactory solutions in practice (see [13]).Therefore, computation of transformed integrals is still an interesting topic.
Another important quadrature rule for HOIs is Levin's method ( [14]).It enjoys a wide application for its being free of computing complex moments and less restrictions to integrands.The spirit of this method is transforming the integration problem into an ODE, and solving this equation by collocation methods.It is well-know that implementing of Levin methods comes down to efficient solutions of linear systems, which are often singular and dense.In [15], an SVD solver for the ill-conditioned system was presented.In [16], Olver developed a moment-free method by using the shifted GMRES.Recently, by employing the property of Chebyshev polynomials and preconditioners, a sparse and well-conditioned Levin method was constructed in [17].
There are many other important methods for calculating HOIs, for example, the homotopy perturbation method ( [18]), the generalized quadrature rule ( [19]), the extrapolation method ( [20]).For simplicity, we omit the details.It is quite unexpected that little attention has been paid to the convolution quadrature rule (CQ) for HOIs ( [21,22]), especially its asymptotic property in the case of high oscillation.In fact, CQ is well-known for its efficient in evaluating convolution integrals, and oscillatory integrals of convolution-type play significant roles in solving oscillatory and evolutionary problems (see [23,24]).Therefore, it is a meaningful issue to study CQ for HOIs.
Consider integral transforms with Bessel kernels as with m ≥ 0 being an integer and ω 1, and oscillatory Volterra integral equations as where f (•) is sufficiently smooth, f (0) = 0, and u(t) is unknown.In this paper, we are devoted to studying convergence property of CQ with respective to the frequency for solving above two problems.The same models have been considered in [25,26], where authors concluded that Filon-type methods enjoyed the property that the higher the oscillation, the better the approximation.In the following, we will find CQ share a similar property as Filon-type methods, and even better when they are applied to solving highly oscillatory integral equations.The remaining parts are organized as follows.
In Section 2, we briefly review CQ and give the convergence analysis.A modified rule is also proposed in this part.Then we study CQ for solving Volterra integral equations with highly oscillatory kernels in Section 3. Some numerical experiments are carried out in Section 4 to verify our given results.

Convergence of the Convolution Quadrature Rule
In this section, we revisit Lubich's convolution quadrature firstly.Then the convergence property with respect to the frequency is studied.In [21,22], Lubich proposed an algorithm for computing the following integral, Let F(•) denote the Laplace transform of f (•) and satisfy By the definition of Laplace transform, it follows that where Γ is a curve locates in the analytic region of F(s) and goes from ∞ • e −i(π−ϕ) to ∞ • e i(π−ϕ) .Substituting ( 5) into (4) gives x 0 e λt g(x − t)dtdλ.
Noting that y(x) = x 0 e λt g(x − t)dt satisfies the initial value problem Defining the grid {t n := nh, n = 0, 1, ...}, we can approximate y(x) by where 8) and summing give where it follows that Since F(s) is analytic in the inside region of the curve Γ, we have, by Cauchy's integral formula (see [27], p. 32), Therefore, it follows that Here the coefficient corresponding to ζ n denotes an approximation to the integral (4) at Then CQ for ( 4) is defined as By the definition of coefficients of Taylor expansions, we get where ρ is sufficiently small such that the disc |z| ≤ ρ falls in the analytic region of F (δ(ζ)/h) .Letting z = ρe iθ , we obtain The last equal sign works due to F δ(ρe iθ ) h e −inθ is 2π−periodic.This leads to Its computation complexity is O(N log N) by FFT.In this paper, we adopt L = 10 N, ρ N = √ , = 10 −16 to guarantee a precision of order O( √ ) in (19).
Remark 1.In [28], convolution quadrature weights were rewritten as Here φ j (t) = e −t t j j! for backward differentiation formula of order 1 (BDF1), and φ j (t) = e −3t/2 for BDF2 with H j (•) denoting the jth Hermite polynomial.Recurrence relations of frequently-used bases can be found in Table 1.

Table 1. Recurrence relations for the CQ basis functions.
Scheme Initial Basis (φ −n = 0, n ≥ 1) Recurrence for Basis Functions (j ≥ 1) Existing convergence analysis of CQ often restricts to the property with respect to the stepsize.For example, setting the convergence rate of CQ is as follows Then we have where the constant C does not depend on the stepsize h.
One point which should be remarkable is that A-stable linear multistep methods for solving the initial value problem (7) are more reliable in actual computation.Therefore, in this paper, we make use of BDF2 in numerical experiments.
In the following parts, we study the convergence property of CQ for Here J m (•) denotes the first kind Bessel function of order m with m being a nonnegative integer.Defining f (s) := f (b − s), we transform (1) into Here the Laplace transform of J m (ωs) can be represented as (see [29], p. 1024) Before elaborating on the convergence property of CQ, the following lemma should be presented first.( Furthermore, if the interval S does not contain origin, then it follows This lemma can be proved by integration by parts and we omit the detail.Now the main theorem of this section follows.
Then there exists a positive constant C independent of ω, such that, as the frequency ω tends to infinity, CQ satisfies Proof.Denote by δ j the coefficients of

and define
Then I[ f ] can be represented as a Dunford-Taylor integral By noting we have This implies Define φ 1 (t) := e −ts h f (x) − f (x − t) and φ 2 (t) := e −ts h f (x).By recurrence relations and derivatives of Bessel functions, we get According to φ 1 (0) = 0 and Lemma 1, it follows where the constant C does not depend on ω.This completes the proof.
This theorem verifies the accuracy of CQ will increase as the frequency tends to infinity.According to the proof of Theorem 2, we can eliminate low order terms in (35) and get a modified convolution quadrature rule.
According to ( [30], p. 681), we have ) Here s (2) µ,ν (x) denotes the Lommel function of the second kind, and can be efficient computed by asymptotic expansions ( [31]).This implies Analogy to the proof of Theorem 2, we immediately arrive at the following results.
Then there exists a positive constant C independent of ω, such that, as the frequency ω tends to infinity, MCQ1 satisfies Remark 2. By Lubich's methodology of eliminating low order terms, we also obtain a modified convolution quadrature rule of the second kind for (1) (MCQ2).Although these two modified quadrature rules share the same convergence rates with respect to the stepsize, their convergence properties are quite different in the case of calculation of HOI.We will illustrate this phenomenon in Section 4.

Application to a Volterra Equation
In literature, convolution quadrature rules are important tools for solving Volterra equations with convolution kernels ( [32]).Many numerical experiments show they are efficient for solving some highly oscillatory Volterra integral equations (HOVIE).Consider HOVIE (3), where f (•) is sufficiently smooth and f (0) = 0. Let be a uniform grid on [0, T] with spacing h := T/N.This equation arises from acoustic scattering problems (see [33]).By applying CQ to (3), we get Here u j denotes the numerical solution ũ(x) at x j and f k := f (x k ).Once the initial value u 0 is known, the numerical solution at the uniform grid I N can be obtained by solving the linear system (42), and the numerical solution on [0, T] can be written as Following the methodology from [24,26], we establish the convergence analysis by expressing the error function in terms of moments with highly oscillatory kernels.So let us consider some integrals involving Bessel kernels.

Lemma 2. Define the functional T
with x ∈ [0, T] and f ∈ C([0, T]).Then ∞−norm of the functional T x satisfies where C is a constant independent of ω.
Proof.It is easy to show By the variable transformation t = ωs, we have According to the asymptotic expansions of Bessel functions (see [34], p. 228), there exists Z > 0 such that for any z > Z, and a constant C * , we have This implies Therefore, there exists a constant C, such that This completes the proof.
Proof.The variable transformation s = ωt gives According to ( [34], p. 242), we have By the asymptotic expansion of Bessel functions, we completes the proof.Now we get the convergence property of CQ for the numerical solution to HOVIE.
Theorem 4. Suppose that f ∈ C 3 [0, T], then the convolution quadrature rule for solving (3) introduces a unique numerical solution ũ, and satisfies where C is a constant independent of ω.
Proof.By noting and with x ∈ I h , we get the error equation where (x) := ũ(x) − u(x).By Remark 1 and Lemma 1, we have Therefore, the remaining work is proving F (s)u(x) − F (s h )u(x) behaves as O(ω −3/2 ).A similar process to the proof of Theorem 2 gives where φ 1 (t) := e −ts h u(x) − u(x − t) and φ 2 (t) := e −ts h u(x).Consider the integrals By using integration by parts we have According to [35], the exact solution to (3) can be written as Furthermore, we have A direct calculation implies as ω goes to infinity.By Lemmas 2 and 3, we obtain With the help of Lemma 1, we have Ĩ = O(ω −1/2 ).It follows that Combining Equations ( 58) and (66) gives where C is a constant independent of ω.This completes the proof.

Numerical Results
In this section, we present some numerical results to verify given estimates in previous sections.All experiments are performed in MATLAB 2013b.
As first example, we consider the following HOI, In Figures 1-3, we show the convergence rates of three convolution quadrature rules.Slowly varied lines in these figures manifest that given asymptotic orders in Section 2 are optimal.Absolute errors of these methods are given in Tables 2 and 3.The numerical results illustrate MCQ1 is much more efficient than other two methods in computing HOI.In the example, we consider application of CQ to solving HOVIE (3).Firstly, by letting T = 2 and N = 20, we give the computed solution in Table 4 with various ω.Absolute errors listed in this table show CQ shares the same property as Filon methods ( [24,26]), that is, the higher the oscillation, the better the approximation.Then we compare these two methods in Figure 4, where we can learn CQ behaves better than Filon methods.

Conclusions
The above theoretical and numerical results contribute to the study on the convergence property of CQ for solving HOP.The theoretical results in Section 2 reveal the convergence rate of CQ with respect to the frequency, that is, CQ converges in negative powers of ω as ω goes to infinity.Among them, the new modified rule (MCQ1) enjoys the fastest convergence rate.When we apply CQ to solving HOVIE, similar phenomenon is detected and analyzed.The numerical results in Section 4 show given convergence orders in Section 2 are optimal.In addition, this paper merely opens a window to the convergence theory of CQ for HOP, much work on various versions of CQ, such as Runge-Kutta CQ ( [36]), Fourier CQ ( [37]), and so forth, is needed in the future.

Lemma 1 .
Suppose f (•) is continuous on the interval S, where S may be a closed interval on the positive real axis or [a, ∞) for some a ≥ 0. Assume S | f (t)|dt, S | f (t)|dt exist.Then for any v > 0 and sufficiently large ω, there exists a constant C, which does not depend on the frequency ω, such that S f (t)J m (ωt)dt ≤ Cω −1 .

Figure 1 .
Figure 1.Asymptotic convergence rates of CQ for

Figure 2 .
Figure 2.Asymptotic convergence rates of MCQ1 for

Figure 4 .
Figure 4. Comparisons between Filon methods and CQ for solving Volterra equations with f (x) = x 1 + x 2 .

Table 2 .
Comparisons of quadrature rules for

Table 3 .
Comparisons of quadrature rules for