A meshless Galerkin scheme for the approximate solution of nonlinear logarithmic boundary integral equations utilizing radial basis functions

Accepted Manuscript A meshless Galerkin scheme for the approximate solution of nonlinear logarithmic boundary integral equations utilizing radial basis functions Pouria Assari, Mehdi Dehghan

PII: DOI: Reference:

S0377-0427(17)30580-0 https://doi.org/10.1016/j.cam.2017.11.020 CAM 11392

To appear in:

Journal of Computational and Applied Mathematics

Received date : 3 March 2017 Revised date : 29 July 2017 Please cite this article as: P. Assari, M. Dehghan, A meshless Galerkin scheme for the approximate solution of nonlinear logarithmic boundary integral equations utilizing radial basis functions, Journal of Computational and Applied Mathematics (2017), https://doi.org/10.1016/j.cam.2017.11.020 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Manuscript Click here to view linked References

A Meshless Galerkin Scheme for the Approximate Solution of Nonlinear Logarithmic Boundary Integral Equations Utilizing Radial Basis Functions Pouria Assari1,a , Mehdi Dehghanb a

Department of Mathematics, Faculty of Sciences, Bu-Ali Sina University, Hamedan 65178, Iran b Department of Applied Mathematics, Faculty of Mathematics and Computer Sciences, Amirkabir University of Technology, No. 424, Hafez Ave., Tehran 15914, Iran July 29, 2017

Abstract. This paper presents a computational scheme to solve nonlinear singular-logarithmic boundary integral equations. These types of integral equations result from boundary value problems of Laplace’s equations with nonlinear Robin boundary conditions. The discrete Galerkin method together with the (inverse) multiquadric radial basis functions established on scattered points are utilized to approximate the solution. The discrete Galerkin method for solving boundary integral equations results from the numerical integration of all integrals in the method. The proposed scheme utilizes a special accurate quadrature formula via the the nonuniform Gauss-Legendre integration rule to compute logarithm-like singular integrals appeared in the scheme. Since the numerical method developed in the current paper does not require any mesh generations on the boundary of the domain, it is meshless and does not depend to the domain form. We also investigate the error analysis of the proposed method. Illustrative examples show the reliability and efficiency of the new scheme and confirm the theoretical error estimates. MSC (2010): 45G05; 74S25; 65G99; 35J65 Keywords: Nonlinear boundary integral equation; Laplace’s equation; Logarithmic singular kernel; Discrete Galerkin method; Radial basis function; Error analysis.

1

Introduction

The radial basis functions (RBFs) are effective techniques for approximating an unknown function over a set of scattered points which have been used in the past few decades. RBFs include an independent variable without consideration of the problem dimensions, so applying them in higher dimensions does not raise the difficulties. Also it should be noted that the RBF method does not need any domain elements, so it is independent of the domain form. Firstly, Hardy has investigated RBFs as a multi-dimensional interpolation method, specially multiquadrics and inverse multiquadrics in 1971 to make a model of the earth’s gravitational field [33]. The thin plate splines as a type of free shape parameter RBFs have been developed by Meinguet [48] and have been utilized for smoothing multi-dimensional noisy data by Wahba [58]. A review paper [28] by Franke has compared scattered data approximations available in the early 1980. In recent years, the RBFs have been applied for the numerical solution of partial differential equations (PDEs) instead of scattered data interpolation. Kansa in 1990 [39, 40] introduced a meshless scheme for the numerical solution of PDEs based on the composition of collocation method and RBFs specially multiquadrics so-called Kansa’s method. Several types of PDEs have been solved 1

Corresponding author. E-mail addresses: [email protected] (P. Assari), [email protected] (M. Dehghan).

2 by Kansa’s idea, such as the nonlinear one-dimensional Burgers equation with shock wave [35], shallow water equations for tide and currents simulation [36], heat transfer equations [56, 61], parabolic equation with nonlocal boundary conditions [55], financial mathematic problems [37, 46], KdV equation [21], Klein-Gordon equation [23], Boussinesq equation [53] and sine-Gordon equation [22]. Afterwards, Kansa’s method has been expanded to Hermite-type interpolation method by Fasshauer [27] for the solvability of the collocation matrix. We would like to survey some newer numerical methods for solving integral equations utilizing meshless approximations. The meshless discrete collocation schemes based on the RBFs have been investigated for solving linear and nonlinear integral equations on non-rectangular domains with smooth kernels [1, 2] and weakly singular kernels [6, 7]. The RBFs have been applied to approximate the solution of one-dimensional Fredholm and Volterra integral equations [31] and nonlinear mixed Volterra-Fredholm integral equations [51]. The meshless product integration (MPI) method [5] has been proposed to solve one-dimensional linear weakly singular integral equations. The MLS methodology as a local meshless method has been used for solving linear and nonlinear integral equations on two-dimensional domains [4, 50] and integro-differential equations [20]. The paper [3] has described a computational method for solving Fredholm integral equations with logarithmic kernels so-called the meshless discrete Galerkin (MDG) method. Authors of [43, 42] have introduced a MLS-based meshless Galerkin method for numerically solving boundary integral equations and obtained the error analysis [44]. Consider the following nonlinear singular-logarithmic boundary integral equation: − πu(x) +

Z

u(y) ∂D

∂ ln kx − yk dsy + ∂ny

Z

∂D

g(y, u(y)) ln kx − ykdsy =

Z

∂D

f (y) ln kx − ykdsy , x ∈ ∂D, (1)

where D is an open, bounded, simply connected domain in R2 and its boundary is indicated by ∂D, nx denotes the outward unit normal vector on ∂D, the given function g(x, u) is assumed to be continuous on ∂D × R and f (x) is a known continuous function on ∂D and the function u(x) ∈ 2 (D) must be determined. These types of integral equations deduce from recomposition of ¯ C 1 (D)∩C the boundary value problem for two-dimensional Laplace’s equation with nonlinear Robin boundary conditions [9, 29] which arise in various branches of engineering and applied sciences such as heat transfer, electromagnetism, fluid and solid mechanics, scattering of waves, transmission line and so on [13, 14, 18, 49]. Since solving analytically boundary integral equations is mostly difficult, it is significant to give their numerical solutions. The Galerkin and collocation methods [10, 52] are useful schemes for approximating solutions of boundary integral equations which use a family of functions as basis. The piecewise polynomial collocation methods [8, 12] have been investigated for solving boundary integral equations. A numerical solution procedure [19] has been studied for solving two-dimensional Laplace’s equation based on the non-linear boundary conditions by mechanical quadrature methods. An iterative quadrature method [60] has been presented for the boundary integral equation deduced as reformulation of some two-dimensional Helmholtz equations. Legendre wavelets [24], spline wavelets [57], biorthogonal wavelets [32], trigonometric wavelet [16, 30] and Daubechies interval wavelets [54, 62] have been used to get the numerical solutions of boundary integral equations. A spectral Galerkin method [47] has been applied to solve boundary integral equations which arise from the two-dimensional Dirichlet problems. A useful research work conducted by authors of [45], has investigated a cell structure with logarithmical Gaussian quadrature schemes for solving boundary integral equations.

3 This article develops a computational scheme using RBFs to solve logarithmic singular boundary integral equations (1). The shape functions of (inverse) multiquadrics established on some distributed nodal points on ∂D are applied in the discrete Galerkin method. We will establish a special integration rule utilizing the composite non-uniform Gauss-Legendre quadrature formula to estimate singular integrals appeared in the new scheme. Because of the fact that the proposed method requires no mesh generation on the boundary of domain, it is a meshless method. This property results that the method is independent of the domain form. We also prove an error bound for the presented scheme. The implementation of method on computers is simple and also computationally attractive. The technique can be easily expanded to most classes of integral equations. The precision and convergence of the new approach are tested in various singular boundary integral equations. All results obtained in numerical examples verify the theoretical error estimates. In the following we consider the outline of the current paper. In Section 2, some basic formulations and error estimate of the multiquadrics and inverse multiquadrics are presented. In Section 3, we set up a new approximate method for solving the boundary integral equation (1) by combining the (inverse) multiquadric RBFs and the discrete Galerkin method. The error analysis of the proposed technique is obtained in Section 4. We consider some illustrative examples in Section 5. Finally, the conclusion of article is given in Section 6.

