Numerical analysis of the spectrum for the highly oscillatory integral equation with weak singularity

Abstract

In this paper, we focus on the spectra numerical computation for the integral equation with the absolute oscillation and power-law or logarithmic singularity. Finite section method is applied to transform the integral equation to an algebraic eigenvalue problem. The entries of the coefficient matrix appearing in the bivariate highly oscillatory singular integrals can be represented explicitly in Gamma or the exponential integral functions. The decay rate of the entries is established to construct the truncation scheme. Then the infinite algebraic eigenvalue problem can be simplified to be the finite one. The corresponding error of the infinite and finite algebraic systems is also bounded. Finally, the numerical experiments are provided to illustrate the theoretical analysis.

Introduction

We focus on the spectra of the integral equation

$$\mathcal{F}\left[f\right]≔{\int }_{-1}^{1}f\left(y\right)k\left(|x-y|\right)\text{d}y=\lambda f\left(x\right),x\in [-1,1],$$

where the kernel function $k\left(|x-y|\right)$ is assumed to be even function, such as

$$k\left(|x-y|\right)=\frac{{\mathrm{e}}^{\mathrm{i}\omega |x-y|}}{{|x-y|}^{\alpha }},0\le \alpha <\frac{1}{2},$$$$k\left(|x-y|\right)={\mathrm{e}}^{\mathrm{i}\omega |x-y|}ln|x-y|,$$

in which $\omega \ge 0$ is an oscillatory parameter. Once $\omega =0$, the integral equation becomes one with weak singularity. The integral operator $\mathcal{F}$ has arisen in the singular integrals of the lower dimensional varieties, singular Radon transform [19], [20], and the application of the boundary integral equation of acoustic, electromagnetic wave scattering [3], [4]. In particular, its numerical computation is also proposed as an open problem in Section 6 of Brunner et al. [1]. The integral operator $\mathcal{F}$ of (1.1) with (1.2) or (1.3) in ${\mathbb{L}}_2[-1,1]$ is compact which inspires the computational way of finite projection operator approximation [14], [21]. The compactness of the operators also indicates that the discrete spectrum of the operator accumulate to the origin.

Highly oscillatory integral arises in many engineering and numerical analysis fields, for instance, acoustic scattering, magnetics, quantum mechanics and so on. In the last decades, much efforts are devoted to the quadrature of the highly oscillatory integral. The asymptotic expansion method is fundamental to the numerical computations but only working for the sufficiently large oscillatory parameter [5]. Filon type method is applied extensively with high accuracy when the moment function can be calculated exactly [9], [12]. Avoiding the computation of moment, Levin method is constructed for the oscillatory integral without the stationary point [15], [16]. Numerical steepest descent method turns the oscillation integral to the non-oscillatory one by changing the integral path for analytic function [11]. The oscillatory integral with Bessel oscillation is explored in [22]. Furthermore, the asymptotic expansion and Filon type methods have been extended to the highly oscillatory singular integral recently [8].

Weak singularity and high oscillations hinder the accurately fast numerical computation for the integral equation. Since the analytic solution is difficult to obtain, numerical computation for (1.1) with (1.2) or (1.3) is a better option. The theoretical analysis of the numerical approximation for a compact operator has been studied in [18]. In [2], several numerical computation methods are presented which are the first step to the numerical computation for the Fox–Li operator. The asymptotic analysis of the spectra of the integral operator with the absolute oscillator out of the singularity is formulated in [1]. In [7], one finite section method with the modified Fourier basis has been presented to analyse the spectra for the Fox–Li integral operator ${\mathrm{e}}^{\mathrm{i}\omega {(x-\gamma y)}^2}$, where $0<\gamma \le 1$ and $\omega \gg 1$. For the integral equation (1.1) with (1.2), the modified Fourier basis is firstly used to approximate the unknown in which the entries of the matrix are represented in the hypergeometric functions [10]. Here, we are interested in introducing the finite section method to analyse the numerical spectrum (1.1), especially the influence of the singular parameter $\alpha$ for the distribution of the spectrum. For the absolute oscillator integral equation with power-law or logarithmic singularities, we construct a unified numerical method which can be implemented easily by considering the affection of two oscillatory factors $m$, $n$ and $\omega$.

