    # Bayesian Approach to Inverse Acoustic Scattering from Disks and Line Cracks with Phaseless Data

This paper is concerned with inverse acoustic scattering problem of inferring the position and shape of a sound-soft obstacle from phaseless far-field data. We propose the Bayesian approach to recover sound-soft disks and line cracks through properly chosen incoming waves in two dimensions. Given the Gaussian prior measure, the well-posedness of the posterior measure in the Bayesian approach is discussed. The Markov Chain Monte Carlo (MCMC) method is adopted in the numerical approximation and the preconditioned Crank-Nicolson (pCN) algorithm with random proposal variance is utilized to improve the convergence rate. Numerical examples are provided to illustrate effectiveness of the proposed method.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Time-harmonic inverse scattering problems have attracted extensive attention due to their numerous applications in many areas such as radar and sonar detection, geophysical prospection, medical imaging, nondestructive testing and so on. In this paper, we are interested in the inverse problem of reconstructing the location and shape of an acoustically sound-soft obstacle using phaseless far-field data.

The propagation of a time-harmonic incident field in a homogeneous and isotropic medium is governed by the Helmholtz equation

 Δuin+k2uin=0inR2, (1)

where is the wavenumber. Let be a sound-soft scatterer, which occupies a bounded subset with -smooth boundary such that the exterior of is connected. In this paper maybe a domain or a curve, which represents an extended obstacle or a crack in acoustics. The forward scattering problem is to find the scattered (perturbed) field to the Helmholtz equation

 Δusc+k2usc=0inR2∖¯D, (2)

which satisfies the Dirichlet boundary condition

 usc=−uinon∂D, (3)

 limr→∞√r(∂usc∂r−ikusc)=0,r=|x|, (4)

uniformly in all directions , . The total field is defined as in . The Sommerfeld radiating solution has an asymptotic behavior of the form

 usc(x)=eik|x|√|x|{u∞(^x)+O(1√|x|)},|x|→∞, (5)

where is called the far-field pattern at the observation direction . Note that is an analytic function with phase information. The above model also appears in the TE polarization of time-harmonic electromagnetic scattering from infinitely long and perfectly conducting cylinders.

The uniqueness, stability and inversion algorithms for recovering from phased far-field patterns have been extensively studied with one or many incoming plane and point source waves. We refer to the monographs [6, 8, 23, 22, 36] for historical remarks, an overview of recent progresses and the comparison between different approaches. In many practical applications, the phase information of the far-field pattern cannot be measured accurately compared with its modulus or intensity. For instance, in optics it is not trivial to measure the phase of electromagnetic waves incited at high frequencies. One of the essential difficulties in using phaseless far-field data lies in the translation invariance property for plane wave incidence, which we state as follows. Let be the far-field pattern corresponding to the incident plane wave ( is the incident direction) and the sound-soft obstacle . For the shifted obstacle , the corresponding far-field pattern is given by (see )

 u∞(^x;Dz,d)=eikz⋅(d−^x)u∞(^x;D,d)for all^x∈S. (6)

Hence, we get

 |u∞(^x;Dz,d)|=|u∞(^x;D,d)|for alld,^x∈S.

This implies that it is impossible to recovery the location of from the phaseless far-field pattern of a plane wave.