2

Multiquadrics RBFs

In this section, the (inverse) multiquadric RBFs as powerful tools for scattered data interpolation are introduced and utilized to approximate a function in any dimensions. Consider the radial function Φ : Rd → R defined as follows: Φ(x) = (kxk2 + c2 )β , c > 0, β 6= Z, where k.k is a common norm on Rd such as the Euclidean norm. The function Φ is called as multiquadrics when β = 12 and inverse multiquadrics when β = − 21 . It should be noted that the multiquadrics are positive definite functions and the inverse multiquadrics are conditionally positive definite functions of order one [59]. We require a set of scattered distinct nodal points such as {x1 , ..., xN } ⊂ Ω ⊂ Rd to estimate a function u(x) utilizing the (inverse) multiquadrics. To apply the RBF method for interpolating a function u(x) ∈ C(Ω), the following linear combination can be considered as [25, 59]: u(x) ≈ PN u(x) =

N X k=1

ck (kx − xk k2 + c2 )β =

N X k=1

ck Φ(x − xk ),

x ∈ Ω,

(2)

where the coefficients {ck }N k=1 are obtained by the interpolation conditions PN u(xi ) = u(xi ),

i = 1, ..., N.

(3)

In the following theorem, we discuss on the well-posedness of the interpolation problem (3). Theorem 2.1. [59] Suppose Φ is a (inverse) multiquadric RBF. Then for any distinct nodal points {x1 , ..., xN } ⊂ Ω ⊂ Rd , the corresponding interpolation matrix based on the expansion (2) is invertible. Therefore the interpolation problem (3) has a unique solution.

4 We now introduce another scheme to estimate a function u(x) ∈ L2 (Ω) via Galerkin’s idea. For general u(x) ∈ L2 (Ω), define QN to be solution of the following minimization problem [10]: ku − QN uk∞ =

min ku − zk∞ ,

(4)

z∈VN (Ω)

where VN (Ω) = span{Φ1 , ..., ΦN } is a subspace of L2 (Ω) with the dimension dN . Since VN (Ω) is finite dimensional, it can be shown that this problem has a solution; and by VN (Ω) being an inner product space, the solution can be shown to be unique [10] and is obtained by u(x) ≈ QN u(x) =

N X k=1

ck (kx − xk k2 + c2 )β =

N X k=1

ck Φ(x − xk ),

x ∈ Ω,

(5)

where {x1 , ..., xN } ⊂ Ω ⊂ Rd is a set of scattered distinct nodal points and Φ(x) = φ(kxk) is a (inverse) multiquadric RBF. In this situation, the coefficients {ck }N k=1 can be obtained by the following linear system: < u, φ(kx − xk k) >=

N X i=1

ci < φ(kx − xi k), φ(kx − xk k) >,

k = 1, ..., N,

where <, > is the inner product with respect to L2 (Ω), i.e. Z f (x)g(x)dx, f, g ∈ L2 (Ω). < f, g >=

(6)

(7)

Ω

Here, we want to represent the error estimate of (inverse) multiquadric method in the sense of the fill distance parameter [15, 59]. To establish this error analysis, we need to introduce some Hilbert spaces with respect to the (inverse) multiquadrics which are named native Hilbert spaces. Definition 2.1. [59] Suppose Φ is a (inverse) multiquadric RBF. Then the real native Hilbert space respect to Φ is defined as follows:

with inner product

fˆ ℵΦ (Ω) = {f ∈ C(Ω) ∩ L2 (Ω) : p ∈ L2 (Ω)}, ˆ Φ

Z ˆ g (w) 1 f (w)ˆ fˆ gˆ 1 p < f, g >ℵΦ (Ω) = √ < p , p >L2 (Ω) = √ dw, 2π ˆ ˆ 2π Ω ˆ Φ Φ Φ

(8)

(9)

where fˆ indicates Fourier transform of f .

From Definition 2.1, we concluded that the native spaces can be regarded as an extension of the standard Sobolev spaces. In other words, if the Fourier transform of Φ decays only algebraically, then the function Φ has a corresponding Sobolev space as its native space [25, 59]. Here some definitions from [15, 59] are represented which are significant to analyze the error of (inverse) multiquadrics method.

5 Definition 2.2. For a set of points X = {x1 , ..., xN } subset of a bounded domain Ω, the fill distance of X is described by hX,Ω = sup min kx − xj k2 . x∈Ω 0≤j≤N

Definition 2.3. The separation distance of X = {x1 , ..., xN } is given by qX =

1 min kxi − xj k2 . 2 i6=j

Also we call the set X as a quasi-uniform set corresponding to a constant c > 0 if qX ≤ hX,Ω ≤ cqX [25]. We are ready to represent the convergence theorem for approximating a function u ∈ ℵΦ (Ω) by a (inverse) multiquadric RBF from [59]. Before that we define the interior cone condition on a set in Rd as follows: Definition 2.4. [59] A set Ω ⊂ Rd confirms an interior cone condition if there exist a radius r > 0 and an angle θ ∈ (0, π/2) so that a unit vector ϑ(x) for every x ∈ Ω exists and the set Ω includes the cone C(x, ϑ(x), θ, r) = {x + µy : y ∈ Rd , kyk2 = 1, yT ϑ(x) ≥ cos θ, µ ∈ [0, r]}.

(10)

Theorem 2.2. [59] Let Φ be a (inverse) multiquadric RBF. Suppose Ω ⊂ Rd is a bounded and open set and also it satisfies an interior cone condition. We represent the approximation of a function u ∈ ℵΦ (Ω) using Φ and the distinct set X = {x1 , ..., xN } based on the collocation method by PN u. Then there exist constants h0 , C > 0 such that ( h−C )

ku − PN uk∞ ≤ e

X,Ω

|u|ℵΦ (Ω) ,

x ∈ Ω,

(11)

provided hX,Ω ≤ h0 . Remark 1. Based on the assumptions of Theorem 2.2, since PN u ∈ VN (Ω), it is clear that ku − QN uk∞ =

( h−C )

min ku − zk∞ ≤ e

X,Ω

z∈VN (Ω)

|u|ℵΦ (Ω) ,

u ∈ ℵΦ (Ω).

(12)

Therefore the convergence rates will be arbitrarily high for (inverse) multiquadrics [59].

3

Galerkin RBF Method

We consider the following Laplace’s equation ∆u(x) = 0,

x ∈ D ⊂ R2 ,

(13)

with the nonlinear Robin boundary condition ∂u(x) + g(x, u(x)) = f (x), ∂nx

x ∈ ∂D,

(14)

6 where D is an open, bounded, simply connected region in R2 , nx is the outward unit normal vector on ∂D, the known given function g(x, u) is assumed to be continuous on ∂D × R and f (x) is ¯ ∩ C 2 (D) must be a given continuous function on ∂D and the unknown function u(x) ∈ C 1 (D) determined [12]. By Green’s formula for harmonic functions, it is concluded that the solution u(x) is reformulated by [9, 12] Z Z 1 ∂u(y) ∂ ln kx − yk 1 u(x) = g(y, u(y)) dsy − ln kx − ykdsy , x ∈ D. (15) 2π ∂D ∂ny 2π ∂D ∂ny Applying the boundary condition of u(x) on ∂D for every x ∈ ∂D, the following logarithmic boundary integral equation is obtained Z Z Z ∂ ln kx − yk f (y) ln kx−ykdsy , x ∈ ∂D. u(y) g(y, u(y)) ln kx−ykdsy = −πu(x)+ dsy + ∂ny ∂D ∂D ∂D (16) Suppose {x1 , ..., xN } is N nodal points randomly selected on the boundary ∂D. We estimate the unknown function u(x) by the (inverse) multiquadrics as follows: u(x) ≈ u ¯N (x) =