In this paper, a finite section method based on the modified Fourier basis is applied to compute the spectrum of the highly oscillatory integral equation with algebraic singularity. The accurate approximation of the non-periodic function facilitates the approximation of the modified Fourier basis for the integral equation [13]. During transforming the problem of the operator spectrum to the algebraic problem, one infinite coefficient matrix is derived whose entries are in the form of bivariate highly oscillatory singular integral. By some algebra computation, the double integral is simplified to be the explicit expression. The asymptotic properties of the entries are proved to construct the truncation strategy. Finally, the trajectories of the spectra are shown for the highly oscillatory singular integral operators.

The paper is organised as follows. In Section 2, for the integral equation (1.1) with (1.2), the modified Fourier basis is applied to obtain the explicit expressions for the resulting matrix. The asymptotic orders are proved for the entries with increasing oscillatory parameters. We construct the truncation scheme based on the asymptotic estimates. Then the numerical results are provided to show how the singular and oscillatory parameters affect the spectra. The similar analysis for the integral (1.1) with (1.3) is explored in Section 3. Section 4 concludes the paper.

Before giving the development of the method in details, we present a basic frame of the method for the integral equation (1.1) with (1.2) or (1.3). Applying an orthogonal system ${\left\{{\psi }_n\left(x\right)\right\}}_{n\ge 0}$ as the trial basis to approximate the unknown eigenfunctions $f\left(x\right)$,

$$f\left(x\right)=\sum _{n=0}^{\infty }{\hat{f}}_n{\psi }_n\left(x\right),$$

with the same test function as the trial basis, results in an algebraic system

$$\mathcal{Q}\hat{\mathbf{f}}=\lambda \hat{\mathbf{f}},$$

where

$$\mathcal{Q}≔{\left(q_{m,n}\right)}_{m,n\ge 0},q_{m,n}={\int }_{-1}^1{\psi }_m\left(x\right)\left({\int }_{-1}^1{\psi }_n\left(y\right)k\left(|x-y|\right)\mathrm{d}y\right)\mathrm{d}x.$$

For the inner integral, it holds that

$${\int }_{-1}^1{\psi }_n\left(y\right)k\left(|x-y|\right)\mathrm{d}y={\int }_{-1}^x{\psi }_n\left(y\right)k\left(x-y\right)\mathrm{d}y+{\int }_x^1{\psi }_n\left(y\right)k\left(y-x\right)\mathrm{d}y={\int }_0^{1+x}{\psi }_n(x-t)k\left(t\right)\mathrm{d}t+{\int }_0^{1-x}{\psi }_n(x+t)k\left(t\right)\mathrm{d}t.$$

Thus, we have

$$q_{m,n}={\int }_{-1}^1{\int }_0^{1+x}{\psi }_m\left(y\right){\psi }_n(x-t)k\left(t\right)\mathrm{d}t\mathrm{d}x+{\int }_{-1}^1{\int }_0^{1-x}{\psi }_m\left(y\right){\psi }_n(x+t)k\left(t\right)\mathrm{d}t\mathrm{d}x={\int }_0^2\left({\int }_{t-1}^1{\psi }_m\left(x\right){\psi }_n(x-t)\mathrm{d}x+{\int }_{-1}^{1-t}{\psi }_m\left(x\right){\psi }_n(x+t)\mathrm{d}x\right)k\left(t\right)\mathrm{d}t.$$