There has been tremendous interest in inverse scattering with phaseless data or in phase retrieval problems in optics and other physical and engineering areas (see, e.g. [1, 2, 12, 17, 20, 24, 25, 33, 34] and the references therein). In a deterministic setting where randomness are not taken into account, Kress & Rundell and Ivanyshyn & Kress proposed a Newton-type iterative approach to reconstruct the shape of sound-soft obstacles from only the modulus of the far-field pattern in [17, 20, 28]. The approach of [17, 20] was based on a pair of nonlinear and ill-posed integral equations motivated by an inverse boundary value problem for the Laplace equation with phase information [18, 29]; see also [12, 19, 27] for inverse scattering from sound-soft cracks using a single far-field pattern with phase or phaseless information. Klibanov proved unique determination of a compactly supported potential of the stationary three-dimensional Schrdinger equation from the phaseless near-field data incited by an interval of frequencies . This was later extended in  to the reconstruction of a smooth wave speed in the three-dimensional Helmholtz equation. To broke the translation invariance property, it was recently prosed in  that, phaseless far-field patterns generated by infinitely many sets of superpositions of two plane waves with different directions can be used to uniquely determine a penetrable or impenetrable scatterer. Inspired by this idea, in this paper we propose to firstly generate the phaseless data using the following superposition of two plane waves

 uinℓ(x):=eikx⋅d0+eikx⋅dℓ,ℓ=0,1,2,⋯,L, (7)

and then to recover a sound-soft disk or a line crack through the Bayesian approach. In (7), we fix and change as incident directions, due to the a priori information of the obstacle; see Theorem 3.1 for a uniqueness proof for sound-soft disks.

In recent years, the Bayesian method has received increasing attention for inverse problems [10, 16, 21, 30, 38], which also has been applied to the inverse scattering problems [3, 5, 14, 31, 32, 39] with phase far-field data. In particular, the authors of  adopt the Bayesian framework of  to shape identification problems in inverse scattering and establish a framework for proving well-posedness of the Bayesian formulation using a suitable shape parametrization and the regularity of shape derivatives. The aim of this paper is to propose the Bayesian method using a single far-field pattern without phase information.

The Bayesian method provides us a new perspective to view the inverse scattering problem in the form of statistical inferences. In this statistical approach, all parameters are random variables and the key issue is to estimate the posterior distribution of the unknown quantities based on the Bayes’ formula

 and the known prior distribution. Compared with the deterministic iterative approaches [12, 17, 20, 28], the statistic characteristics of this posterior distribution are used to estimate unknown parameters and provide a quantification of the uncertainties arising from the corresponding model predictions. Moreover, the theoretical analysis and numerical methods in the Bayesian framework are only based on the deterministic forward model. The Bayesian method could be an alternative method to overcome the challenges in the deterministic inverse problems, although it most likely leads to expensive computational cost. In this paper, the Markov chain Monte Carlo (MCMC) method [4, 11, 13] is proposed to accomplish the characterization of the posterior distribution, while the preconditioned Crank-Nicolson (pCN) algorithm  is adopted to improve the convergence rate in the iteration of MCMC method.

This paper is organized as follows. In section 2, we adapt the Bayesian framework to inverse scattering problems with phaseless data. In section 3, we exhibit numerical results for recovering a disk and a line crack. Conclusions are given in section 4.

## 2 Bayesian Framework

In this paper we want to recover an unknown sound-soft obstacle from phaseless far-field patterns corresponding to a set of superposition of two plane waves. We propose the Bayesian approach to solve this inverse scattering problem. First of all, we set a suitable parameterization of the position and the shape of an obstacle. Then, in

the Bayesian framework, we estimate the posterior distribution of the unknown obstacle parameters. By the Bayes’ theorem

[30, 38], the posterior distribution of these parameters can be obtained from the prior distribution and the likelihood function to be specified in this section. Numerically, we will adopt the Markov chain Monte Carlo method (MCMC) to get an approximation of the posterior distribution.

### 2.1 Parameterization of the obstacle

Since the boundary of the underlying obstacle is a -smooth curve, we can represent or approximate its geometrical shape by a finite set of variables

 Z:=(z1,z2,⋯,zN)⊤∈RN,N∈N0. (8)

For example, we can use four parameters to represent a line segment, where and denote respectively the two ending points, or we can use , , , , , , to approximate a star-shaped closed curve where stand for the Fourier coefficients in the truncated Fourier expansion.

Recalling the incident waves in the form of a set of superpositions of two plane waves (7), we express the (phased) far-field patterns of the scattering model (2)-(4) by

 (9)