N X j=1

c¯j φ(kx − xj k) =

N X

c¯j Φj (x),

j=1

x ∈ ∂D,

(17)

where Φj (x) = φ(kx − xj k) = (kx − xj k2 + c2 )β , are the shape functions of the (inverse) multiquadric Φ corresponding to the set X. We will find the coefficients {¯ cj , ..., c¯N } in the following. Assume L2 (∂D) is the square integrable function space on ∂D with respect to the inner product Z f (x)g(x)dsx , f, g ∈ L2 (∂D). (18) hf, gi = ∂D

By replacing the expansion (17) in the boundary integral equation (16) instead of u(x) and taking inner product h., Φi i upon both sides, we obtain −π

N X j=1

c¯j hΦj (x)Φi (x)i +

N X j=1

c¯j



  ∂ ln kx − yk , Φj (y) , Φi (x) ∂ny

+ hhg(y, u¯N (y)), ln kx − yki , Φi (x)i = hhf (y), ln kx − yki , Φi (x)i ,

(19)

or equivalently, we can rewrite −π

N X j=1

+

Z

c¯j

Z

N X

Z

c¯j Φj (x)Φi (x)dsx + ∂D j=1 {z } | | ∂D

| ∂D

Z

∂D

I1

Z

∂ ln kx − yk Φj (y)Φi (x) dsy dsx ∂ny ∂D {z } Z

I2

Φi (x)g(y, u¯N (y)) ln kx − ykdsy dsx = {z } | ∂D I3

Z

∂D

Φi (x)f (y) ln kx − ykdsy dsx . (20) {z } I4

7 The discrete Galerkin methods result from the numerical integration of all integrals in the system (20) associated with the Galerkin method. There are four forms of integrals to be estimated in this system. For approximating these integrals, assume that the boundary ∂D is a twice continuously differentiable, simple and closed curve in R2 [10]. Assume ∂D is parameterized by 0 ≤ t ≤ 1,

r(t) = (α(t), β(t)),

(21)

with r ∈ C 2 [0, 1] and kr′ (t)k = 6 0. Moreover, let the parametrization of ∂D traverse in a counterclockwise direction. We define the interior unit normal n(t) which is orthogonal to ∂D at r(t) as (−β ′ (t), α′ (t)) . (22) n(t) = p α′ (t)2 + β ′ (t)2

By the definition of r(t), we obtain

dsy = and dsx =

p

α′ (s)2 + β ′ (s)2 ds = kr′ (s)kds,

p

(23)

α′ (t)2 + β ′ (t)2 dt = kr′ (t)kdt.

(24)

Using these representations, we can rewrite the following integral as: Z 1 Φj (r(t))Φi (r(t))kr′ (t)kdt, I1 = hΦj (x)Φi (x)i =

(25)

0

and I2 =



  Z 1Z 1 ∂ ln kx − yk , Φj (y) , Φi (x) = q(t, s)Φj (r(s))Φi (r(t))kr′ (t)kkr′ (s)kdsdt, ∂ny 0 0

where q(t, s) =

 ′ ′ (s)[β(t)−β(s)]  √ ′−β2(s)[α(t)−α(s)]+α ,  ′ 2   α (s) +β (s) ([α(t)−α(s)]2 +[β(t)−β(s)]2 )    

′′

′′

′ (t)α (t)+α′ (t)β (t) √ −β , 2 α′ (s)2 +β ′ (s)2 (α′ (t)2 +β ′ (t)2 )

(26)

s 6= t, (27) s = t.

By applying the mN -point Gauss-Legendre integration formula corresponding to the coefficients {νk } and weights {wk } in the interval [−1, 1] to estimate the integrals (25) and (26), we have I1 = hΦj (x)Φi (x)i = where H1 [µk , q, i, j] = Φj and I2 =





M mN 1 XX 1 wk H1 [µk , q, i, j] + O( 2m ), 2M q=1 M N

(28)