Furthermore, if the orthogonal system ${\left\{{\psi }_n\left(x\right)\right\}}_{n\ge 0}$ is constituted of the even basis ${\left\{{\varphi }_n\left(x\right)\right\}}_{n\ge 0}$ and odd basis ${\left\{{\phi }_n\left(x\right)\right\}}_{n\ge 0}$, the eigenfunction $f$ in (1.1) is expanded as

$$f\left(x\right)=\sum _{n=0}^{\infty }{\hat{f}}_n^e{\varphi }_n\left(x\right)+\sum _{n=0}^{\infty }{\hat{f}}_n^o{\phi }_n\left(x\right).$$

Substituting the expansion into (1.1), we get an infinite-dimensional algebraic eigenvalue problem

$$\mathcal{Q}\hat{\mathbf{f}}=\left(\begin{array}{cc}\hfill \mathcal{A}\hfill & \hfill \mathcal{B}\hfill \\ \hfill \mathcal{C}\hfill & \hfill \mathcal{D}\hfill \end{array}\right)\left(\begin{array}{c}\hfill {\hat{\mathbf{f}}}^e\hfill \\ \hfill {\hat{\mathbf{f}}}^o\hfill \end{array}\right)=\lambda \left(\begin{array}{c}\hfill {\hat{\mathbf{f}}}^e\hfill \\ \hfill {\hat{\mathbf{f}}}^o\hfill \end{array}\right),$$

where the elements of the matrices are

$$\mathcal{A}:={\left(a_{m,n}\right)}_{m,n\ge 0},a_{m,n}={\int }_{-1}^1{\int }_{-1}^1{\varphi }_m\left(x\right){\varphi }_n\left(y\right)k\left(\left|x-y\right|\right)\mathrm{d}y\mathrm{d}x,$$$$\mathcal{B}:={\left(b_{m,n}\right)}_{m,n\ge 0},b_{m,n}={\int }_{-1}^1{\int }_{-1}^1{\varphi }_m\left(x\right){\phi }_n\left(y\right)k\left(\left|x-y\right|\right)\mathrm{d}y\mathrm{d}x,$$$$\mathcal{C}:={\left(c_{m,n}\right)}_{m,n\ge 0},c_{m,n}={\int }_{-1}^1{\int }_{-1}^1{\phi }_m\left(x\right){\varphi }_n\left(y\right)k\left(\left|x-y\right|\right)\mathrm{d}y\mathrm{d}x,$$$$\mathcal{D}:={\left(d_{m,n}\right)}_{m,n\ge 0},d_{m,n}={\int }_{-1}^1{\int }_{-1}^1{\phi }_m\left(x\right){\phi }_n\left(y\right)k\left(\left|x-y\right|\right)\mathrm{d}y\mathrm{d}x,$$

and ${\hat{\mathbf{f}}}^e={\left({\hat{f}}_0^e,{\hat{f}}_1^e,\dots \right)}^{\top }$, ${\hat{\mathbf{f}}}^o={\left({\hat{f}}_0^o,{\hat{f}}_1^o,\dots \right)}^{\top }$. The parities of two sets ${\left\{{\varphi }_n\right\}}_{n\ge 0}$ and ${\left\{{\phi }_n\right\}}_{n\ge 0}$ determine that the integrals in $\mathcal{B}$ and $\mathcal{C}$ are zeros. Therefore, the system (1.5) reduces to two uncoupled sub-systems

$$\mathcal{A}{\hat{\mathbf{f}}}^e=\lambda {\hat{\mathbf{f}}}^e,$$$$\mathcal{D}{\hat{\mathbf{f}}}^o=\lambda {\hat{\mathbf{f}}}^o.$$

The definitions of (1.6), (1.7) imply that $a_{m,n}=a_{n,m}$ and $d_{m,n}=d_{n,m}$ which means the matrices $\mathcal{A}$ and $\mathcal{D}$ are (complex) symmetric. In the next section, we explore the exact expressions and asymptotics of $a_{m,n}$ and $d_{m,n}$ by introducing the modified Fourier basis.