Correspondingly, the phaseless far-field pattern are denoted by

 (10)

where is the modulus of a complex number.

### 2.2 Prior distribution

By (8), the prior distribution of the obstacle parameters depends on the distribution of . Let be independent variables with the prior density and prior measure . Then the prior density and prior measure of are respectively given by

 πpr(Z) = N∏n=1πnpr(zn), (11) μpr(dZ) = N∏n=1μnpr(dzn). (12)

In this paper we assume that

are random variables with the Gaussian distribution, that is,

 μnpr=N(mn,σn),n=1,2,⋯,N.

For simplicity, we assume that and , implying that , where

is the identity matrix.

### 2.3 Observation of far-field pattern

To bridge the parameterization (8) of the obstacle and the associated phaseless far-field data (10), we define an operator as

 |u∞(^x)|=F(Z),^x∈S, (13)

which can be regarded an abstract map from the space of obstacle parameters to the space of observation data in the continuous sense. From the well-posedness of forward scattering, is continuous but highly non-linear.

Let be a bounded linear observation operator with given by

 gm(|u∞(^x)|)=|u∞(^xm)|,m=1,2,⋯,M,

where is the set of discrete observation directions. Then the observation at the observation direction can be rephrased as

 ym=gm(|u∞(^x)|)+ηm=|u∞(^xm)|+ηm, (14)

where represents the noise polluting the observation data at the direction .

Set and . Denote by the map from the obstacle parameter space to observation space , that is,

 (15)

Our inverse problem in this paper is to determine the obstacle parameters from the observation data with the noise pollution .

### 2.4 Likelihood

We assume the observation pollution is independent of and drawn from the Gaussian distribution with the density , where is a self-adjoint positive matrix. By the observation of the phaseless data with noise (15), we can get the relationship . Define the model-data misfit function as

 Φ(Z;Y)=12|Y−G(Z)|2Ση, (16)

where . Hence, the likelihood function is given by

 ρ(Y−G(Z))=1((2π)Mdet(Ση))1/2e−Φ(Z;Y).

Furthermore, the posterior density and the posterior measure are connected to the prior measure through the Radon-Nikodym derivative , given by

 dμpostdμpr(Z)∝e−Φ(Z;Y). (17)

### 2.5 Well-posedness of Bayesian framework

The well-posedness arguments of [5, 38] can be applied to deal with our inverse scattering problem with the Bayesian approach. In our phaseless case, we are required to justify the following Assumption 2.5, relying on regularity properties of the forward operator .

###### Assumption

The map satisfies

• For every , there is an such that, for all ,

 |G(Z)|Ση≤eε∥Z∥22+^M,

where is the Euclidean norm.

• For every , there is a such that, for all with
, it holds that

 |G(Z1)−G(Z2)|Ση≤K∥Z1−Z2∥2.

We remark that there is no essential difference in proving Assumption 2.5 (i) between the phased and phaseless inverse scattering problems. The Assumption 2.5 (ii) follows directly from the triangle inequality for complex numbers and the corresponding assumption for phased inverse scattering problems. Hence, when is sound-soft disk, Assumption 2.5 can be proved following the phased arguments of ; see also [31, 39] for the proofs in the case of limited aperture data and for interior scattering problems. If is sound-soft crack, the same results can be verified by applying the Frchet differentiability with respect to the boundary of the far field operator ([26, 27]). The above assumptions together with the choice of the Gaussian prior measure (which satisfies ) lead to well-posedness of the Bayesian inverse problem, which is a result of application of Lemma 2.8, Theorem 4.1, Theorem 4.2 and Theorem 6.31 in . Before stating the well-posedness (see Theorem 2.5 below), we recall the Hellinger distance defined by

 dHell(μ1,μ2):= ⎷12∫(√dμ1dμ0−√dμ2dμ0)2dμ0, (18)

where are two measures that are absolutely continuous with respect to .