k=1

   µk + 2q − 1 µk + 2q − 1 µk + 2q − 1 ) Φi r( ) kr′ ( )k, r( 2M 2M 2M

  mN M mN M X X ∂ ln kx − yk 1 1 XX wr H2 [µk , µr , q, i, j] + O( 2m ), , Φj (y) , Φi (x) = wk 2 ∂ny 4M M N q=1 k=1

p=1 r=1

(29)

8 where H2 [µk , µr , q, i, j] = q



µk + 2q − 1 µr + 2p − 1 , 2M 2M



Φj



   µk + 2q − 1 µr + 2p − 1 ) Φi r( ) r( 2M 2M µk + 2q − 1 µr + 2p − 1 ×kr′ ( )kkr′ ( )k. 2M 2M

Also, we can represent the dual singular integrals which appeared in the system (20) as follows: I3 = hhg(y, u¯N (y)), ln kx − yki , Φi (x)i =

Z

1

0

Z

1

0

Φi (r(t))g(r(s), u ¯N (r(s))) ln kr(t) − r(s)kkr′ (t)kkr′ (s)kdsdt, (30)

and I4 = hhf (y), ln kx − yki , Φi (x)i =

Z

0

1Z 1 0

Φi (r(t))f (r(s)) ln kr(t) − r(s)kkr′ (t)kkr′ (s)kdsdt.

(31)

Since these integrals are singular along the diagonal t = s, they cannot be computed by the classical numerical integration rule. Also the function ln kr(t) − r(s)k is periodic with period 1, so these integrals are also singular at points (0, 1) and (1, 0). Thus a particular quadrature formula is required. For approximating these integrals, the composite mN -point Gauss-Legendre rule with M non-uniform subdivisions is utilized [24]. Let h(t, s) be defined on [0, 1] × [0, 1] and satisfy 2m ∂ N h −ǫ−2mN , ∂t2mN < C1 t

2m ∂ N h −ǫ ∂s2mN < C2 t ,

(32)

for some ǫ ∈ (0, 1). Let α1 , α2 be functions in C 2mN [0, 1]. Then, for any given integer M , we have [24, 38] Z

0

1 Z α2 (t)

h(t, s)dtds =

α1 (t)

q=1

where θkq = with tq = ∆s(θkq )

mN M X ∆tq X

2

k=1

Mp,r mN ∆s(θkq ) X X 1 wp h(θkq , ηpr (θkq )) + O( 2m ), wk 2 M N

(33)

r=1 p=1

∆tq tq + tq−1 vk + t¯q , ∆tq = tq − tq−1 and t¯q = , 2 2

 q s , M

s=

α2 (θkq ) − α1 (θkq ) = Mp,r

2mN + 1 , 1−ǫ and

Mp,r = 1 + [M (α2 (θkq ) − α1 (θkq ))],

ηpr (θkq )

∆s(θkq ) 1 sp + α1 (θkq ) + (r − )∆s(θkq ). = 2 2

It is clear that the integrals (30) and (31) are not accessible for numerical computations by the quadrature formula (33), since the singularity takes place along the diagonal. To improve this problem, the following change of variables are most helpful [24, 38]: σ = t − s,

ν = s + t.

9 In fact, with the change of variables, the unit square [0, 1] × [0, 1] is transformed to the diamond {(σ, ν) : |σ| + |ν − 1| ≤ 1}, and so Z 1 h1 (σ)dσ, (34) I3 = hhg(y, u¯N (y)), ln kx − yki , Φi (x)i = −1

where h1 (σ) =

1 2

Z

α2 (σ)

Φi (r(

α1 (σ)

ν −σ ν −σ σ+ν ν −σ σ+ν ν −σ σ+ν ))g(r( ), u¯N (r( ))) ln kr( )−r( )kkr′ ( )kkr′ ( )kdν, 2 2 2 2 2 2 2

and I4 = hhf (y), ln kx − yki , Φi (x)i = where h2 (σ) =

1 2

Z

α2 (σ)

Φi (r( α1 (σ)

Z

1

h2 (σ)dσ,

(35)

−1

ν−σ σ+ν ν−σ σ+ν σ−ν σ+ν ))f (r( )) ln kr( ) − r( )kkr′ ( )kkr′ ( )kdν. 2 2 2 2 2 2

The integrals (34) and (35) have weakly singularity at u = 0, ±1 and are sufficiently smooth for every v. To approximate these singular integrals via the quadrature rule (33), we rewrite I3 = hhg(y, u¯N (y)), ln kx − yki , Φi (x)i

Z

=

−1 2

h1 (σ)du +

−1

Z

=

1

[h1 (

0

Z

0 −1 2

h1 (σ)du +

Z

1 2

h1 (σ)du + 0

Z

1 1 2

h1 (σ)dσ,

σ σ σ −σ ) + h1 ( ) + h1 ( − 1) + h1 (1 − )]dσ, 2 2 2 2

(36)

and similarly I4 = hhf (y), ln kx − yki , Φi (x)i =

Z

1

0

[h2 (

−σ σ σ σ ) + h2 ( ) + h2 ( − 1) + h2 (1 − )]dσ. 2 2 2 2

(37)

Therefore the integrands of (36) and (37) are singular in u = 0 and satisfy the condition (32) for any positive integers M and for any small positive numbers ǫ [24]. Now, using the quadrature rule (33), we compute I3 = hhg(y, u¯N (y)), ln kx − yki , Φi (x)i ≈

where

mN M X ∆σq X q=1

q q q θq ˜ 1 ( θk ) + ˜ ˜ 1 (1 − θk )], ˜ 1 ( −θk ) + h h1 ( k − 1) + h wk [h 2 2 2 2 } | {z 2 k=1 q H3 [θk ,ηpr ,i,¯ c1 ,...,¯ cN ]

(38)

Mp,r mN ηpr (θkq ) − σ ηpr (θkq ) − σ σ + ηpr (θkq ) ∆ν(θkq ) X X ˜ ))g(r( ), u¯N (r( ))) wp Φi (r( h1 (σ) = 4 2 2 2 r=1 p=1

ln kr(

σ + ηpr (θkq ) ηpr (θkq ) − σ σ + ηpr (θkq ) σ − ηpr (θkq ) ) − r( )kkr′ ( )kkr′ ( )k, (39) 2 2 2 2

with θkq =

 q s 2mN + 1 σq + σq−1 ∆σq , s= µk +¯ σq , ∆tq = tq −tq−1 , σ ¯q = , σq = , Mp,r = 1+[M (α2 (θkq )−α1 (θkq ))], 2 2 M 1−ǫ

10 α2 (θkq ) − α1 (θkq ) Mp,r

∆ν(θkq ) =

and ηpr (θkq ) =

∆ν(θkq ) 1 µp + α1 (θkq ) + (r − )∆ν(θkq ). 2 2

Also we have I4 = hhf (y), ln kx − yki , Φi (x)i ≈

mN M X ∆σq X

θkq θkq θkq −θkq ˜ ˜ ˜ ˜ ) + h2 ( ) + h2 ( − 1) + h2 (1 − )], wk [h2 ( 2 2 2 2 } | {z 2 k=1

q=1

where

H4 [θkq ,ηpr ,i]

(40)

Mp,r mN σ + ηpr (θkq ) ηpr (θkq ) − σ ∆ν(θkq ) X X ˜ wp Φi (r( ))f (r( )) h2 (σ) = 4 2 2 r=1 p=1

ln kr(

σ + ηpr (θkq ) ηpr (θkq ) − σ σ + ηpr (θkq ) σ − ηpr (θkq ) ) − r( )kkr′ ( )kkr′ ( )k. 2 2 2 2

(41)

Utilizing the numerical integration schemes (28), (29), (38) and (40) in the system (20) yields the linear system of algebraic equations −π

N X j=1

mN M mN M mN M X N X X 1 XX 1 XX cˆj wr H2 [µk , µr , q, i, j] cˆj wk wk H1 [µk , q, i, j] + 2M 4M 2 q=1 k=1

q=1 k=1

j=1

+

M X q=1

∆σq 2

mN X k=1

p=1 r=1

wk H3 [θkq , ηpr , i, cˆ1 , ..., cˆN ]

=

mN M X ∆σq X q=1

2

k=1

wk H4 [θkq , ηpr , i], (42)

for the unknowns Cˆ = [ˆ c1 , ..., cˆN ]. The solution of this system eventually leads to the following numerical solution as N X cˆj Φj (x), x ∈ ∂D. (43) u ¯N (x) = j=1

4

Error estimate

The operator K : C(∂D) → C(∂D) is introduced as  Z  ∂ ln kx − yk + g(y, u(y)) ln kx − yk dsy , Ku(x) = u(y) ∂ny ∂D

x ∈ ∂D,

(44)

by parameterizing ∂D by r(t) = (α(t), β(t)) for 0 ≤ t ≤ 1, we can rewrite the operator (44) as Z 1 {q(t, s)u(r(s)) + g(r(s), u(r(s))) ln kr(t) − r(s)k} kr′ (s)kds, 0 ≤ t ≤ 1. (45) Ku(x) = 0

Therefore, we can rewrite the boundary integral equation (16) in the operator form as

where g(x) =

Z

∂D

(−π + K)u = g,

(46)

Z

(47)

f (y) ln kx − ykdsy =

0

1

f (r(s)) ln kr(t) − r(s)k]kr′ (s)kds.

11 Also, the nonlinear integral operator F : C(∂D) → C(∂D) is introduced as follows: Fu(x) =

1 [Ku(x) − g(x)], π

x ∈ ∂D.

(48)

Therefore solving the integral equation (16) is equivalent to obtain the fixed points of operator F. Definition 4.1. A nonlinear operator on L2 (∂D) is recognized compact if every bounded set transforms into a precompact set in L2 (∂D), that’s mean an operator K : L2 (∂D) → L2 (∂D) is called compact if the set {Fu : kuk ≤ 1} has compact closure in L2 (∂D). Theorem 4.1. [41] The logarithmic integral operator (48) is compact on L2 (∂D). We define QN : L2 (∂D) → VN (∂D) as a Galerkin projection operator by QN u(x) =

N X

cj Φj (x) =

j=1

N X k=1

cj (kx − xj k2 + c2 )β ,

x ∈ ∂D,

(49)

where the space VN (∂D) = span{Φ1 , ..., ΦN } ⊂ L2 (∂D) with the dimension dN and the coefficients {c1 , ..., cN } obtained by solving the system < u, Φj >=

N X

cj < Φi , Φj >,

i = 1, ..., N,

(50)

j=1

where h., .i is the inner product on L2 (∂D) defined as Z 1 Z f (r(t))g(r(t))kr′ (t)kdt, f (x)g(x)dsx = hf, gi = 0

∂D

f, g ∈ L2 (∂D).

(51)

To obtain a better understanding of QN , we give an explicit formula for QN u. We introduce a new basis {Ψ1 , ..., ΨN } for VN (∂D) by using the Gram-Schmidt process to create an orthonormal basis from {Φ1 , ..., ΦN }. The element Ψi is a linear combination of {Φ1 , ..., ΦN } and moreover < Ψi , Ψj >= δij ,

i, j = 1, ..., N.

(52)

With this new basis, it is clear that QN u(x) =

N X

< u, Ψi > Ψi (x),

i=1

x ∈ ∂D.

(53)

Therefore the operator QN is called an orthogonal projection operator [10]. By using this operator, the system (19) is written in the operator form as QN F u ¯N = u ¯N .

(54)

Using the composite mN -point Gauss-Legendre integration rule over [−1, 1] with M subdivisions corresponding to the coefficients {µr } and weights {wk }, we define a discrete semidefinite inner product mN M X 1 X f (r(θkq ))g(r(θkq ))kr′ (θkq )k, f, g ∈ L2 (∂D), (55) wk hf, giN = 2M k=1

q=1

12 where θkq =

µk +2q−1 . 2M

Now we can introduce a discrete seminorm as √ kgkN = < g, g >N , g ∈ L2 (∂D).

(56)

As before, we know that < f, g >=< f, g >N +O(

1 M 2mN

f, g ∈ L2 (∂D).

),

We present the discrete projection operator as RN u(x) =

N X

x ∈ ∂D,

ck Φk (x),

k=1

(57)

where the coefficients {c1 , ..., cN } are found by solving the linear system hu, Φj iN =

N X k=1

ck hΦk , Φj iN ,

j = 1, ..., N.

(58)

Now, the following theorem is presented with reference to the discrete projection operator with the (inverse) multiquadric RBF shape functions as the basis. Lemma 4.1. [4] Suppose that the assumptions of Theorem 2.2 hold. Let RN , N ≥ 1 be the discrete orthogonal projections of the (inverse) multiquadric RBF relative to nodal points X = {x1 , ..., xN } ⊂ ∂D. Suppose {RN : N ≥ 1} is a uniformly bounded family on L2 (∂D), say kRN k ≤ m < ∞. If u ∈ ℵΦ (∂D) then RN u converges to u as N → ∞. Moreover, there exist constants h0 , C > 0 such that ( −C ) kRN u − uk∞ ≤ (1 + m)e hX,∂D |u|ℵΦ (∂D) , x ∈ ∂D, (59) provided hX,∂D ≤ h0 . Proof. Since PN u ∈ VN (∂D), from the properties of projection operators [10], we know that RN (PN u) = PN u. By the assumptions, we have sup kRN k∞ ≤ m.

(60)

N ≥1

Thus, we obtain u − RN u = u − PN u + PN u − RN u

(61)

= u − PN u + RN PN u − RN u = (I − RN )(u − PN u).

(62)

Therefore, kRN u − uk∞ ≤ (1 + m)ku − PN uk∞ .

(63)

Using Theorem 2.2, for every u ∈ ℵΦ (∂D), we have ( h −C )

kRN u − uk∞ ≤ (1 + m)e

X,∂D

|u|ℵΦ (∂D) ,

(64)

for all x ∈ ∂D, provided hX,∂D ≤ h0 . Since, hX,∂D → 0 as N → ∞ (justified by the quasi-uniform condition on X), it yields RN u → u.

13 Remark 2. The proof of uniformly bounded for the operators {QN } has been investigated at same length in Atkinson and Bogomolny [11]. Based on the use of the composite mN -point Gauss-Legendre rule with M subdivisions, a sequence of numerical integral operators FN on L2 (∂D) is given by mN M X 
 1X w ¯k u(r(θkq ))q(t, θkq ) + g(r(θkq ), r(θkq ))) − f (r(θkq ) ln kr(t) − r(θkq )k kr′ (θkq )k. FN u(x) = π k=1

q=1

(65) 1 We know that {FN } is a collectively compact family and pointwise convergent of O( M 2m ) to F N 2m 2 N (∂D), we have [24] on L (∂D) [24, 38]. Furthermore, it is shown that for every u ∈ C kFu − FN uk∞ ≤

CN sup |u(2mN ) (x)|, u ∈ C (2mN ) (∂D), M 2mN x∈∂D

(66)

where CN > 0 is a constant. Now, if we utilize the numerical integration operators (57) and (65) , then we can rewrite the system (42) as RN FN u ˆN = u ˆN . (67) Define the iterated solution by u ¯N = FN u ˆN ,

(68)

RN u ¯N = u ˆN ,

(69)

FN RN u ¯N = u ¯N .

(70)

then it is easily seen that and so To establish the error analysis of the new method, we present the convergence theorem of iterated Galerkin method from [11] as follows: Theorem 4.2. [11] Let the nonlinear boundary integral equation (16) have a unique solution u0 and assume that 1 is not an eigenvalue of F ′ (u0 ), where F ′ (u0 ) indicates the Frechet derivative of F at u0 . There are ε, N¯ > 0 such that the operator FN RN has a unique fixed point u ¯N in the ball ¯ . Moreover B(u0 , ε) for every N > N k¯ uN − u0 k∞ ≤ γ1 kFu0 − FN RN u0 k∞ ,

¯, N >N

(71)

where C > 0 is a constant. Now, we are ready to give the convergence theorem of the presented scheme. Theorem 4.3. Suppose that u0 ∈ ℵΦ (∂D) ∩ C 2mN (∂D) is the unique solution of the nonlinear boundary integral equation (16). Then based on the assumption of Theorem 2.2 and Lemma 4.1, for sufficiently large N the proposed method has a unique solution u ˆN on a ball around u0 . In addition ( h −C )

kˆ uN − u0 k∞ ≤ (1 + mγ1 γ2 )(1 + m)e where γ1 , γ2 , m, CN , C > 0 are constants.

X,∂D

|u|ℵΦ (∂D) +

CN mγ1 (2m ) sup |u0 N (x)|, 2m N M x∈∂D

(72)

14 Proof. This can be immediately obtained from Theorem 4.2 that there exists a unique iterated solution u ¯N ∈ B(u0 , ε) such that k¯ uN − u0 k∞ ≤ γ1 kFu0 − FN RN u0 k∞

≤ γ1 [kFu0 − FN u0 k∞ + kFN (u0 − RN u0 )k∞ ].

(73)

Since the family FN is the pointwise convergence, from the principle of uniform boundedness [10], it can be supposed that kFN k ≤ γ2 , so k¯ uN − u0 k∞ ≤ γ1 kFu0 − FN u0 k∞ + γ1 γ2 ku0 − RN u0 k∞ .

(74)

If u ˆ N = RN u ¯N then the presented method has the unique solution u ˆN . We consider the following decomposition u0 − u ˆN = u0 − RN u ¯N = (u0 − RN u0 ) + RN (u0 − u ¯N ). (75) Since kRN k ≤ m, by the inequality (74), we obtain kˆ uN − u0 k∞ ≤ ku0 − RN u0 k∞ + kRN kku0 − u ¯N k∞ .

≤ (1 + mγ1 γ2 )ku0 − RN u0 k∞ + mγ1 kFu0 − FN u0 k∞ .

(76)

1 ) and applying the error bound (59), With convergence order FN to F which is of order O( M 2m N we have ( h −C )

kˆ uN − u0 k∞ ≤ (1 + mγ1 γ2 )(1 + m)e

X,∂D

|u|ℵΦ (∂D) +

CN mγ1 (2m ) sup |u0 N (x)|. 2m N M x∈∂D

Since hX,∂D vanishes as N → ∞, justified by the quasi-uniform condition on X, kˆ uN − u0 k∞ → 0 ˆ > 0 so that for every N ≥ N ˆ , kˆ where N → ∞. Therefore, there is an integer N uN − u0 k∞ < εˆ. Now the proof is completed. Corollary 4.1. From Theorem 4.3, we find that the error of the proposed scheme investigated in the current paper appertains to the error of the mN -point numerical integration method and the (inverse) multiquadrics error. It is concluded that for sufficiently large integration nodes mN the error of the (inverse) multiquadrics is dominated over the error of Gauss-Legendre. Therefore increasing mN has no important effect on the global error and the approximation order of the proposed method will be arbitrarily high.

5

Numerical examples

In this section, we solve four nonlinear boundary integral equations on various kinds of boundaries to show accuracy and efficiency of the new method. Here we apply multiquadrics and inverse multiquadrics as basis in the scheme which are globally functions. This property eventuates a full and ill-conditioned coefficient matrix [17]. In this situation, the selection of the parameter c, called as the shape parameter, influences heavily on the accuracy of the method [59]. In fact, if we choose a bigger shape parameter with a fixed number N , the approximation will be more accurate while the condition number of associated coefficient matrixes will be large. Therefore, the selection of optimal value for c is still under the study and some values are suggested by many authors. For

15

Table 1: Some numerical results of Example 5.1 N 10 20 30 40 50 60

c

MQ keN k2 2.72 × 10−4 4.28 × 10−6 5.49 × 10−8 6.30 × 10−10 6.58 × 10−12 9.12 × 10−14

0.904 0.654 0.538 0.468 0.420 0.384

keN k∞ 4.02 × 10−4 9.30 × 10−6 1.23 × 10−7 1.48 × 10−9 1.61 × 10−11 2.79 × 10−13

Ratio − 5.43 10.65 15.38 20.25 22.24

Time 2.34 7.08 14.48 23.66 27.82 55.08

IM Q keN k2 2.54 × 10−4 3.97 × 10−6 5.68 × 10−8 5.95 × 10−10 1.09 × 10−12 9.24 × 10−14

keN k∞ 3.76 × 10−4 8.51 × 10−6 1.11 × 10−7 1.32 × 10−9 1.45 × 10−11 2.65 × 10−13

Ratio − 5.46 10.69 15.41 20.22 21.95

Time 2.32 7.25 14.52 23.71 36.31 55.80

P example, Hardy [33] proposed the use of c = 0.815d, where d = (1/N ) N i=1 di and di is the distance from xi to its nearest neighbor, for (inverse) multiquadric method in R2 . We choose an average value for c, i.e., c = √3N [5, 26]. We utilize 10-points composite Gauss-Legendre quadrature formula with non-uniform subdivisions M = 10 to compute singular logarithmic integrals in the scheme. We have measured the accuracy of presented technique by the maximum error keN k∞ and the mean error keN k2 which can be defined as follows: 1 Z 2 2 |uex (x) − u ˆN (x)| dx , keN k∞ = max {|uex (x) − u ˆN (x)|}, keN k2 = x∈∂D

∂D

where the exact solution uex (x) is estimated by the numerical solution u ˆN (x) presented in the current paper. The convergence rate of the presented scheme has been also reported by Ratio =

ln(keN k∞ ) − ln(keN ′ k∞ ) . ln(N ′ ) − ln(N )

We have written all routines in ”Maple” software with the ”Digits” 30 (Digits environment variable controls the number of digits in Maple) and a Laptop with 2.10 GHz of Core 2 CPU and 4 GB of RAM has been used to run these. To solve the final nonlinear system of algebraic equations the ”Fsolve” command has been employed based on the use of floating-point arithmetic. Example 5.1. Consider the following nonlinear boundary integral equation: −πu(x) +

where and

Z

u(y)

∂D

∂ ln kx − yk dsy + ∂ny

Z

∂D

25y13 + 25y22 + 101 ln kx − ykdsy = 1 + u2 (y)

Z

∂D

f (y) ln kx − ykdsy ,

 ∂D = (x1 , x2 ) : x21 + x22 − 2x2 = 1 ,

1.5(x3 − 8x21 ) + 0.2(2x22 − 12x2 ) + 25(x31 + x22 + 4), f (x) = f (x1 , x2 ) = p 1 3 2 (1 + 3(x2 − 6))(x1 + x2 + 4)

with the exact solution uex (x) = uex (x1 , x2 ) = √ 5

1 . x31 +x22 +4

We present the numerical results in terms of kek∞ , kek2 and the values of ratio at different numbers of N utilizing multiquadrics and inverse multiquadrics in Table 1. Also, the CPU times (sec.) for each RBFs are reported which show the short run-times needed for the new scheme. The spreading of nodes, selected randomly on the domain ∂D, is drawn in Figure 1. The obtained errors for different numbers of N are described in Figure 2. It seen that the obtained numerical results converge to the exact solution when the number of nodes increases with arbitrarily high order.

16

Figure 1: Node distribution for Example 5.1

17

Figure 2: Absolute error distributions of Example 5.1

Table 2: Some numerical results of Example 5.2 N 10 20 30 40 50 60

c 0.904 0.654 0.538 0.468 0.420 0.384

MQ keN k2 6.59 × 10−6 1.11 × 10−8 1.23 × 10−10 3.06 × 10−12 1.30 × 10−13 7.13 × 10−15

keN k∞ 1.88 × 10−5 3.24 × 10−8 3.73 × 10−10 9.52 × 10−12 4.06 × 10−13 9.80 × 10−15

Ratio − 9.09 11.14 12.75 14.14 20.42

Time 1.65 5.27 10.66 17.23 26.71 39.16

IM Q keN k2 1.33 × 10−5 3.50 × 10−8 4.72 × 10−10 1.34 × 10−11 6.27 × 10−13 3.85 × 10−14

keN k∞ 3.82 × 10−5 1.08 × 10−7 1.44 × 10−9 4.18 × 10−11 1.96 × 10−12 1.19 × 10−13

Ratio − 8.46 10.65 12.31 13.71 15.23

Time 1.59 5.22 10.80 17.56 26.79 39.01

18

Figure 3: Node distribution for Example 5.2

19

Figure 4: Absolute error distributions of Example 5.2

Example 5.2. Consider the following nonlinear boundary integral equation: −πu(x) +

Z

u(y)

∂D

∂ ln kx − yk dsy + ∂ny

where and

Z

∂D

(1 − u2 (y)) ln kx − yk p dsy = (2y1 y2 + y22 )2 + (y1 + 2y1 y2 )2

Z

∂D

f (y) ln kx − ykdsy ,

 ∂D = (x1 , x2 ) : (x1 − x2 )3 = x1 x2 ,

  cos(x1 + x2 + 1) x22 − x21 + cos(x1 + x2 + 1) p , f (x) = f (x1 , x2 ) = (2x1 x2 + x22 )2 + (x1 + 2x1 x2 )2

with the exact solution uex (x) = uex (x1 , x2 ) = sin(x1 + x2 + 1). Previous numerical methods have difficulties to solve these types of integral equations, but we can easily compute the approximate solution for this problem utilizing the meshless method presented in this work based on some random nodes over the ∂D which is shown in Figure 3. Table 2 reports results in terms of kek2 , kek∞ , the CPU times (sec.) and the values of ratio at different numbers of N utilizing multiquadrics and inverse multiquadrics. It is worth pointing out that the CPU times of the method are low which confirm that the new approach is fast. We have compared the obtained errors for different numbers of N in Figure 4. It is easy to see from Corollary 4.1 that increasing mN has no important effect on the global error and because for sufficiently large integration points mN , the error of the (inverse) multiquadrics is dominated over the error of Gauss-Legendre method. Example 5.3. Consider the following nonlinear boundary integral equation: −πu(x) +

where

Z

∂D

u(y)

∂ ln kx − yk dsy + ∂ny

Z

∂D

u(y) ln kx − ykdsy = 1 + u(y)

Z

∂D

  (x1 − 2)2 2 ∂D = (x1 , x2 ) : + (x2 − 6) = 1 , 4

f (y) ln kx − ykdsy ,

20

Figure 5: Node distribution for Example 5.3

21

Table 3: Some numerical results of Example 5.3 N 10 20 30 40 50 60

c 0.904 0.654 0.538 0.468 0.420 0.384

MQ keN k2 4.01 × 10−7 3.87 × 10−10 3.35 × 10−12 6.71 × 10−14 2.32 × 10−15 1.01 × 10−16

keN k∞ 9.69 × 10−7 1.19 × 10−9 1.01 × 10−11 2.07 × 10−13 7.22 × 10−15 2.77 × 10−16

Ratio − 9.66 11.75 13.53 15.03 17.88

Time 2.46 7.64 16.51 27.32 41.39 60.37

IM Q keN k2 3.69 × 10−7 5.86 × 10−10 5.57 × 10−12 1.04 × 10−13 2.24 × 10−15 1.14 × 10−16

keN k∞ 1.06 × 10−6 1.79 × 10−9 1.67 × 10−11 3.13 × 10−13 6.65 × 10−15 2.03 × 10−16

Ratio − 9.21 11.53 13.81 17.40 18.94

Time 2.51 8.03 16.97 28.09 42.31 61.12

Ratio − 7.14 10.86 19.33 29.80 24.33

Time 2.16 6.21 12.75 20.84 32.01 47.12

Figure 6: Absolute error distributions of Example 5.3

Table 4: Some numerical results of Example 5.4 N 10 20 30 40 50 60

c 0.904 0.654 0.538 0.468 0.420 0.384

MQ keN k2 3.65 × 10−4 3.19 × 10−6 2.59 × 10−8 1.04 × 10−11 1.89 × 10−14 6.85 × 10−16

keN k∞ 9.57 × 10−4 5.20 × 10−6 7.41 × 10−8 3.17 × 10−10 9.77 × 10−13 2.06 × 10−15

Ratio − 7.52 10.48 18.95 25.91 33.79

Time 2.08 6.03 12.55 20.71 31.86 46.81

IM Q keN k2 3.04 × 10−4 3.23 × 10−6 2.44 × 10−8 4.33 × 10−11 1.25 × 10−13 1.84 × 10−15

keN k∞ 7.68 × 10−4 5.43 × 10−6 7.13 × 10−8 2.90 × 10−10 2.75 × 10−13 4.44 × 10−15

22 and f (x) = f (x1 , x2 ) =

2x1 − 4x2 − 8x22 + 48x2 , (x21 + 3x2 + π)2

with the exact solution uex (x) = uex (x1 , x2 ) = x1 +2x1 2 +10 . Table 3 reports kek∞ , kek2 , the CPU times (sec.) and the values of ratio at different numbers of N utilizing multiquadrics and inverse multiquadrics. The dispersion of nodes, selected randomly on the domain ∂D, is drawn in Figure 5. We have compared the obtained errors for different numbers of N in Figure 6. Obviously, the method supplies the accurate solutions for the logarithmic singular boundary integral equation. Theorem 4.3 confirms that the results converge to the exact solution together with the increase of the scattered points of high order. Example 5.4. Consider the following boundary integral equation of the second-kind: −πu(x) +

Z

∂D

u(y)

∂ ln kx − yk dsy + ∂ny

where and

Z

∂D

cos(u(y)) p ln kx − ykdsy = 4y12 y2 − y12 + 1

 ∂D = (x1 , x2 ) : x41 + x22 − 4x2 + 2x21 x2 = 0 , 2

2

Z

∂D

f (y) ln kx − ykdsy ,

2

2

(4x41 + 4x21 x2 + 21 x21 + 12 x2 − 1)e2x1 +x2 −π + cos(e2x1 +x2 −π ) p f (x) = f (x1 , x2 ) = , 4x21 x2 − x21 + 1 2

2

with the exact solution uex (x) = uex (x1 , x2 ) = e2x1 +x2 −π . Numerical results are presented in Table 4 in terms of kek∞ , kek2 and the values of ratio at different numbers of N utilizing multiquadrics and inverse multiquadrics. Moreover, the used CPU times (sec.) have been reported for the proposed methods which show that the new approach is fast. The spreading of nodes, selected randomly on the domain ∂D, is drawn in Figure 7. We have compared the obtained errors for different numbers of N in Figure 8. All numerical results in this example have properly verified the theoretical error analysis obtained in Theorem 4.3.

6

Conclusion

In this paper, we have expanded a simple and numerically interesting method for solving nonlinear logarithmic boundary Fredholm integral equations using multiquadrics and inverse multiquadrics in the discrete Galerkin method. The logarithmic singular integrals occurred in this method have been computed utilizing a composite non-uniform Gauss-Legendre quadrature scheme. Since the new method has been constructed on a set of scattered points, it does not need any background mesh and so is independent of the geometry of the domain. We have also obtained the error bound for the scheme and found that the convergence rate of proposed method is arbitrarily high. As demonstrated by the computational results reported in Section 5, the proposed technique is able to produce high accurate solutions with low CPU time for various types of logarithmic boundary integral equations. All results in numerical examples have confirmed the theoretical error estimates.

Acknowledgments The authors are very grateful to the reviewers for their valuable comments and suggestions which have improved the paper.

23

Figure 7: Node distribution for Example 5.4

24

Figure 8: Absolute error distributions of Example 5.4

References [1] P. Assari, H. Adibi and M. Dehghan. A meshless method for solving nonlinear two-dimensional integral equations of the second kind on non-rectangular domains using radial basis functions with error analysis. J. Comput. Appl. Math., 239(1):72–92, 2013. [2] P. Assari, H. Adibi and M. Dehghan. A numerical method for solving linear integral equations of the second kind on the non-rectangular domains based on the meshless method. Appl. Math. Model., 37(22):9269–9294, 2013. [3] P. Assari, H. Adibi and M. Dehghan. A meshless discrete Galerkin (MDG) method for the numerical solution of integral equations with logarithmic kernels. J. Comput. Appl. Math., 267:160–181, 2014. [4] P. Assari, H. Adibi and M. Dehghan. A meshless method based on the moving least squares (MLS) approximation for the numerical solution of two-dimensional nonlinear integral equations of the second kind on non-rectangular domains. Numer. Algor., 67(2):423–455, 2014. [5] P. Assari, H. Adibi and M. Dehghan. The numerical solution of weakly singular integral equations based on the meshless product integration (MPI) method with error analysis. Appl. Numer. Math., 81:76–93, 2014. [6] P. Assari and M. Dehghan. A meshless method for the numerical solution of nonlinear weakly singular integral equations using radial basis functions. Eur. Phys. J. Plus., 132:1–23, 2017. [7] P. Assari and M. Dehghan. The numerical solution of two-dimensional logarithmic integral equations on normal domains using radial basis functions with polynomial precision. Eng. Comput., 2017. https://doi.org/10.1007/s00366-017-0502-5 [8] K. E. Atkinson. The numerical solution of a non-linear boundary integral equation on smooth surfaces. IMA J. Numer. Anal., 14(4):461–483, 1994. [9] K. E. Atkinson and G. Chandler. Boundary integral equation methods for solving Laplaces equation with nonlinear boundary conditions: The smooth boundary case. Math. Comput., 55(192):451–472, 1990. [10] K. E. Atkinson. The Numerical Solution of Integral Equations of the Second Kind. Cambridge University Press, Cambridge, 1997. [11] K. E. Atkinson and A. Bogomolny. The discrete Galerkin method for integral equations. Math. Comp., 48, 1987.

25 [12] K. E. Atkinson and D. Chien. Piecewise polynomial collocation for boundary integral equations. SIAM J. Sci. Comput., 16(3):651–681, 1995. [13] J. P. Bardhan. Numerical solution of boundary-integral equations for molecular electrostatics. J. Chem. Phys., 130(9), 2009. [14] Rokhlin V. Bremer, J. and I. Sammis. Universal quadratures for boundary integral equations on two-dimensional domains with corners. J. Comput. Physics, 229:8259–8280, 2010. [15] M. D. Buhmann. Radial Basis Functions: Theory and Implementations. Cambridge University Press, Cambridge, 2003. [16] W. Chen and W. Lin. Galerkin trigonometric wavelet methods for the natural boundary integral equations. Appl. Math. Comput., 121(1):75–92, 2001. [17] A. H. Cheng, M. A. Golberg, E. J. Kansa and Q. Zammito. Exponential convergence and h-c multiquadric collocation method for partial differential equations. Numer. Meth. Partial. Differ. Eqs., 19(5):571–594, 2003. [18] P. Cheng and J. Huang. Extrapolation algorithms for solving nonlinear boundary integral equations by mechanical quadrature methods. Numer. Algor., 58(4):545–554, 2011. [19] P. Cheng, J. Huang and Z. Wang. Mechanical quadrature methods and extrapolation for solving nonlinear boundary Helmholtz integral equations. Appl. Math. Mech., 32(12):1505–1514, 2011. [20] M. Dehghan and R. Salehi. The numerical solution of the non-linear integro-differential equations based on the meshless method. J. Comput. Appl. Math., 236(9):2367–2377, 2012. [21] M. Dehghan and A. Shokri. A numerical method for kdv equation using collocation and radial basis functions. Nonlinear Dyn., 50(1-2):111–120, 2007. [22] M. Dehghan and A. Shokri. A numerical method for solution of the two-dimensional sine-Gordon equation using the radial basis functions. Math. Comput. Simulat., 79(3):700–715, 2008. [23] M. Dehghan and A. Shokri. Numerical solution of the nonlinear Klein-Gordon equation using radial basis functions. J. Comput. Appl. Math., 230(2):400–410, 2009. [24] W. Fang, Y. Wang and Y. Xu. An implementation of fast wavelet Galerkin methods for integral equations of the second kind. J. Sci. Comput., 20(2):277–302, 2004. [25] G. E. Fasshauer. Meshfree methods. In Handbook of Theoretical and Computational Nanotechnology. American Scientific Publishers, 2005. [26] G. E. Fasshauer and J. G. Zhang. On choosing ”optimal” shape parameters for RBF approximation. Numer. Algor., 45(1-4):345–368, 2007. [27] G. E. Fasshauer. Solving partial differential equations by collocation with radial basis functions. In In: Surface Fitting and Multiresolution Methods A. Le M’ehaut’e, C. Rabut and L.L. Schumaker (eds.), Vanderbilt, pp. 131–138. University Press, 1997. [28] R. Franke. Scattered data interpolation: Tests of some methods. Math. Comput., 38(157):181–200, 1982. [29] M. Ganesh and O. Steinbach. Nonlinear boundary integral equations for harmonic problems. J. Integral Equ. Appl., 11(4):437–459, 1999. [30] J. Gao and Y. Jiang. Trigonometric Hermite wavelet approximation for the integral equations of second kind with weakly singular kernel. J. Comput. Appl. Math., 215(1):242–259, 2008. [31] A. Golbabai and S. Seifollahi. Numerical solution of the second kind integral equations using radial basis function networks. Appl. Math. Comput., 174(2):877–883, 2006. [32] H. Harbrecht and R. Schneider. Wavelet Galerkin schemes for boundary integral equations - implementation and quadrature. SIAM J. Sci. Comput., 27(4):1347–1370, 2006. [33] R. L. Hardy. Multiquadric equations of topography and other irregular surfaces. J. Geophys. Res., 176(8):1905– 1915, 2006. [34] J. Helsing and A. Karlsson. An explicit kernel-split panel-based Nystrom scheme for integral equations on axially symmetric surfaces. J. Comput. Phys., 272:686–703, 2014. [35] Y. C. Hon and X. Z. Mao. An efficient numerical scheme for Burgers’ equation. Appl. Math. Comput., 95(1):37– 50, 1998.

26 [36] Y.C. Hon, K.F. Cheung, X. Mao and E.J. Kansa. Multiquadric solution for shallow water equations. ASCE J. Hydraul. Eng., 125(5):524–533, 1999. [37] Y.C. Hon and X.Z. Mao. A radial basis function method for solving options pricing model. Financ. Eng., 8:31–49, 1999. [38] H. Kaneko and Y. Xu. Gauss-type quadratures for weakly singular integrals and their application to Fredholm integral equations of the second kind. Math. Comp., 62(206):739–753, 1994. [39] E. J. Kansa. Multiquadrics-a scattered data approximation scheme with applications to computational fluiddynamics-i surface approximations and partial derivative estimates. Comput. Math. Appl., 19(8-9):127–145, 1990. [40] E. J. Kansa. Multiquadrics-a scattered data approximation scheme with applications to computational fluiddynamics-ii solutions to parabolic, hyperbolic and elliptic partial differential equations. Comput. Math. Appl., 19(8-9):147–161, 1990. [41] B. Kress. Linear Integral Equations. Springer-Verlag, Berlin, 1989. [42] X. Li. Meshless Galerkin algorithms for boundary integral equations with moving least square approximations. Appl. Numer. Math., 61(12):1237–1256, 2011. [43] X. Li and J. Zhu. A Galerkin boundary node method and its convergence analysis. J. Comput. Appl. Math., 230(1):314–328, 2009. [44] X. Li and J. Zhu. A Galerkin boundary node method for biharmonic problems. Eng. Anal. Bound. Elem., 33(6):858–865, 2009. [45] X. Li and J. Zhu. A meshless Galerkin method for stokes problems using boundary integral equations. Comput. Methods Appl. Mech. Engrg., 198(37-40):2874–2885, 2009. [46] M. D. Marcozzi, S. Choi and C. S. Chen. On the use of boundary conditions for variational formulations arising in financial mathematics. Appl. Math. Comput., 124(2):197–214, 2001. [47] W. McLean. A spectral Galerkin method for a boundary integral equation. Math. Comput., 47(176):597–607, 1986. [48] J. Meinguet. Multivariate interpolation at arbitrary points made simple. Z. Angew. Math. Phys., 30(2):292–304, 1979. [49] D. Mirzaei and M. Dehghan. Implementation of meshless LBIE method to the 2d non-linear SG problem. Int. J. Numer. Methods Eng., 79(13):1662–1682, 2009. [50] D. Mirzaei and M. Dehghan. A meshless based method for solution of integral equations. Appl. Numer. Math., 60(3):245–262, 2010. [51] K. Parand and J. A. Rad. Numerical solution of nonlinear Volterra-Fredholm-Hammerstein integral equations via collocation method based on radial basis functions. Appl. Math. Comput., 218(9):5292–5309, 2012. [52] K. Ruotsalainen and J. Saranen. On the collocation method for a nonlinear boundary integral equation. J. Comput. Appl. Math., 28(C):339–348, 1989. [53] A. Shokri and M. Dehghan. A not-a-knot meshless method using radial basis functions and predictor-corrector scheme to the numerical solution of improved Boussinesq equation. Comput. Phys. Comm., 181(12):1990–2000, 2010. [54] X. Tang, Z. Pang, T. Zhu and J. Liu. Wavelet numerical solutions for weakly singular Fredholm integral equations of the second kind. Wuhan Univ. J. Nat. Sci., 12(3):437–441, 2007. [55] M. Tatari and M. Dehghan. On the solution of the non-local parabolic partial differential equations via radial basis functions. Appl. Math. Model., 33(3):1729–1738, 2009. [56] M. Tatari and M. Dehghan. A method for solving partial differential equations via radial basis functions: Application to the heat equation. Eng. Anal. Bound. Elem., 34(3):206–212, 2010. [57] T. Von Petersdorff and C. Schwab. Wavelet approximations for first kind boundary integral equations on polygons. Numer. Math., 74(4):479–519, 1996. [58] G. Wahba. Convergence rate of ”thin plate” smoothing splines when the data are noisy (preliminary report). Springer Lecture Notes in Math., 757, 1979.

27 [59] H. Wendland. Scattered Data Approximation. Cambridge University Press, New York, 2005. [60] Y. Yan. Fast numerical solution for a second kind boundary integral equation with a logarithmic kernel. SIAM J. Numer. Anal., 31(2):477–498, 1994. [61] M. Zerroukat, H. Power and C.S. Chen. A numerical method for heat transfer problems using collocation and radial basis functions. Internat. J. Numer. Methods Engrg., 42(7):1263–1278, 1998. [62] P. Zhang and Y. Zhang. Wavelet method for boundary integral equations. J. Comput. Math., 18(1):25–42, 2000.