If the operator satisfies the Assumption 2.5 and the prior measure satisfies , then the posterior measure

is a well-defined probability measure on

and absolutely continuous with respect to prior measure . What’s more, the posterior measure is Lipschitz in the data , with respect to the Hellinger distance: if and are two posterior measures corresponding to data and , then there exists such that,

 dHell(μ1post,μ2post)≤C∥Y1−Y2∥2,

for all , with .

### 2.6 Preconditioned Crank-Nicolson (pCN) algorithm with random proposal variance

This subsection is devoted to the numerical approximation of the posterior distribution. We adopt the Markov chain Monte Carlo method (MCMC) [4, 11, 13] to generate a large number of samples subject to the posterior distribution. The numerical approximation of the posterior distribution of unknown obstacle parameters can be obtained by statistic analysis on these samples. The Metropolis-Hastings [15, 35] algorithm will be used to construct MCMC samples. To improve the convergence rate of the MCMC method, we apply the preconditioned Crank-Nicolson algorithm .

According to the pCN algorithm, the new obstacle parameter can be iteratively updated by the old parameter (initial guess) through the formula

 X=(1−β2)1/2Z+βω, (19)

where is the proposal variance coefficient and

is a zero-mean normal random vector with covariance matrix

. We remark that it is important and very tricky to select a suitable , because the value of dominates the proposal variance in the pCN algorithm. If is small, the parameter will be updated slightly in the MCMC sequence, leading to a time-consuming iteration process to get the ergodic in the space of obstacle parameters. If is big, the parameter may stay at one state for quite a long time with a huge number of iterations in the MCMC method. Consequently, one cannot get enough number of samples to approximate the posterior distribution, due to the computational cost prohibition. To over come this difficulty, we recommend the pCN algorithm with a random proposal variance  to obtain good MCMC sequences (see below for the description).

###### Algorithm

pCN Algorithm with Random Proposal Variance

• Initialize and .

• Repeat

1. Draw new obstacle parameter from the old state by the pCN algorithm (19) with the proposal variance coefficient as:

 (20)
2. Compute Hasting ratio as:

 α(X,Zj)=min{1,eΦ(Zj;Y)−Φ(X;Y)}; (21)
3. Accept or reject : draw and then update by the criterion

 Zj+1={X,if\; U≤α(X,Zj),Zj,if\, otherwise; (22)
4. Generate new proposal variance coefficient from . First we set

 (23)

with . In our case we choose . Then can be updated by

 βj+1=⎧⎪⎨⎪⎩βnew,ifβnew∈[0,1],−βnew,ifβnew<0,βnew−1,ifβnew>1. (24)

The randomness in the proposal variance has the potential advantage of including the possibility of large and small proposal variance coefficients . The large proposal variance coefficient helps the pCN algorithm explore the state space efficiently, and the small proposal variance coefficient protects the iterations from dropping into a fixed state. Then the random proposal variance gives rise to ergodic Markov chains.

## 3 Numerical Examples

In this section we exhibit numerical examples to demonstrate the effectiveness of the method described in the previous section. To save computational costs, we consider only two types of acoustically impenetrable scatterers in two dimensions:

• Sound-soft disks with unknown center and radius;

• Line cracks with unknown starting and ending points.

Obviously, there are totally three unknown parameters for disks and four parameters for cracks in 2D, implying that our parameters lie in a finite space with low dimensions.

With the definition of the observation (15), we construct two types of observations. In the first case, we consider an ideal model where the phaseless data are approximated by the exact forward solution without noise pollution. In the second case, we consider a practical model with noise-polluted far-field data.

### 3.1 Disk

In this subsection, we assume the underlying obstacle is a sound-soft disk. The inverse problem of recovering disks arises from, for instance, the polarization model of time-harmonic electromagnetic scattering from perfectly conducting cylinders whose cross-section is a disk. The parameterization of a disk is given by

 Z:=(z1,z2,z3)⊤=(x1,x2,r)⊤, (25)

where is the center and is the radius of the disk. Since , we assume that is a lognormal random variable.

Let the incident wave be given by the sum of two plane waves of the form (7). If the disk is located at the origin, it is well-known that the corresponding far-field pattern incited by the plane wave is given by the convergent series

 u∞(^x;Z,dℓ,k)=−e−iπ4√2πk⎡⎣J0(kr)H(1)0(kr)+2∞∑n=1Jn(kr)H(1)n(kr)cos(nθℓ)⎤⎦,ℓ=0,1,⋯,L. (26)

Here, denotes the angle between the observation direction and the incident direction , is the Bessel function of order and is the Hankel function of the first kind of order . If the disk is located at , by the translational formula (6) and the linear superposition principle, the exact far-field pattern with phase information of the scattered waves can be expressed as

 (27)

Note that the first line on the right hand side of (27) denotes the far-field pattern corresponding to the incoming plane wave , while the second line corresponding to .

The following theorem states that our phaseless data set is sufficient to uniquely identify a sound-soft disk. Let and be fixed. Then the data uniquely determine a sound-soft disk (that is, the center and radius of a disk). Suppose that () are two sound-soft disks centered at with the radius (). Set and denote by the far-field data corresponding to and the incident wave (7). Suppose that the phaseless far-field pattern are identical, i.e.,

 |u∞(^x;Z(1),d0,d,k)|=|u∞(^x;Z(2),d0,d,k)|for alld,^x∈S. (28)

In particular, choosing in the previous relation yields

 |u∞(^x;Z(1),d0,d0,k)|=|u∞(^x;Z(2),d0,d0,k)|=2|u∞(^x;Z(j),d0,k)|,j=1,2,

for all , where stands for the far-field pattern corresponding to the plane wave incident onto . This implies that

 |u∞(^x;Z(1),d0,k)|=|u∞(^x;Z(2),d0,k)|for all^x∈S.

Next, we shift the center of the disk to the origin and set . Recalling the translational formula (6), we obtain

 |u∞(^x;Z(1)0,d0,k)|=|u∞(^x;Z(2)0,d0,k)|% for all^x∈S.

Since the shifted disks with the parameters are rotationally invariant, the far-field pattern only depends on the angle between the incident direction and the observation direction . For any

, there exist a orthogonal matrix

such that . It then follows that (see e.g., [8, Chapter 5.1])

 u∞(^x;Z(j)0,d,k)=u∞(^x;Z(j)0,Qd0,k)=u∞(Q^x;Z(j)0,d0,k),∀^x∈S.

Combining the previous two identities yields

 |u∞(^x;Z(1)0,d,k)|=|u∞(^x;Z(2)0,d,k)|for all^x,d∈S,

which together with the translational formula implies

 |u∞(^x;Z(1),d,k)|=|u∞(^x;Z(2),d,k)|for all^x,d∈S. (29)

As a consequence of [40, Theorem 2.2], the relations (28) and (29) lead to the coincidence of and , which proves Theorem 3.1.

To apply the pCN algorithm (19) or (20), we need to set the key parameters and . For sound-soft disks, we set and let be the -by- identity matrix. Then the proposal is given by

 X=√0.99Zj+0.1ω,ω∼N(0,I). (30)

We describe the settings of our computational performance as follows:

• Unless otherwise specified, the wave number is always taken as ;

• The incident directions are ;

• The observation directions are ;

• To compute the far-field pattern (27), we truncate the infinite series of (27) by using the Bessel and first-kind Hankel functions of order ;

• In the setting of the prior distribution , we assume ;

• For the observation corresponding to the incident wave in (7), , we assume the observation pollution is a M-dimensional Gaussian variable, given by

 ηℓm=ση×|u∞(^xm;Z,d0,dℓ,k)|ωℓm, (31)

where , and is the noise coefficient. In our numerical tests we choose . In other words, we take and the diagonal matrix with ;

• Unless otherwise specified, the accurate obstacle parameters are set as , that is, a disk centered at with the radius .

In the first part, we adopt the ideal setting. Noting that in the formula (31) , we can gain a special sample of the observation noise with and . This would help us to investigate the accuracy of the numerical method and to verify the above uniqueness result in Theorem 3.1.

At first, we discuss the accuracy of the numerical solutions for different choice of and . Recall that the parameter denotes the number of incident waves and the number of observation directions. In Table 1

, we show the mean, standard deviation and the relative error of the reconstructed parameters with different choice of

. The numerical solutions of the recovered centers and radii are shown in Figures 1 and 2. The histograms of numerical solutions of the centers and radii are shown in Figures 4, 3 and 5, respectively.

Based on these results, we find that the reconstructed parameters are getting more accurate as the number of incident or observation directions become larger. The numerical solutions with , , , , are relatively inaccurate, since the resulted relative errors are large than in these cases. In contrast, the relative error is less if we choose and large enough such as , , . On the other hand, it can be observed from Table 1 that the standard deviation decreases as or increases. Table 2 and Figure 6 illustrate that small and may lead to unreliable reconstructions.

In the following we suppose that the location of the center is known and the knowledge of the radius needs to be recovered. Since only one parameter of the scatter remains unknown, we make use of minimal number of incident and observation directions by setting . In our tests we set incident directions , observation direction and accurate radius . The numerical approximations of radius vs different wave numbers are exhibited in Figure 7. For each fixed , we plot the phaseless far-field pattern against the radius in Figure 8.

From the numerical results we conclude that an accurate approximation of the radius can be obtained if the wave number is less than a threshold. It is seen from Figure 8 that the function is monotonically increasing in where as . This suggests that for large such as , there are more than one radii corresponding to the measured phaseless far-field pattern at . Hence, the reconstructed radii are inaccurate. These findings are consistent with the uniqueness result of , which states that a sound-soft disk can be uniquely determined from the phaseless far-field pattern at one observation direction, provided the radius is sufficiently small for a fixed wave number. The monotonicity property of the backscattered phaseless data with respect to the radius was rigorously justified in . Figure 7: Mean of the reconstructed radii with different wavenumbers k ranging from 10−10 to 1010 (top), from 10 to 100 (middle) and those from 90 to 100 (bottom). Figure 8: Phaseless far-field pattern |u∞(^x1)| vs radius r for different wavenumbers k.

Having verified the accuracy of our inversion scheme at the ideal setting with a special sample of the observation noise, we now consider the inverse problem with a general sample of the observation noise at the noise coefficient . In the second part, we estimate the obstacle parameters by setting , and . We generate one sample of the observation noise, which is a matrix with elements constructed by the formula (31). The numerical approximations from the polluted observation data are exhibited in Figures 9 and 10. The mean of the numerical center and radius are and , respectively. The standard deviations of these parameters are and the relative errors are . Figure 9: Reconstruction of the center (left) and radium (right) of a disk with one sample of observation noise at ση=3%. The red star ∗ (resp. line) is the accurate center (radius); the green dots ⋅ are the numerical reconstructions; the black circle ∘ (resp. line) is the mean of the centers and radii. We choose k=1, L=32 and M=64. Figure 10: Histograms of the reconstructed parameters x1,x2,r with one sample of the noisy data polluted at the level ση=3%.

To demonstrate the robustness of the numerical scheme, we generate 1000 samples of the observation noise. For each sample, one can calculate the mean solution of the parameters . Hence, we can perform statistical analysis over totally 1000 mean solutions of . In our tests, we pollute the phaseless data at different levels and exhibit the numerics in Table 3, Figure 11 and Figure 12. From these reconstructed parameters we conclude that the mean and relative error are robust against the noise pollution, but the standard deviation is very sensitive to the noisy level. Further, the phaseless data with less noise give rise to a more reliable reconstruction result. Figure 11: Reconstructions of the center (top) and radius (bottom) of a disk at different noise levels ση=3% (left), 6% (middle), 9% (right). The red star ∗ (resp. line) is the accurate center (resp. radius), the green dots ⋅ are the numerical reconstructions with each sample of observation noise and the black ∘ (resp. line) is the mean of reconstructed centers (resp. radii) with 1000 samples of observation noise. Figure 12: Histogram of 1000 numerical reconstructions of x1,x2,r with each sample of observation noise at different noise levels ση=3% (top), ση=6% (middle), ση=9% (bottom).

### 3.2 Line cracks

A crack or an open arc can be used to model the defects inside elastic and solid bodies such as bridge structures, aircraft engines and wings etc. Detection of such scatterers is important in safety and health assessment and is one of the fundamental topics in ultrasonic non-destructive testing. In this subsection, we want to recover a sound-soft crack of line-segment-type with the starting point at and the ending point . Hence, such line cracks can be characterized by parameters:

 Z:=(z1,z2,z3,z4)⊤=(x1,x2,y1,y2)⊤. (32)

Unlike the scattering from disks, we do not have an analytical expression of the far-field pattern corresponding to a line crack. Below we describe the integral equation method to solve the forward scattering problem, following the numerical scheme of  for general cracks. Denote by an open arc of class in 2D, which can be parameterized as

 Γ={z(s):s∈[−1,1]}. (33)

Using the integral equation method, the solution to the Helmholtz equation (2) in can be expressed as a single-layer potential ()

 usc(x)=∫ΓΦ(x,y)φ(y)ds(y),x∈R2∖Γ, (34)

where is the fundamental solution to the Helmholtz equation in two dimensions given by

 Φ(x,y):=i4H(1)0(k|x−y|),x,y∈R2,x≠y. (35)

Due to the Dirichlet boundary condition (3) on , the unknown density function is sought as a solution to the integral equation

 ∫ΓΦ(x,y)φ(y)ds(y)=f(x),x∈Γ,f=−uin. (36)

Once is calculated from (36), the far-field pattern could be expressed in the form

 u∞(^x)=eiπ/4√8πk∫Γe−ik^x⋅yφ(y)ds(y),^x∈S. (37)

To describe the numerical scheme of , we first introduce two functions defined on as follows

 H1(t,τ):={J0(k|z(cost)−z(cosτ)|)−1,t≠τ,0,t=τ, (38)

and

 H2(t,τ):=⎧⎪⎨⎪⎩πiH(1)0(k|z(cost)−z(cosτ)|)−{1+H1(t,τ)}ln(4e2[cost−cosτ]2),t≠τ,πi+2C+2ln{ke4|z′(cost)|},t=τ. (39)

Here, is the Euler’s constant. Then the integral equation (36) can be rephrased as

 (40)

for . Here,

 K(t,τ) = {1+sin2t−τ2K1(t,τ)}ln(4esin2t−τ2)+12H2(t,τ), (41) K1(t,τ) = ⎧⎪⎨⎪⎩H1(t,τ)sin2((t−τ)/2),t≠τ,−k2sin2(t)|z′(cost)|2,t=τ. (42)

The quadrature method  can be employed to discretize the integral equation (40

), based on the trigonometric interpolation with

equidistant nodal points , . Then the unknown solution to the integral equation (40) can be approximated by the discrete nodal values . Since with , it suffices to compute the discrete nodal values from the following algebraic system

 2n−1∑j=0ψj{R|k−j|+F|k−j|K1(tk,tj)+12nK2(tk,tj)}=g(tk),k=0,1,…,n, (43)

with

 Rj := 12n{c0+2n−1∑m=1cmcosmjπn+(−1)jcn},cm:=−1max{1,|m|}, (44) Fj := 12n{γ0+2n−1∑m=1γmcosmjπn+(−1)jγn},γm:=14(2cm−cm+1−cm−1). (45)

Note that there are totally unknown discrete nodal values