- Email: [email protected]
Meshless methods for two-dimensional oscillatory Fredholm integral equations
Accepted Manuscript Meshless methods for two-dimensional oscillatory Fredholm integral equations Siraj-ul-Islam, Zaheer-ud-Din
PII: DOI: Reference:
S0377-0427(17)30581-2 https://doi.org/10.1016/j.cam.2017.11.021 CAM 11393
To appear in:
Journal of Computational and Applied Mathematics
Received date : 30 May 2017 Revised date : 26 October 2017 Please cite this article as: Siraj-ul-Islam , Zaheer-ud-Din, Meshless methods for two-dimensional oscillatory Fredholm integral equations, Journal of Computational and Applied Mathematics (2017), https://doi.org/10.1016/j.cam.2017.11.021 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.
Meshless methods for two-dimensional oscillatory Fredholm integral equations Siraj-ul-Islamb∗, Zaheer-ud-Dina,b a Department
of Basic Sciences, CECOS University of IT and Emerging Sciences, Peshawar, Pakistan.
b Department
of Basic Sciences, University of Engineering and Technology, Peshawar, Pakistan.
Abstract In this paper, a meshless solution procedure incorporating delaminating quadrature method for two-dimensional highly oscillatory Fredholm integral equation is put forward. The proposed method is an extension of our earlier findings of meshless methods for two-dimensional oscillatory Fredhom integral equation having kernel function free of stationary phase point(s) [1]. The current method not only deals with the kernels having no stationary phase point(s) but also with the kernels having stationary phase point(s) in the context of highly oscillatory integral equations. The new method is numerically stable and computationally fast. Numerical experiments are presented to testify the robustness and efficiency of the proposed method. Keywords: Fredholm integral equations, Levin’s quadrature, PDE, ODE, Multi-quadric radial basis function (MQ-RBF), Stationary phase point.
1
Introduction
Integral equation models have immense applications in the field of engineering and physics. Boundary value problems can be converted to Fredholm integral equations and vice versa. These mathematical models are used in the formulation of scattering problems and electromagnetic waves etc. [1–8]. One can write a two-dimensional Fredholm integral equation model of the second kind as [1]; x(t, y) = h(t, y) +
Z bZ a
d
k(t, y, u, v)x(u, v)dvdu,
c
(t, y) ∈ [a, b] × [c, d],
(1)
where x(t, y) is the unknown function. The oscillatory form of (1) will be called in this paper as highly oscillatory Fredholm integral equation (HOFIE) and can be expressed as: x(t, y) = h(t, y) +
Z bZ a
d
f (t, y, u, v)eι ω g(t,y,u,v) x(u, v)dvdu.
(2)
c
where ω is the frequency parameter, g(t, y, u, v), f (t, y, u, v), and u(t, y, u, v) are smooth functions. The function g is called the phase function. ∗
The author to whom all the correspondence should beaddressed. [email protected], [email protected]
1
Email:
[email protected],
The oscillatory integral equation models having kernel with stationary phase point(s) have applications in physics and engineering. Problems in acoustics and electromagnetic waves are modeled as the Helmholtz equation: ∆u(t) + ω 2 u(t) = 0, in Rd \Ω, d = 2 or 3, u = 0, on ∂Ω,
(3)
where Ω is the convex scattering obstacle. To characterize the scattering behaviour, one can reformulate to get the following form of oscillatory integral equation:[12-14] Z 1 ∂G(x, y) ∂uinc (x) ν(x) − + ιγG(x, y) ν(y)ds(y) = + ιγuinc (x), (4) 2 ∂n(x) ∂n(x) ∂Ω ∂u(x) where uinc (x) = eιωx is the incident plane wave, ν(x) = ∂n(x) is the surface current, n(x) is the unit normal vector at point x on ∂Ω, γ is a positive constant called coupling parameter and ι|x−y| G(x, y) = e|x−y| is the Green function. The kernel function K(x, y) = ∂G(x,y) ∂n(x) + ιγG(x, y) is then highly oscillatory. Accordingly, we can introduce an ansatz [9, 10]:
ν(x) = νslow (x)eιωx , to obtain 1 νslow (x) − 2
Z
∂Ω
K(x, y)eιω(y−x) νslow (y)ds(y) = ι(ω · n(x) + γ), x ∈ ∂Ω,
(5)
(6)
having a slow variational function νslow . The numerical solution of the integral equations (2) and (6) is more challenging, specially if its kernel is highly oscillatory having phase function with stationary phase point(s). The literature is deficient in handling such type of multi-dimensional models and the numerical procedure reported so far in the literature is specific to one-dimensional HOFIE. Finding efficient and accurate numerical solution of (6) is not easy due to the oscillatory kernel involved in the integral equation. To the best of our knowledge, no numerical method is so far reported in the literature for multidimensional HOFIE. The main contribution of the current work is to present an efficient and accurate numerical technique for two-dimensional HOFIE with and without stationary phase point(s). A quadrature method [10] designed for one dimensional HOFIE is based upon localized integration around stationary phase point(s) but this method is less efficient and prone to illconditioning. Some other numerical methods designed for such type of integral models are reported in [9, 11–20]. Levin’s type of methods for numerical solution of the oscillatory integrals and integral equations are reported in [17, 18, 21]. The nature of numerical solution of oscillatory integral equations is discussed in [22]. A new method is proposed in [21], for numerical solution of one-dimensional model HOFIE , which might be less accurate on nodal distribution other than Chebeshev-type. In [1], the authors have developed meshless methods based on Levin’s theory for univariate and multi-variate highly oscillatory integral equations. The current research is extension of the method [1] to highly oscillatory one and multi-dimensional integral equations having kernel functions containing stationary phase point(s). The proposed meshless algorithms is based on global and local interpolations. Both the methods is shown to perform well when the number of nodes is less. As the nodal points were increased, the global method was getting ill-conditioned due to the influence of the shape parameter. In contrast the local meshless method is not plagued 2
with the problem of ill-conditioning on dense points. It is also less sensitive to variation of the shape parameter. Some applications of local RBF in the context of numerical solution of PDEs are reported in [23–29]. The rest of the paper is summarized as follows. Techniques based on global and local MQ-RBF are given in section 2. Section 3 describes algorithm of the proposed methods handling the case of kernel function with and without stationary phase point(s). The robustness and effectiveness of the newly proposed algorithm is confirmed from the numerical experiments provided in section 5. Some conclusion are given in the section 6.
2
RBF-based Differentiation Matrices
Meshless procedures based on MQ-RBF namely, the global RBF-based differential matrix procedure (GDMP) and the local RBF-based differential matrix procedure (LDMP) are described in detail [1]. These procedure are currently used to solve differential equations given in (12) and c ), in a domain R2 , one can (23). Accordingly, for a set of N 2 centres (tc1 , y1c ), (tc2 , y2c ), ..., (tcN 2 , yN 2 construct a global dense differentiation matrix of order N 2 × N 2 w.r.t y, denoted by gDy , whose rows are the weights ws , each of order 1×N 2 corresponding to each centre (ts , ys ), s = 1, · · · , N 2 . This matrix can be written in the following form: w1 w2 gDy = . , (7) . . wN 2 where
ws = L (X) X−1 , L =
∂ , ∂y
and X is a global RBF matrix mentioned in [1]. Similarly, we can construct a one-dimensional RBF differentiation matrix w.r.t t, which is denoted by gDt . In the case of GDMP, a global dense matrix X needs to be inverted, which is a daunting task, specially in the case of multi-dimensional PDEs. The matrix X is often ill-conditioned which might lead to computational instability [1]. In order to avoid such complications, a localized meshless collocation procedure, the LDMP is adopted instead, for the correct and stable evaluation of the obtained ODEs or PDEs. In the LDMP, the derivatives of a function f are approximated at (ts , ys ), s = 1, 2, ..., N 2 , using a local stencil having size ns , satisfying {(ts1 , ys1 ), (ts2 , ys2 ), ..., (tsns , ysns )} ⊂ {(t1 , y1 ), (t2 , y2 ), ..., (tN 2 , yN 2 )}, for each ns N 2 , as explained in [1]. Consequently, one can construct a sparse local differentiation matrix lDy of order N 2 × N 2 , whose rows contain differential weights wsns , each of order 1 × N 2 , corresponding to the centres (ts , ys ), s = 1, · · · , N 2 accompanied by N 2 − ns zero entries, i.e w1ns w2n s lDy = .. . (8) . wNn2s 3
The vectors wsns are non-zero differential weights corresponding to the local stencils each having size ns . These weights are used to approximate function derivatives at sth nodal point and are represented by −1 L(Φsns ), wsns = Lsns
where Φ is the MQ-RBF. In similar lines one can construct lDt , a 1D RBF differentiation matrix w.r.t t as well.
3
New quadrature rule
To solve the integral equation given in (2), we start with discretizing it on the prescribed uniform or non-uniform nodes. The uniform nodes (tj , yj ) are generated as: tj = t1 + (j − 1)∆t, yj = y1 + (j − 1)∆y, j = 1, ..., M 2 ,
(9)
such that t = [t1 , t2 , · · · , tM 2 ]T and y = [y1 , y2 , · · · , yM 2 ]T are vectors each of order M 2 × 1. In case the of non-uniform nodes, we define T lir = lir1 , lir2 , · · · , lirM 2 . The vector lir can be achieved by using MATLAB rand command vrand = rand(1, M 2 ), in the following way: lir = s ◦ vrand ,
where s is either t or y, depending upon the direction and ◦ denotes the dot product of the vectors. The intensity of the randomness of the nodal points is shown in Fig. 3, which shows that the points are condensed around the centre of the domain.
3.1
Integral Equation in Discretized form
Discretized form of the integral equation (2) at the nodal points (tj , yj ) leads to; Z bZ d x(tj , yj ) = h(tj , yj ) + f (tj , yj , u, v)eι ω g(tj ,yj ,u,v) x(u, v)dvdu = h(tj , yj ) + Ij . a
where Ij =
(10)
c
Z bZ a
d
f (tj , yj , u, v)eι ω g(tj ,yj ,u,v) x(u, v)dvdu,
(11)
c
for each j = 1, ..., M 2 . Our goal is to compute the integral Ij efficiently and accurately. For this purpose, we choose the modified Levin’s quadrature procedure as given in [1].
3.2
Inner integral quadrature
We can calculate the inner integral about v in (11) by Levin’s [19] approach. Accordingly, a particular function p(u, v) satisfying the following PDE is obtained as: ∂g(tj , yj , u, v) ∂p(u, v) + ιω p(u, v) = f (tj , yj , u, v)x(u, v). ∂v ∂v Substituting (12) into (11), we get Z bZ d ∂g(tj , yj , u, v) ∂p(u, v) Ij = + ιω p(u, v) eι ω g(tj ,yj ,u,v) dvdu, j = 1, ..., M 2 . ∂v ∂v a c 4
(12)
Equivalently,
Z bZ
d
∂ (p(u, v)eι ω g(tj ,yj ,u,v) )dvdu, ∂v a c Z bh id = (p(u, v)eι ω g(tj ,yj ,u,v) ) du.
Ij =
(13) (14)
c
a
With the help of (12) and (14), we have Ij =
Z
b
ι ω g(tj ,yj ,u,d)
p(u, d)e
a
where I1 =
Z
du −
b
Z
b
a
p(u, c)eι ω g(tj ,yj ,u,c) du = I1 − I2 .
(15)
p(u, d)eι ω g(tj ,yj ,u,d) du,
(16)
p(u, c)eι ω g(tj ,yj ,u,c) du.
(17)
a
and I2 =
Z
b
a
To evaluate the integrals given in (16) and (17), we need to solve (12) to get the numerical value of p. The numerical procedure for evaluating p is given as follows. Let ur , r = 1, 2, ..., N 2 be a given u coordinate of the point ((ur , vs ), s = 1, 2, ..., N 2 ), then the numerical solution of the PDE given in (12) on every node on the v− axis ((ur , vs ), s = 1, 2, ..., N 2 ) constitutes a linear system of PDEs as follows:
∂p(ur ,v1 ) ∂v ∂p(ur ,v2 ) ∂v
.. .
∂p(ur ,vN 2 ) ∂v
+ ιω
∂g(tj ,yj ,ur ,v1 ) p(ur , v1 ) ∂v ∂g(tj ,yj ,ur ,v2 ) p(ur , v2 ) ∂v
.. .
∂g(tj ,yj ,ur ,vN 2 ) p(ur , vN 2 ) ∂v
=
f (tj , yj , ur , v1 )x(ur , v1 ) f (tj , yj , ur , v2 )x(ur , v2 ) .. . f (tj , yj , ur , vN 2 )x(ur , vN 2 )
.
(18)
In vector form, the above system can be expressed as:
∂gr ∂Pr +ιω Pr = Fr xr , ∂v ∂v where denotes the Hadamard product. Using the RBF differential matrix D given in section 2, one can write (19) can be simplified as: (Dv + ι ω Σr ) Pr = Fr xr , where T Pr = p(ur , v1 ), p(ur , v2 ), · · · , p(ur , vN 2 ) ,
T Fr = f (tj , yj , ur , v1 ), f (tj , yj , ur , v2 ), · · · , f (tj , yj , ur , vN 2 ) , Σr = diag and
h
∂g(tj ,yj ,ur ,v1 ) ∂g(tj ,yj ,ur ,v2 ) , , ∂v ∂v
··· ,
∂g(tj ,yj ,ur ,vN 2 ) ∂v
5
i
,
(19)
∂Pr ∂v
= Dv Pr . Hence (20)
xr = x(ur , vs ) is a vector, s = 1, 2, ..., N 2 , for each r =, 1, 2, ..., N 2 . Solution of (20) leads to N 2 different numerical vectors Pr (s) as Pr = Γr xr , r =, 1, 2, ..., N 2 , (21) where Γr = diag (Dv + ι ω Σr )−1 Fr . Using (21), we obtained the numerical solution of the PDE in (12) as; P1 P2 PN 2 ×N 2 = . . (22) . . PN 2
3.3
Outer integral quadrature
The numerical solutions obtained from (22) can be utilized to solve the one-dimensional integrals given in (15) at points c and d. The integral in (16) is evaluated by finding a proper function q(u) that satisfies the following ODE: ∂g(tj , yj , u, d) ∂q(u) + ιω q(u) = p(u, d). ∂u ∂u
(23)
Substituting the solution obtained from (23) into (16) we get I1 = q(b) eι ω g(tj ,yj ,b,d) − q(a) eι ω g(tj ,yj ,a,d) .
(24)
Incorporating Du of the GDMP or the LDMP (given in section 2) in (23) leads to (Du + ιωΣd ) qj = P1 ; where Σd = diag
h
∂g(tj ,yj ,u1 ,d) ∂g(tj ,yj ,u2 ,d) , , ∂u ∂u
··· ,
∂g(tj ,yj ,uN 2 ,d) ∂u
i
(25)
is a diagonal matrix.
To solve the outer integral in (24), we need to replace q(b) and q(a) by the first and the last components of the vector qj respectively and thus obtain I1 = Qdj qj = Qdj (Du + ιωΣd )−1 P1 = Qdj (Du + ιωΣd )−1 Γ1 xr , (26) where Qdj = eι ω g(tj ,yj ,b,d) , 0, · · · , 0, −eι ω g(tj ,yj ,a,d) . In similar lines, one can evaluate the integral I2 given in (17). Finally, the obtained numerical values of the integrals I1 and I2 are used in (15) to evaluate the integral Ij as:
or Assume
Ij = Qdj (Du + ιωΣd )−1 P1 − Qcj (Du + ιωΣc )−1 PN 2 ,
(27)
h i Ij = Qdj (Du + ιωΣd )−1 Γ1 − Qcj (Du + ιωΣc )−1 ΓN 2 xr .
(28)
Bj = Qdj (Du + ιωΣd )−1 Γ1 − Qcj (Du + ιωΣc )−1 ΓN 2 ,
(29)
6
we get the following form of (28), Ij = Bj xr , j = 1, ..., N 2 ,
(30)
where Bj are vectors of complex entries each of order 1 × N 2 , and xr is a complex vector of order N 2 × 1. Substituting the numerical value of the integral given in (30) into (10) we get; x(tj , yj ) = h(tj , yj ) + Ij , j = 1, ..., M 2 .
(31)
If the nodes in the (t, y) direction are specified to be similar to those in (u, v) direction (M 2 = N 2 and (tj , yj ) = (uj , vj ), j = 1, 2, 3, . . . , N 2 ) then (31) in matrix form can be written as: xn = Hn + Bn xn , n = 1, ..., N 2 , where
(32)
T xn = x(t1 , y1 ), x(t2 , y2 ), · · · , x(tN 2 , yN 2 ) ,
Hn is an N 2 × 1 vector of complex numbers and Bn is a matrix of order N 2 × N 2 having complex entries. Symbolically, these matrices can be expressed as; B1 h(t1 , y1 ) B2 h(t2 , y2 ) Bn = . , Hn = . .. .. . h(tN 2 , yN 2 )
BN 2
The total computation complexity of the algorithm given in (30) and (32) adds up to (N 2 + 1)(N 2 × N 2 ). In the local procedure the total computational cost of N 2 × N 2 global dense matrix reduces to N 2 × N 2 local sparse matrix which is solved by an efficient sparse matrix solver.
3.4
Integral equation with Stationary point(s)
The procedure when the kernel function is free of stationary point is straightforward as given in subsections 3.2 and 3.3. In the current section we introduce an algorithm of solving an oscillatory integral equation having kernel function with stationary phase point(s). The challenging aspect is to handle the existence of the stationary phase point(s) in such type of models. Oscillatory integrals involving stationary phase point(s) create difficulties for the meshless algorithm described in the previous section. Levin’s and modefied Levin’s type of methods fail at the stationary phase point(s) [9]. Such problem is addressed earlier in [1]. Accordingly, the subdivision of the domain interval can be done in two ways. In the first case, to divide the domain sub-interval at the position of the stationary phase point(s), the oscillatory integral having kernel function with stationary point say (ζu , ζv ) ∈ [a, b] × [c, d], looks like: Z bZ a
c
d
Υj dzdy =
Z
b
ζu
Z
d
ζv
Υj dvdu +
Z
b
ζu
Z
ζv
Υj dvdu +
c
Z
a
ζu
Z
d
ζv
Υj dvdu +
Z
a
ζu
Z
ζv
Υj dvdu, (33)
c
where Υj = f (tj , yj , u, v)eι ω g(tj ,yj ,u,v) x(u, v), for each j = 1, ..., N 2 . Using this mode of subdivision, one has to be careful about selecting number of nodal points. 7
In this case, (xj , tj ) 6= (ζu , ζv ). Thus, one should select dense points around stationary phase point(s) to obtain good accuracy. Therefore, (11) can be written in form of (33) as: Ij = Ij1 + Ij2 + Ij3 + Ij4 , where Ij1 Ij2
Z
=
ζu
Z
=
Z
Z
b
=
Υj dvdu,
(35)
Υj dvdu,
(36)
Υj dvdu,
(37)
Υj dvdu.
(38)
ζv
c
Z
ζu
a
Z
d
ζv
ζu
Ij3 = Ij4
Z
b
(34)
d
ζv
Z
ζu
a
ζv
c
In the second case, one can divide the domain sub-interval in such a way that the stationary phase point(s) is enclosed in small sub-interval separating it from the other sub-intervals where the kernel function is free of stationary phase point(s). The u− and v− components of the stationary-point(s) (ζu , ζv ) satisfies ζu ∈ (ζu − ηu , ζu + ηu ) and ζv ∈ (ζv − ηv , ζv + ηv ) where ηu , ηv > 0 are assumed to be relatively small numbers. This technique of sub-division can be preferred to use, as it leads to good accuracy irrespective of uniform or non-uniform nodal distribution. Consequently, the integral part of the integral equation (2) can be written as: Z bZ a
d
Υj dvdu =
c
+ + =
Z
b
Z
d
Υj dvdu +
ζu +ηu ζv +ηv Z ζu +ηu Z d ζu −ηu Z ζu −ηu a Ij1
ζv +ηv Z d
Υj dvdu + Υj dvdu +
Z
Z
b
+ Ij2 + · · · + Ij9 ,
where Ij1
=
ζu −ηu Z ζu −ηu
Ij3 Ij4 Ij5
=
Z
Z
b
Z
b
ζu +ηu
Z
b
ζu +ηu
= =
ζv −ηv Z ζv +ηv
Z
Υj dvdu + Υj dvdu +
ζv −ηv
Z
b
ζu +ηu c Z ζu +ηu ζu −ηu Z ζu −ηu a
ζv −ηv
Z
Υj dvdu
ζv −ηv
Υj dvdu
c
Z
ζv −ηv
Υj dvdu,
c
(39)
ζu +ηu
Ij2 =
Υj dvdu +
ζu +ηu ζv −ηv Z ζu +ηu Z ζv +ηv a
ζv +ηv
ζv +ηv
Z
Z
ζu +ηu
ζu −ηu
Υj dvdu,
(40)
Υj dvdu,
(41)
Υj dvdu,
(42)
Υj dvdu,
(43)
Υj dvdu,
(44)
ζv +ηv
Z
ζv +ηv
ζv −ηv
Z
ζu +ηu
ζu −ηu
d
ζv −ηv
c
Z
d
ζv +ηv
Z
ζv +ηv
ζv −ηv
8
Ij6
=
Z
ζu +ηu
ζu −ηu
Ij7
=
Z
Z
Z
Ij9 =
Z
and
Υj dvdu,
(45)
Υj dvdu,
(46)
Υj dvdu,
(47)
Υj dvdu.
(48)
c
Z
ζu −ηu
a
Ij8 =
ζv −ηv
d
ζv +ηv
ζu −ηu
Z
ζv +ηv
ζv −ηv
a ζu −ηu
a
Z
ζv −ηv
c
The integrals Ij1 − Ij9 shown above are evaluated on the local nodes in their respective subdomains in the following sequel: Dom Ij1 = [ζu + ηu , b] × [ζv + ηv , b], Dom Ij2 = [ζu + ηu , b] × [ζv − ηv , ζv + ηv ], Dom Ij3 = [ζu + ηu , b] × [c, ζv − ηv ], Dom Ij4 = [ζu − ηu , ζu + ηu ] × [ζv + ηv , d], Dom Ij5 = [ζu − ηu , ζu + ηu ] × [ζv − ηv , ζv + ηv ], Dom Ij6 = [ζu − ηu , ζu + ηu ] × [c, ζv − ηv ], Dom Ij7 = [a, ζu −ηu ]×[ζv +ηv , d], DomIj8 = [a, ζu −ηu ]×[ζv −ηv , ζv +ηv ], and DomIj9 = [a, ζu −ηu ]×[c, ζv −ηv ]. The nodes of these sub-domains are generally different from the global nodes occurring in the main interval Dom Ij = [a, b] × [c, d]. We can evaluate all these integrals using the procedure of delaminating quadrature mentioned in subsections 3.2 and 3.3, which transforms a two- dimensional problem into one-dimensional problem. Consequently, one can easily handle the case of stationary phase point(s) in a two-dimensional integral equation by using solution approach of a one-dimensional integral equation. To do this, as a priori requirement, we establish connectivity between the local and global nodes through barycentric interpolation [30, 31]. Let x(vj ) be the given function values for vj , j = 1, 2, ..., N 2 then the function value at inner point can be interpolated as [31]: x(v) = where the weights are defined as: $j = Q
PN 2
$j j=1 v−vj x(vj ) , PN 2 $j j=1 v−vj
1 , j, k = 1, 2, ..., N 2 . (v − v ) j k k6=j
For the uniform nodes, the weights are simplified as: $j = (−1)j
N 2! , − j)!j!
(N 2
whereas for Chebyshev-type nodes, the simplified weights are: 1/2, j = 1 or j = N 2 ; j $j = (−1) δj , δj = 1, otherwise. (49) can be expressed in matrix form as: T x(ζ) = l1 (ζ), l2 (ζ), · · · , lN 2 (ζ) x(v1 ), x(v2 ), · · · , x(vN 2 ) , 9
(49)
where coefficients are lj =
$j v−vj PN 2 $ j=1 v−vj
. r2 local nodes ξk , k = 1, 2, · · · , r2 can be expressed as: l2 (ξ1 ) · · · lN 2 (ξ1 ) x(v1 ) l2 (ξ2 ) · · · lN 2 (ξ2 ) x(v2 ) = M l X. .. .. .. .. . . . .
Consequently, the function values on l1 (ξ1 ) x(ξ1 ) x(ξ2 ) l1 (ξ2 ) = .. .. . . x(ξr2 ) l1 (ξr2 )
l2 (ξr2 ) · · · lN 2 (ξr2 )
x(vN 2 )
For each i = 1, 2, · · · , 9, and ti = 1, 2, · · · , Ni2 ; (uiti , vtii ) be the new nodes in the sub-intervals Dom Iji . Connectivity between the global nodes vk , k = 1, 2, · · · , N 2 and local ones are given in the v− direction as: P11 l1 (v11 ) l2 (v11 ) · · · lN 2 (v11 ) P1 P12 l1 (v21 ) l2 (v21 ) · · · lN 2 (v21 ) P2 P1N 2 ×N 2 = .. = .. .. .. .. = Ol1 PN 2 ×N 2 , .. . . . . . . 1 1 ) l (v 1 ) · · · l l1 (vN 2 N2 2 N 2 (vN 2 )
P1N 2 1
P2N 2 ×N 2 =
P21 P22 .. . P2N 2 2
and so on
P9N 2 ×N 2
=
1
P91 P92 .. . P9N 2 9
l1 (v12 ) l1 (v22 ) .. .
PN 2
1
1
l2 (v12 ) l2 (v22 ) .. .
··· ··· .. .
lN 2 (v12 ) lN 2 (v22 ) .. .
P1 P2 .. .
= 2 2 2 l1 (vN 2 ) l2 (vN 2 ) · · · lN 2 (vN 2 ) PN 2 2
2
2
lN 2 (v19 ) P1 9 l 2 (v2 ) P2 N .. = .. . . 9 9 2 9 l1 (vN 2 ) l2 (vN 2 ) · · · lN (vN 2 ) PN 2
l1 (v19 ) l1 (v29 ) .. . 9
l2 (v19 ) l2 (v29 ) .. .
··· ··· .. .
9
9
= Ol2 PN 2 ×N 2 .
= Ol9 PN 2 ×N 2 ,
where the superscript 1 − 9 indicates the solution on the respective sub-domains; DomIji . The connectivity between global nodes uk , k = 1, 2, · · · , N 2 and local nodes in the u− direction can be expressed as: l1 (u11 ) l2 (u11 ) · · · lN 2 (u11 ) x(u11 ) x(u1 ) x(u12 ) l1 (u12 ) l2 (u12 ) · · · lN 2 (u12 ) x(u2 ) Ψ1 = .. .. .. = Ml1 xr , = .. .. ... . . . . . l1 (u1N 2 ) l2 (u1N 2 ) · · · lN 2 (u1N 2 ) x(u1N 2 ) x(uN 2 ) 1
Ψ2 =
and so on
x(u21 ) x(u22 ) .. .
1
l1 (u21 ) l1 (u22 ) .. .
1
1
l2 (u21 ) l2 (u22 ) .. .
··· ··· .. .
lN 2 (u21 ) lN 2 (u22 ) .. .
x(u1 ) x(u2 ) .. .
= 2 2 2 2 x(uN 2 ) l1 (uN 2 ) l2 (uN 2 ) · · · lN 2 (uN 2 ) x(uN 2 ) 2
2
2
2
10
= Ml2 xr .
Ψ9 =
lN 2 (u91 ) x(u1 ) lN 2 (u92 ) x(u2 ) .. = = Ml9 xr . .. . . 2 (u9 ) l1 (u9N 2 ) l2 (u9N 2 ) · · · lN x(u ) x(u9N 2 ) 2 2 N N x(u91 ) x(u92 ) .. .
l1 (u91 ) l1 (u92 ) .. . 9
l2 (u91 ) l2 (u92 ) .. . 9
··· ··· .. .
9
After establishing connectivity between global and local nodes, the numerical results of the integrals Ij1 − Ij9 can be calculated in the similar way as shown in subsections 3.2 and 3.3 of the same section. The integrals Ij1 , j = 1, 2, ..., N 2 can be written as: Du + ιωΛ1u Ij1 = Q1u j j
where
−1
−1 Ol1 P1 − Q1l Du + ιωΛ1l Ol1 PN 2 , j j
−1 1u −1 l 1l 1l l = Q1u D + ιωΛ O Γ − Q D + ιωΛ O Γ Ml1 xr 2 u 1 u j j 1 j j 1 N
(50)
(51)
ι ω g(t ,u ,b,d) T j j Q1u , 0, · · · , 0, −eι ω g(tj ,uj ,ζu +ηu ,d) , j = e
ι ω g(t ,u ,b,ζ +η ) T v v , j j Q1l 0, · · · , 0, −eι ω g(tj ,uj ,ζu +ηu ,ζv +ηv ) , j = e h i ∂g(tj ,uj ,uN 2 ,d) ∂g(tj ,uj ,u1 ,d) ∂g(tj ,uj ,u2 ,d) Λ1u = diag , and , , · · · , j ∂u ∂u ∂u h i ∂g(tj ,uj ,uN 2 ,ζv +ηv ) ∂g(tj ,uj ,u1 ,ζv +ηv ) ∂g(tj ,uj ,u2 ,ζv +ηv ) Λ1l = diag . , , · · · , j ∂u ∂u ∂u
Similar numerical procedure can be adopted for other integrals; Ij2 − Ij9 to obtain the numerical value of Ij given in (39), over the entire domain as: Ij = + + + + + + + +
−1 Ol1 Γ1 − Q1l Du + ιωΛ1l Ol1 ΓN 2 Ml1 xr j j −1 −1 l Du + ιωΛ2u O2 Γ1 − Q2l Du + ιωΛ2l Ol2 ΓN 2 Ml2 xr j j j −1 −1 l l l 3l 3l Du + ιωΛ3u O Γ − Q D + ιωΛ O Γ M 2 u 1 j 3 N 3 xr j 3 j −1 −1 l 4l 4l l l Du + ιωΛ4u O Γ − Q D + ιωΛ O Γ M 2 u j 4 1 j j 4 N 4 xr −1 −1 l Du + ιωΛ5u O5 Γ1 − Q5l Du + ιωΛ5l Ol5 ΓN 2 Ml5 xr j j j −1 −1 l Du + ιωΛ6u O6 Γ1 − Q6l Du + ιωΛ6l Ol6 ΓN 2 Ml6 xr j j j −1 −1 l Du + ιωΛ7u O7 Γ1 − Q7l Du + ιωΛ7l Ol7 ΓN 2 Ml7 xr j j j −1 −1 l 8l 8l l l Du + ιωΛ7u O Γ − Q D + ιωΛ O Γ M 2 1 u j 8 j j 8 N 8 xr −1 −1 l 9l 9l l l Du + ιωΛ9u O Γ − Q D + ιωΛ O Γ M 2 u j 9 1 j j 9 N 9 xr .
Q1u Du + ιωΛ1u j j
Q2u j Q3u j Q4u j Q5u j Q6u j Q7u j Q8u j Q9u j
−1
11
(52)
Letting Uj =
−1 Ol1 Γ1 − Q1l Du + ιωΛ1l Ol1 ΓN 2 Ml1 j j −1 −1 l 2l 2l l Du + ιωΛ2u O Γ − Q D + ιωΛ O Γ Ml2 2 u j 2 1 j j 2 N −1 −1 l l 3l 3l Ml3 O Γ Du + ιωΛ3u O Γ − Q D + ιωΛ u 3 N2 j 3 1 j j −1 −1 l Du + ιωΛ4u O4 Γ1 − Q4l Du + ιωΛ4l Ol4 ΓN 2 Ml4 j j j −1 −1 l Du + ιωΛ5u O5 Γ1 − Q5l Du + ιωΛ5l Ol5 ΓN 2 Ml5 j j j −1 −1 l Du + ιωΛ6u O6 Γ1 − Q6l Du + ιωΛ6l Ol6 ΓN 2 Ml6 j j j −1 −1 l 7l 7l l Du + ιωΛ7u O Γ − Q D + ιωΛ O Γ Ml7 2 u j 7 1 j j 7 N −1 −1 l l 8l 8l Ml8 O Γ Du + ιωΛ7u O Γ − Q D + ιωΛ u 8 N2 j 8 1 j j −1 −1 l Du + ιωΛ9u O9 Γ1 − Q9l Du + ιωΛ9l Ol9 ΓN 2 Ml9 . j j j
Q1u Du + ιωΛ1u j j
+ Q2u j + Q3u j + Q4u j + Q5u j + Q6u j + Q7u j + Q8u j + Q9u j
−1
(53)
The interpolation matrices Mli , i = 1, . . . , 9 depend on the nature of stationary phase point(s). These matrices Ml1 − Ml9 remain unchanged for fixed stationary phase point(s). By substituting the value of Uj in (52) we have; Ij = Uj xr ,
(54)
where Uj are complex vectors each of order 1 × N 2 and xr is a vector of order N 2 × 1. Putting the value of Ij from (54) into (10) and following the steps given in subsection 3.3 the required solution of the integral equation having stationary-point is obtained.
4
Convergence Analysis
Lemma 4.1. (V ander corput, Stein [32], p.332). Suppose g(y) is real valued and smooth in (a, b) and that g k (y) ≥ 1 ∀ x ∈ (a, b) for fixed value of k then Z b −1 ι ω g(y) e dy ≤ c(k)ω k
(55)
a
holds when 1. k ≥ 2, or
2. k = 1 and g 0 (y) is monotonic. The bound c(k) is independent of g and ω, c(k) = 5.2k−1 − 2. Lemma 4.2. (Stein [32], p.334). Under the assumption on g(y) in Lemma 4.1, it has been concluded that Z b Z b 0 ι ω g(y) ≤ c(k)ω −1 ϕ (b) dy k e ϕ(y)dy |ϕ(b)| + (56) a
a
12
In the procedure of a delaminating quadrature a 2D integral is evaluated by transforming it into a 1D inner and outer integrals. Let a general 1D oscillatory integral mentioned in sections 3.2 and 3.3 is of the form Z b
Ix =
K(y)eι ω g(y) dy,
(57)
a
where K(y) = f (t, y)x(t, y), g(y) = g(t, y). The approximate value of (57) using the Levin’s quadrature theory and the procedures of GDMP and/or LDMP at descretized nodes tk , k = 1, 2, ..., N 2 is given by Z b K(y)eι ω g(y) dy, I n xn = a Z b ∂ = p(y)eι ω g(y) dy, a ∂y h ib = p(y)eι ω g(y) , y=a
= p(b)e
ι ω g(b)
− p(a)eι ω g(a) ,
(58)
where K(y) = K(tk , y) = f (tk , y)x(tk , y), g(y) = g(tk , y) for each fixed k. The error bound of the discretized form xn − I n xn = hn are discussed in the following theorems. Theorem 4.3. Suppose that f (y) and g(y) are suitably smooth functions and g 0 (y) 6= 0∀y ∈ [a, b]. Then (57) and (58) at arbitrary interpolation nodes y0 , y1 , · · · , yM satisfy
3(M + 2) F M +1 (y) ∞ (b − a)M +1 E = |Ix − I n xn | ≤ . (59) ω(M + 1)!
where F (y) = K(y) − L(y), L(y) ≡ p0 (y) + ιωg 0 (y)p(y). More specifically, when the shape parameter tends to zero, then (59) reduces to
3(M + 2)(b − a)M +1 K M +1 (y) ∞ E = |Ix − I n xn | ≤ . ω(M + 1)!
(60)
Proof. Let p(y) =
M X
λj Φ (r).
(61)
p 1 + 2 r2 ,
(62)
j=0
where Φ (r) is MQ-RBF defined as: Φ (r) =
and r2 = (y − yjc )2 , j = 0, 1, ..., M. According to Levin’s quadrature theory we have
L(y) ≡ p0 (y) + ιωg 0 (y)p(y), then L(y) ≡
M X j=0
λj Φ0 (r)
0
+ ιωg (y)
M X j=0
13
λj Φ (r), j = 0, 1, · · · , M,
(63)
(64)
where the interpolation nodes y0 , y1 , · · · , yM ∈ [a, b]. We have Z b ιωg(y) dy , E = |Ix − I n xn | = (K(y) − L(y))e a Z b ιωg(y) = F (y)e dy ,
(65) (66)
a
where F (y) = K(y) − L(y) satisfies the interpolation condition F (yj ) = 0, j = 0, 1, · · · , M . One can find dj ∈ (yj , yj+1 ) 3 F 0 (dj ) = 0, j = 0, 1, · · · , M − 1, then F (y) and F 0 (y) can be written in terms of Lagrange interpolation form as: M F M +1 (ξ1 ) Y F (y) = (y − yjc ), (M + 1)! j=0
and F 0 (y) =
M −1 F M +1 (ξ2 ) Y (y − yjc ), (M )! j=0
for some ξ1 , ξ2 ∈ [a, b]. By Lemma 4.1 and 4.2, for k = 1, c(k) = 3, we have
≤
Z
F 0 (y) dy, a
M +1 M +1
F 3(M + 2)(b − a) (y) ∞ . ω(M + 1)!
E = |Ix − I n xn | ≤ c(k)ω
−1/k
(|F (b)| +
b
(67)
When the shape parameter → 0, the interpolant L(y) converge to Lagrange interpolating polynomial [33, 34] i.e
LM +1 (y) ≡ 0 ⇒ F M +1 (y) ∞ = K M +1 (y) ∞ and above (67), then becomes
3(M + 2)(b − a)M +1 K M +1 (y) ∞ E = |Ix − I n xn | ≤ . ω(M + 1)!
(68)
Definition 1. (Kress [35],p.288). A linear operator I : X → Y from a normed space X into a normed space Y is called compact if for each bounded sequence (xn ) in X the sequence (Ixn ) contains a convergent subsequence in Y , i.e., if each sequence from the set {Ix : x ∈ X, kxk ≤ 1} contains a convergent subsequence. Definition 2. (Kress [36],p.5). A subset A of a normed space X is called complete if every Cauchy sequence of elements in A converges to an element in A. A normed space X is called a Banach space if it is complete.
14
Theorem 4.4. (Kress [35],p.292). Suppose that I : X → X be a compact linear operator on a Banach space X such that I − I is injective. Assume that the sequence I n : X → X of bounded
linear operators is norm convergent i.e I n − I → 0, n → ∞. Then for sufficiently large n the inverse operators (I − I n )−1 : X → X exists and are uniformly bounded. For solution of the equations x − Ix = f and xn − I n xn = fn we have an error estimate kxn − xk ≤ C{kI n x − Ixk + kfn − f k},
(69)
for some constant C. Theorem 4.5. Let I : [a, b] → [a, b] be a compact linear operator such that it satisfies the statement of the Theorem 4.4. The numerical approximation of the operator equation x−Ix = f is given by xn − I n xn = fn solved through meshless methods, either by the GDMP or the LDMP. Then the error bound at M + 1 collocation points satisfies (b − a)M +1 kxn − xk = O . (70) ω Proof. By Theorem 4.4, we have an error estimate kxn − xk ≤ C1 {kI n x − Ixk + kfn − f k},
(71)
for some constant C1 . Since f is the known function therefore both f and fn agree at the given interpolation points. (71) then becomes kxn − xk ≤ C1 {kI n x − Ixk ,
≤ C1 kI n − Ik kxk .
(72) (73)
Using the definition (1), for kxk ≤ 1 we have kxn − xk ≤ C1 kI n − Ik .
(74)
Incorporating the error bound of the integral operator given in Theorem 4.3, (74) simplified to
3(M + 2)(b − a)M +1 K M +1 (y) ∞ kxn − xk ≤ C1 , ω(M + 1)!
(b − a)M +1
K M +1 (y) , = C ∞ ω
1 (M +2) where C = 3C(M +1)! . By the uniform boundedness of K M +1 [11, 37], we obtained (b − a)M +1 kxn − xk = O . ω
This completes the proof.
15
(75) (76)
(77)
5
Numerical Experiments
In this section we have provided examples of oscillatory integral equations having kernel function with and without stationary phase point(s). The benchmark problems having kernels free of stationary phase point(s) are taken from [1]. Performance of the proposed methods is checked on both uniform and non-uniform points distribution. Furthermore, we have also checked efficiency of the proposed methods on non-uniform nodes. The non-uniform nodes are condensed around stationary-point(s) [1] which are useful to alleviate the negative effects of stationary phase point(s) on the accuracy of the method. Other nodal points distribution with uniform boundary points and scattered interior points is shown in Fig. 4. Similarly, Chebeshev-type of nodal points distribution is shown in Fig. 2. MQ-RBF depends on the value of the shape parameter. In general, we have taken this value to be = 0.5 in this paper. The accuracy is measured in terms of maximum relative error norm which is defined as follows: Ljab = |x(tj , yj ) − x ˜(tj , yj )|, Lre =
max |{z}
j=1,...,N 2
j = 1, · · · , N 2 ,
Ljab , |x(tj , yj )|
where the exact solution is represented by x and the computed solution by x ˜. Computations are executed using Core I5. The technique of the LDMP depends upon stencil size, which has been taken as 3 throughout the computation. Figure 1 is for Test Problem 1. Figures 2 and 3 are of Test Problem 2. Figures 4 and 5 are related to Test Problem 3. Figures 6-9 are related to Test Problem 4. Test Problem 1. Consider an oscillatory integral equation [1]; x(t, y) = h(t, y) +
Z
1
−1
Z
1
3 ev
−1
√ eι ω (3 (u + ecos( u + v)
7y 5
) − 5 t v) x(u, v) dvdu,
where h(t, y) = e
√ cos( t + y)
21
+
2e5
ιωy−5ιωt
sin(ω)(4 cos2 (ω) − 1) (− e− 1 + 10 ι ω t + e) . ω (− 1 + 5 ι ω t) √
This example has the exact solution x(t, y) = ecos( t+y) . The integral equation has an oscillatory kernel function free of stationary phase point(s). Fig. 1 (left) shows convergence in terms of relative error Lre of the proposed approach against increasing frequencies. The Lre is decreasing exponentially in the range 10−6 − 10−16 for fixed mesh-points N 2 = 16 and frequency range 104 − 109 . Fig. 1 (right) shows behavior of the Lre where N 2 varies in the range 16 − 324 for fixed frequency ω = 106 . The error norm Lre is not much influenced by increase in the nodal points and is maintained around O(10−11 ). It is witnessed that the performance of proposed methods: the GDMP and the LDMP is very good on a small number of points. Test Problem 2. Consider the integral equation [1]; x(t, y) = h(t, y) +
Z
1
−1
Z
1
−1
et + u ι ω (t u + v) e x(u, v) dvdu, cos(u) 16
−6
−9
10
10 LDMP GDMP
−8
10
−10
Lr e
Lr e
10
−10
10
−12
10
−14
10
−11
10
−16
10 10000
100000
1000000
10000000
100000000
1000000000
16
36
64
100
144
196
256
324
N2
ω
N 2 = 16
ω = 106
Figure 1: Test Problem 1, accuracy in terms of Lre versus ω, (left) and N 2 (right).
where
ι et −ι ω − ι ω t e−1 − e1 + 2ι ω t h(t, y) = cos(t) − ω (1 + ι ω t)
−1 + e2 ι ω
.
Exact solution of the integral equation is x(t, y) = cos(t). This problem is solved using Chebeshev-type of nodal points distribution as shown in Fig. 2. In Fig. 3 (left), accuracy of proposed algorithms the GDMP and the LDMP is shown against the increasing frequencies for fixed N 2 = 100. One can observe from this figure that as long as the frequencies are increasing in the range of 104 − 109 , accuracy of the proposed method increases in the range 10−4 − 10−14 . It is also witnessed from Fig. 3 (right) that the accuracy is retained with increasing nodal points. We have also observed that the proposed algorithms perform well in the case of Chebeshev-type nodes as well. Test Problem 3. Let an integral equation be of the form: Z 1Z 1 5 7v cos v x(t, y) = h(t, y) + eι ω (5 u − 3 − 2 ) x(u, v) dvdu, −1 −1 v + cos u where h(t, y) = y+cos t−2
7e
−5ιω 3
ω − 7e
10ιω 3
ω − 7ω cos(1)e
−31ιω 6
+ 7ω cos(1)e 5ω(49ω 2 − 4)
−ιω 6
− 2ι sin(1)e
−31ιω 6
+ 2ι sin(1)e
The exact solution of the problem is x(t, y) = y + cos t. We have considered this example to justify the results over another type of scattered nodal points as shown in Fig. 4 for N 2 = 100. Accuracy of the proposed algorithms GDMP and LDMP in terms of Lre is in the range of 10−7 −10−16 against increasing frequencies as shown in Fig 5 (left). This figure also confirms that whatsoever the arrangements of nodal points might be, the new proposed algorithms converge and hence are accurate. In Fig. 5 (right), we observed some mild fluctuations in the accuracy against increasing nodes for fixed frequency ω = 106 . 17
−ιω 6
;
1 0.8 0.6 0.4
y−axis
0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1
−0.8
−0.6
−0.4
−0.2
0 x−axis
0.2
0.4
0.6
0.8
1
N 2 = 324 Figure 2: Test Problem 2, model domain of a Chebeshev-type nodal points.
−4
−10
10
10 LDMP GDMP
−6
10
−8
Lr e
Lr e
10
−11
10
−10
10
−12
10
−12
10
−14
10 10000
100000
1000000
10000000
100000000
1000000000
16
36
64
100
144
196
256
324
N2
ω
N 2 = 100
ω = 106
Figure 3: Test Problem 2, accuracy in terms of Lre versus ω, (left) and N 2 (right) incorporating Chebyshev-type nodal points.
18
1 0.8 0.6 0.4
y−axis
0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1
−0.8
−0.6
−0.4
−0.2
0 x−axis
0.2
0.4
0.6
0.8
1
N 2 = 100 Figure 4: Test Problem 3, model domain of a non-uniform data.
−7
−10
10
10 LDMP GDMP
−8
10
−9
10
−10
10
−11
Lr e
Lr e
10
−11
10
−12
10
−13
10
−14
10
−15
10
−12
10
−16
10 10000
100000
1000000
10000000
100000000
16
36
64
100
144
196
256
324
N2
ω
N 2 = 100
ω = 106
Figure 5: Test Problem 3, accuracy in terms of Lre versus ω, (left) and N 2 (right) incorporating non-uniform data.
19
Test Problem 4. Let the Fredholm integral equation be Z 1Z 1 (sin t + e2v + y)eu+v ι ω (u2 + v2 ) x(t, y) = h(t, y) + e x(u, v) dvdu, cos u + e2u+v −1 −1 where ι √ √ √ ι ι ιω − 1 ιω + 1 −ιω π −2 sin t erf −ιω + e ω erf √ − 2y erf −ιω + e ω erf √ h(t, y) = x(t, y)+ erf 2ω −ιω −ιω 2t+y
The exact solution of the problem is x(t, y) = cos(t)+e . In this example the kernel function et+y contains the stationary phase point at 0. To handle such type of problem, we follow the procedure explained in section 3.4. According to the first strategy of domain sub-division, the oscillatory integral having kernel with stationary point 0 can be written as: Z 1Z 1 Z 1Z 1 Z 1Z 0 Z 0Z 1 Z 0Z 0 χ dvdu = χ dvdu + χ dvdu + χ dvdu + χ dvdu, −1
−1
0
0
(sin t+e2v +y)eu+v
0
2
−1
−1
−1
0
−1
2
where χ = cos u+e2u+v eι ω (u + v ) x(u, v). We have tested this mode of domain subdivision both on uniform and non-unform nodal points. Accuracy in terms of Lre increases with the increase of frequency on fixed non-uniform scattered nodal point(s) that are condensed around the stationary point 0 as shown in Fig. 6 (left). In Fig. 6 (right), we observed that for the same number of uniform nodal points, the Lre falls behind the scattered data points for increasing frequencies. In the second case, one can divide the domain interval in such a way that the stationary-point 0 ∈ (−η, η) (η > 0 is a small real number). Consequently, the integral in hand can be written as: Z 1Z 1 Z 1Z 1 Z 1Z η Z 1 Z −η Z ηZ 1 Z ηZ χ dvdu = χ dvdu + χ dvdu + χ dvdu + χ dvdu + −1
−1
+
η Z η
−η
η
Z
η
−η
−1
χ dvdu +
Z
−η −η Z 1
−1
η
χ dvdu +
Z
−1 −η
−1
η
Z
η
−η
−η
χ dvdu +
Z
η −η
−1
Z
−η
−η
χ dvdu.
−1
For this scheme of intervals sub-division, the proposed meshless procedure GDMP and LDMP produces accurate results as shown in Fig. 7. Accuracy in terms of Lre varies between 10−4 to 10−9 against increasing frequencies ω and fixed nodal points N 2 = 36, both for uniform and non-unform nodal points. Performance of the proposed meshless methods, GDMP and LDMP are equally good for uniform and non-uniform nodal points. Furthermore, we have observed from Fig. (9) that both the methods GDMP and LDMP are efficient on coarse grid as well.
6
Conclusion
The proposed solution algorithms are applied on different models of oscillatory integral equations. These models are implemented on uniform and non-uniform nodal points. The models given in test problems 1-3 are examples of Fredholm integral equations having oscillatory kernel function free of stationary phase point(s). The Fredholm integral equations model given in test problem 4 is an example where the oscillatory kernel function is having stationary phase point(s). The 20
η
−η
χ dvdu
−3
−2
10
10 LDMP GDMP
LDMP GDMP
−4
10
−3
10
−5
10
−4
Lr e
Lr e
10 −6
10
−5
10 −7
10
−6
10
−8
10
−9
10 10000
−7
100000
1000000
10000000
100000000
10 10000
1000000000
100000
1000000
10000000
ω
ω
N 2 = 16
N 2 = 16
100000000
1000000000
Figure 6: Test Problem 4, accuracy in terms of Lre versus ω, incorporating non-uniform data condensing around the stationary phase point (left) and uniform data (right) after subdivision of the interval at the position of stationary point.
−3
−3
10
10 LDMP GDMP
LDMP GDMP
−4
−4
10
10
−5
−5
10
Lr e
Lr e
10
−6
10
−7
−7
10
10
−8
−8
10
10
−9
10 10000
−6
10
−9
100000
1000000
10000000
100000000
1000000000
10 10000
100000
1000000
10000000
ω
ω
N 2 = 36
N 2 = 36
100000000
1000000000
Figure 7: Test Problem 4, accuracy in terms of Lre versus ω, incorporating non-uniform data condensing around the stationary point (left) and uniform data (right) after subdivision of the interval into sub-intervals inclosing stationary point.
21
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8 −1
−0.8
−0.6
−0.4
−0.2
N2
0
0.2
0.4
0.6
0.8
1
= 36
Figure 8: Test Problem 4, non-uniform data condensing around the stationary phase point model.
4
10
Uniform nodes−LDMP Scattered nodes−LDMP Uniform nodes−GDMP Scattered nodes−GDMP 3
time
10
2
10
1
10
36
81
144
N2
Figure 9: Test Problem 4, CPU time(s) versus N 2 .
22
proposed methods GDMP and LDMP give better accuracy for larger frequencies on a coarse grid for integral equation free of stationary phase point(s). In case when the kernel function incorporate stationary phase point(s), the domain interval is divided into sub-intervals according to the position of stationary phase point(s). The existence of stationary phase point(s) in the kernel of oscillatory integral equation offer difficulties in its solution. Levin’s type methods [19] fail in such scenario. To handle such type of issues in oscillatory integral equation models, we implemented the proposed methods on two types of domain sub-division. In the first type of domain sub-division, the proposed methods are converging fast when the nodes are non-uniform and condensing around the stationary phase point(s) but diverging after initial convergence when the nodes are uniform. In the second case of domain sub-division, the proposed methods work accurately both for uniform and non-uniform nodal points discretization. Thus, the proposed algorithms GDMP and LDMP are accurate in both the cases, with and without stationary phase point(s).
References [1] Siraj-ul-Islam, I. Aziz, Zaheer-ud-Din, Meshless methods for multivariate highly oscillatory Fredholm integral equations, Eng. Anal. Bound. Elem. 53 (2015) 100–112. [2] G. P. Galdi, K. Pileckas, A. L. Silvestre, On the unsteady poiseuille flow in a pipe., Z. Angew. Math. Phys. 58 (2007) 994–1007. [3] M. A. Bartoshevich, A heat-conduction problem., J. Eng. Phy. Thermophy. 28 (1975) 240– 244. [4] P. Baratella, A Nystrom interpolant for some weakly singular linear Volterra integral equations, Comput. Appl. Math. 231 (2009) 725–734. [5] A. M. Wazwaz, Linear and Non Linear Integral Equations: Methods and Appications, Higher Education Press, Beijing and Springer-Verlog Berlin Heidekberg, Germany, 2011. [6] K. Atkinson, The numerical solution of integral equations of the second kind, Cambridge University Press, Cambridge, 1997. [7] A. J. Jerri, Introduction to Integral Equations with Applications, John Wiley and Sons, INC, 1999. [8] H. Brunner, Y. Ma, Y. Xu, The oscillation of solutions of volterra integral and integrodifferential equations with highly oscillatory kernels, J. Integ. Eq. Appl. 27 (4) (2015) 455– 487. [9] C. Geuzaine, O. Bruno, F. Reitich, On the O(1) solution of multiple-scattering problems, IEEE Trans. Magn. 41 (5) (2005) 1488–1491. [10] O. Bruno, C. Geuzaine, J. Monro, F. Reitich, Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case, R. Soc. Lond. Trans. Ser. 362 (1816) (2004) 629–645. [11] Siraj-ul-Islam, S. Zaman, New quadrature rules for highly oscillatory integrals with stationary points, J. Comput. Appl. Math. 278 (2015) 75–89.
23
[12] H. Guoqiang, W. Jiong, Extrapolation of Nystrom solution for two-dimensional nonlinear Fredholm integral equations, J. Comput. Appl. Math. 134 (2001) 259–268. [13] W. J. Xie, F. R. Lin, A fast numerical solution of two-dimensional Fredholm integral equations of the second kind, Appl. Math. Comput. 59 (2009) 1709–1719. [14] A. Iserles, S. P. Norsett, Efficient quadrature of highly oscillatory integrals using derivatives, Proc. R. Soc. 461 (2005) 1388–1399. [15] D.Huybrechs, S.Vandewalle, A sparse discretisation for integral equation formulations of high frequency scattering problems, Tech. Rep., M.U.Leuven, 2006. [16] D. Huybrechs, S. Vandewalle, A sparse discretization for integral equation formulations of high frequency scattering problems, SIAM J. Sci. Comput. 29 (6) (2007) 2305–2328. [17] J. Li, X. Wang, T. Wang, A universal solution to one-dimensional highly oscillatory integrals, Sci. China. Ser. 51 (10) (2008) 1614–1622. [18] J. Li, X. Wang, T. Wang, C. Shen, Delaminating quadrature method for multi-dimensional highly oscillatory integrals, Appl. Math. Comput. 209 (2) (2009) 327–338. [19] D. Levin, Procedures for computing one and two-dimensional integrals of functions with rapid irregular oscillations, Math. Comput. 38 (1982) 531–538. [20] Siraj-ul-Islam, I. Aziz, W. Khan, Numerical integration of multi-dimensional highly oscillatory, gentle oscillatory and non-oscillatory integrands based on wavelets and radial basis functions., Eng. Anal. Bound. Elem. 36 (2012) 1284–1295. [21] J. Li, X. Wang, S. Xiao, T. Wang, A rapid solution of a kind of 1D Fredholm oscillatory integral equation, J. Comput. Appl. Math. 236 (2012) 2696–2705. [22] F. Ursell, Integral equations with a rapidly oscillating kernel, J. Lond. Math. Soc. 44 (1969) 449–459. [23] Q. Shen, Local rbf-based differential quadrature collocation method for the boundary layer problem, Eng. Anal. Bound. Elem. 34 (2010) 213–228. [24] S. Sarra, A local radial basis function method for advection-diffusion-reaction equations on complexly shaped domains, Appl. Math. Comput. 212 (2012) 9853–9865. [25] G. Liu, Meshfree Methods: Moving Beyond the Finite Element Method, CRC Press, 2003. [26] C. Shu, H. Ding, K. Yeo, Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations, Comp. Meth. Appl. Mech. Eng. 192 (2003) 941–954. ˇ [27] B. Sarler, R. Vertnik, Meshfree explicit local radial basis function collocation method for diffusion problems, Comput. Math. Appl. 51 (2006) 1269–1282. ˇ [28] G. Kosec, B. Sarler, Meshless approach to solving freezing with natural convection, Mater. Sci. Forum. 22 (2010) 205–210. ˇ [29] Siraj-ul-Islam, R. Vertnik, B. Sarler, Radial basis function collocation method for the numerical solution of the two-dimensional transient nonlinear coupled Burger’s equations, App. Math. Model. 36 (3) (2012) 1148–1160. 24
[30] D. Kincaid, W. Cheney, Numerical Analysis, Vol. 2, Math sci comp, Brooks/Cole, Pacific Grove, 2002. [31] J. Berrut, L. Trefethen, Barycentric lagrange interpolation, SIAM Review 46 (3) (2004) 501–517. [32] E. Stein, Harmonic Analysis: Real-rariable Methods, Orthogonality, and Oscillatory Integrals, Princeton University Press, Princeton, NJ. [33] T. Driscoll, B. Fornberg, Interpolation in the limit of increasingly flat radial basis function, Comput. Math. Appl. 43 (2000) 3–5. [34] G. Song, J. Riddle, G. E. Fasshauer, F. J. Hickernell, Multivariate interpolation with increasing flat radial basis functions of finite smoothness, Ad. Comput. Math. 36 (2012) 485–501. [35] R. Kress, Numerical Analysis, Springer-Verlag, New York, Inc., 1998. [36] R. Kress, Linear Integral Equations, Springer-Verlag, Berlin, 1989. [37] S. Xiang, Efficient quadrature for highly oscillatory integrals involving critical points, J. Comput. Appl. Math. 206 (2007) 688–698.
25
PII: DOI: Reference:
S0377-0427(17)30581-2 https://doi.org/10.1016/j.cam.2017.11.021 CAM 11393
To appear in:
Journal of Computational and Applied Mathematics
Received date : 30 May 2017 Revised date : 26 October 2017 Please cite this article as: Siraj-ul-Islam , Zaheer-ud-Din, Meshless methods for two-dimensional oscillatory Fredholm integral equations, Journal of Computational and Applied Mathematics (2017), https://doi.org/10.1016/j.cam.2017.11.021 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.
Meshless methods for two-dimensional oscillatory Fredholm integral equations Siraj-ul-Islamb∗, Zaheer-ud-Dina,b a Department
of Basic Sciences, CECOS University of IT and Emerging Sciences, Peshawar, Pakistan.
b Department
of Basic Sciences, University of Engineering and Technology, Peshawar, Pakistan.
Abstract In this paper, a meshless solution procedure incorporating delaminating quadrature method for two-dimensional highly oscillatory Fredholm integral equation is put forward. The proposed method is an extension of our earlier findings of meshless methods for two-dimensional oscillatory Fredhom integral equation having kernel function free of stationary phase point(s) [1]. The current method not only deals with the kernels having no stationary phase point(s) but also with the kernels having stationary phase point(s) in the context of highly oscillatory integral equations. The new method is numerically stable and computationally fast. Numerical experiments are presented to testify the robustness and efficiency of the proposed method. Keywords: Fredholm integral equations, Levin’s quadrature, PDE, ODE, Multi-quadric radial basis function (MQ-RBF), Stationary phase point.
1
Introduction
Integral equation models have immense applications in the field of engineering and physics. Boundary value problems can be converted to Fredholm integral equations and vice versa. These mathematical models are used in the formulation of scattering problems and electromagnetic waves etc. [1–8]. One can write a two-dimensional Fredholm integral equation model of the second kind as [1]; x(t, y) = h(t, y) +
Z bZ a
d
k(t, y, u, v)x(u, v)dvdu,
c
(t, y) ∈ [a, b] × [c, d],
(1)
where x(t, y) is the unknown function. The oscillatory form of (1) will be called in this paper as highly oscillatory Fredholm integral equation (HOFIE) and can be expressed as: x(t, y) = h(t, y) +
Z bZ a
d
f (t, y, u, v)eι ω g(t,y,u,v) x(u, v)dvdu.
(2)
c
where ω is the frequency parameter, g(t, y, u, v), f (t, y, u, v), and u(t, y, u, v) are smooth functions. The function g is called the phase function. ∗
The author to whom all the correspondence should beaddressed. [email protected], [email protected]
1
Email:
[email protected],
The oscillatory integral equation models having kernel with stationary phase point(s) have applications in physics and engineering. Problems in acoustics and electromagnetic waves are modeled as the Helmholtz equation: ∆u(t) + ω 2 u(t) = 0, in Rd \Ω, d = 2 or 3, u = 0, on ∂Ω,
(3)
where Ω is the convex scattering obstacle. To characterize the scattering behaviour, one can reformulate to get the following form of oscillatory integral equation:[12-14] Z 1 ∂G(x, y) ∂uinc (x) ν(x) − + ιγG(x, y) ν(y)ds(y) = + ιγuinc (x), (4) 2 ∂n(x) ∂n(x) ∂Ω ∂u(x) where uinc (x) = eιωx is the incident plane wave, ν(x) = ∂n(x) is the surface current, n(x) is the unit normal vector at point x on ∂Ω, γ is a positive constant called coupling parameter and ι|x−y| G(x, y) = e|x−y| is the Green function. The kernel function K(x, y) = ∂G(x,y) ∂n(x) + ιγG(x, y) is then highly oscillatory. Accordingly, we can introduce an ansatz [9, 10]:
ν(x) = νslow (x)eιωx , to obtain 1 νslow (x) − 2
Z
∂Ω
K(x, y)eιω(y−x) νslow (y)ds(y) = ι(ω · n(x) + γ), x ∈ ∂Ω,
(5)
(6)
having a slow variational function νslow . The numerical solution of the integral equations (2) and (6) is more challenging, specially if its kernel is highly oscillatory having phase function with stationary phase point(s). The literature is deficient in handling such type of multi-dimensional models and the numerical procedure reported so far in the literature is specific to one-dimensional HOFIE. Finding efficient and accurate numerical solution of (6) is not easy due to the oscillatory kernel involved in the integral equation. To the best of our knowledge, no numerical method is so far reported in the literature for multidimensional HOFIE. The main contribution of the current work is to present an efficient and accurate numerical technique for two-dimensional HOFIE with and without stationary phase point(s). A quadrature method [10] designed for one dimensional HOFIE is based upon localized integration around stationary phase point(s) but this method is less efficient and prone to illconditioning. Some other numerical methods designed for such type of integral models are reported in [9, 11–20]. Levin’s type of methods for numerical solution of the oscillatory integrals and integral equations are reported in [17, 18, 21]. The nature of numerical solution of oscillatory integral equations is discussed in [22]. A new method is proposed in [21], for numerical solution of one-dimensional model HOFIE , which might be less accurate on nodal distribution other than Chebeshev-type. In [1], the authors have developed meshless methods based on Levin’s theory for univariate and multi-variate highly oscillatory integral equations. The current research is extension of the method [1] to highly oscillatory one and multi-dimensional integral equations having kernel functions containing stationary phase point(s). The proposed meshless algorithms is based on global and local interpolations. Both the methods is shown to perform well when the number of nodes is less. As the nodal points were increased, the global method was getting ill-conditioned due to the influence of the shape parameter. In contrast the local meshless method is not plagued 2
with the problem of ill-conditioning on dense points. It is also less sensitive to variation of the shape parameter. Some applications of local RBF in the context of numerical solution of PDEs are reported in [23–29]. The rest of the paper is summarized as follows. Techniques based on global and local MQ-RBF are given in section 2. Section 3 describes algorithm of the proposed methods handling the case of kernel function with and without stationary phase point(s). The robustness and effectiveness of the newly proposed algorithm is confirmed from the numerical experiments provided in section 5. Some conclusion are given in the section 6.
2
RBF-based Differentiation Matrices
Meshless procedures based on MQ-RBF namely, the global RBF-based differential matrix procedure (GDMP) and the local RBF-based differential matrix procedure (LDMP) are described in detail [1]. These procedure are currently used to solve differential equations given in (12) and c ), in a domain R2 , one can (23). Accordingly, for a set of N 2 centres (tc1 , y1c ), (tc2 , y2c ), ..., (tcN 2 , yN 2 construct a global dense differentiation matrix of order N 2 × N 2 w.r.t y, denoted by gDy , whose rows are the weights ws , each of order 1×N 2 corresponding to each centre (ts , ys ), s = 1, · · · , N 2 . This matrix can be written in the following form: w1 w2 gDy = . , (7) . . wN 2 where
ws = L (X) X−1 , L =
∂ , ∂y
and X is a global RBF matrix mentioned in [1]. Similarly, we can construct a one-dimensional RBF differentiation matrix w.r.t t, which is denoted by gDt . In the case of GDMP, a global dense matrix X needs to be inverted, which is a daunting task, specially in the case of multi-dimensional PDEs. The matrix X is often ill-conditioned which might lead to computational instability [1]. In order to avoid such complications, a localized meshless collocation procedure, the LDMP is adopted instead, for the correct and stable evaluation of the obtained ODEs or PDEs. In the LDMP, the derivatives of a function f are approximated at (ts , ys ), s = 1, 2, ..., N 2 , using a local stencil having size ns , satisfying {(ts1 , ys1 ), (ts2 , ys2 ), ..., (tsns , ysns )} ⊂ {(t1 , y1 ), (t2 , y2 ), ..., (tN 2 , yN 2 )}, for each ns N 2 , as explained in [1]. Consequently, one can construct a sparse local differentiation matrix lDy of order N 2 × N 2 , whose rows contain differential weights wsns , each of order 1 × N 2 , corresponding to the centres (ts , ys ), s = 1, · · · , N 2 accompanied by N 2 − ns zero entries, i.e w1ns w2n s lDy = .. . (8) . wNn2s 3
The vectors wsns are non-zero differential weights corresponding to the local stencils each having size ns . These weights are used to approximate function derivatives at sth nodal point and are represented by −1 L(Φsns ), wsns = Lsns
where Φ is the MQ-RBF. In similar lines one can construct lDt , a 1D RBF differentiation matrix w.r.t t as well.
3
New quadrature rule
To solve the integral equation given in (2), we start with discretizing it on the prescribed uniform or non-uniform nodes. The uniform nodes (tj , yj ) are generated as: tj = t1 + (j − 1)∆t, yj = y1 + (j − 1)∆y, j = 1, ..., M 2 ,
(9)
such that t = [t1 , t2 , · · · , tM 2 ]T and y = [y1 , y2 , · · · , yM 2 ]T are vectors each of order M 2 × 1. In case the of non-uniform nodes, we define T lir = lir1 , lir2 , · · · , lirM 2 . The vector lir can be achieved by using MATLAB rand command vrand = rand(1, M 2 ), in the following way: lir = s ◦ vrand ,
where s is either t or y, depending upon the direction and ◦ denotes the dot product of the vectors. The intensity of the randomness of the nodal points is shown in Fig. 3, which shows that the points are condensed around the centre of the domain.
3.1
Integral Equation in Discretized form
Discretized form of the integral equation (2) at the nodal points (tj , yj ) leads to; Z bZ d x(tj , yj ) = h(tj , yj ) + f (tj , yj , u, v)eι ω g(tj ,yj ,u,v) x(u, v)dvdu = h(tj , yj ) + Ij . a
where Ij =
(10)
c
Z bZ a
d
f (tj , yj , u, v)eι ω g(tj ,yj ,u,v) x(u, v)dvdu,
(11)
c
for each j = 1, ..., M 2 . Our goal is to compute the integral Ij efficiently and accurately. For this purpose, we choose the modified Levin’s quadrature procedure as given in [1].
3.2
Inner integral quadrature
We can calculate the inner integral about v in (11) by Levin’s [19] approach. Accordingly, a particular function p(u, v) satisfying the following PDE is obtained as: ∂g(tj , yj , u, v) ∂p(u, v) + ιω p(u, v) = f (tj , yj , u, v)x(u, v). ∂v ∂v Substituting (12) into (11), we get Z bZ d ∂g(tj , yj , u, v) ∂p(u, v) Ij = + ιω p(u, v) eι ω g(tj ,yj ,u,v) dvdu, j = 1, ..., M 2 . ∂v ∂v a c 4
(12)
Equivalently,
Z bZ
d
∂ (p(u, v)eι ω g(tj ,yj ,u,v) )dvdu, ∂v a c Z bh id = (p(u, v)eι ω g(tj ,yj ,u,v) ) du.
Ij =
(13) (14)
c
a
With the help of (12) and (14), we have Ij =
Z
b
ι ω g(tj ,yj ,u,d)
p(u, d)e
a
where I1 =
Z
du −
b
Z
b
a
p(u, c)eι ω g(tj ,yj ,u,c) du = I1 − I2 .
(15)
p(u, d)eι ω g(tj ,yj ,u,d) du,
(16)
p(u, c)eι ω g(tj ,yj ,u,c) du.
(17)
a
and I2 =
Z
b
a
To evaluate the integrals given in (16) and (17), we need to solve (12) to get the numerical value of p. The numerical procedure for evaluating p is given as follows. Let ur , r = 1, 2, ..., N 2 be a given u coordinate of the point ((ur , vs ), s = 1, 2, ..., N 2 ), then the numerical solution of the PDE given in (12) on every node on the v− axis ((ur , vs ), s = 1, 2, ..., N 2 ) constitutes a linear system of PDEs as follows:
∂p(ur ,v1 ) ∂v ∂p(ur ,v2 ) ∂v
.. .
∂p(ur ,vN 2 ) ∂v
+ ιω
∂g(tj ,yj ,ur ,v1 ) p(ur , v1 ) ∂v ∂g(tj ,yj ,ur ,v2 ) p(ur , v2 ) ∂v
.. .
∂g(tj ,yj ,ur ,vN 2 ) p(ur , vN 2 ) ∂v
=
f (tj , yj , ur , v1 )x(ur , v1 ) f (tj , yj , ur , v2 )x(ur , v2 ) .. . f (tj , yj , ur , vN 2 )x(ur , vN 2 )
.
(18)
In vector form, the above system can be expressed as:
∂gr ∂Pr +ιω Pr = Fr xr , ∂v ∂v where denotes the Hadamard product. Using the RBF differential matrix D given in section 2, one can write (19) can be simplified as: (Dv + ι ω Σr ) Pr = Fr xr , where T Pr = p(ur , v1 ), p(ur , v2 ), · · · , p(ur , vN 2 ) ,
T Fr = f (tj , yj , ur , v1 ), f (tj , yj , ur , v2 ), · · · , f (tj , yj , ur , vN 2 ) , Σr = diag and
h
∂g(tj ,yj ,ur ,v1 ) ∂g(tj ,yj ,ur ,v2 ) , , ∂v ∂v
··· ,
∂g(tj ,yj ,ur ,vN 2 ) ∂v
5
i
,
(19)
∂Pr ∂v
= Dv Pr . Hence (20)
xr = x(ur , vs ) is a vector, s = 1, 2, ..., N 2 , for each r =, 1, 2, ..., N 2 . Solution of (20) leads to N 2 different numerical vectors Pr (s) as Pr = Γr xr , r =, 1, 2, ..., N 2 , (21) where Γr = diag (Dv + ι ω Σr )−1 Fr . Using (21), we obtained the numerical solution of the PDE in (12) as; P1 P2 PN 2 ×N 2 = . . (22) . . PN 2
3.3
Outer integral quadrature
The numerical solutions obtained from (22) can be utilized to solve the one-dimensional integrals given in (15) at points c and d. The integral in (16) is evaluated by finding a proper function q(u) that satisfies the following ODE: ∂g(tj , yj , u, d) ∂q(u) + ιω q(u) = p(u, d). ∂u ∂u
(23)
Substituting the solution obtained from (23) into (16) we get I1 = q(b) eι ω g(tj ,yj ,b,d) − q(a) eι ω g(tj ,yj ,a,d) .
(24)
Incorporating Du of the GDMP or the LDMP (given in section 2) in (23) leads to (Du + ιωΣd ) qj = P1 ; where Σd = diag
h
∂g(tj ,yj ,u1 ,d) ∂g(tj ,yj ,u2 ,d) , , ∂u ∂u
··· ,
∂g(tj ,yj ,uN 2 ,d) ∂u
i
(25)
is a diagonal matrix.
To solve the outer integral in (24), we need to replace q(b) and q(a) by the first and the last components of the vector qj respectively and thus obtain I1 = Qdj qj = Qdj (Du + ιωΣd )−1 P1 = Qdj (Du + ιωΣd )−1 Γ1 xr , (26) where Qdj = eι ω g(tj ,yj ,b,d) , 0, · · · , 0, −eι ω g(tj ,yj ,a,d) . In similar lines, one can evaluate the integral I2 given in (17). Finally, the obtained numerical values of the integrals I1 and I2 are used in (15) to evaluate the integral Ij as:
or Assume
Ij = Qdj (Du + ιωΣd )−1 P1 − Qcj (Du + ιωΣc )−1 PN 2 ,
(27)
h i Ij = Qdj (Du + ιωΣd )−1 Γ1 − Qcj (Du + ιωΣc )−1 ΓN 2 xr .
(28)
Bj = Qdj (Du + ιωΣd )−1 Γ1 − Qcj (Du + ιωΣc )−1 ΓN 2 ,
(29)
6
we get the following form of (28), Ij = Bj xr , j = 1, ..., N 2 ,
(30)
where Bj are vectors of complex entries each of order 1 × N 2 , and xr is a complex vector of order N 2 × 1. Substituting the numerical value of the integral given in (30) into (10) we get; x(tj , yj ) = h(tj , yj ) + Ij , j = 1, ..., M 2 .
(31)
If the nodes in the (t, y) direction are specified to be similar to those in (u, v) direction (M 2 = N 2 and (tj , yj ) = (uj , vj ), j = 1, 2, 3, . . . , N 2 ) then (31) in matrix form can be written as: xn = Hn + Bn xn , n = 1, ..., N 2 , where
(32)
T xn = x(t1 , y1 ), x(t2 , y2 ), · · · , x(tN 2 , yN 2 ) ,
Hn is an N 2 × 1 vector of complex numbers and Bn is a matrix of order N 2 × N 2 having complex entries. Symbolically, these matrices can be expressed as; B1 h(t1 , y1 ) B2 h(t2 , y2 ) Bn = . , Hn = . .. .. . h(tN 2 , yN 2 )
BN 2
The total computation complexity of the algorithm given in (30) and (32) adds up to (N 2 + 1)(N 2 × N 2 ). In the local procedure the total computational cost of N 2 × N 2 global dense matrix reduces to N 2 × N 2 local sparse matrix which is solved by an efficient sparse matrix solver.
3.4
Integral equation with Stationary point(s)
The procedure when the kernel function is free of stationary point is straightforward as given in subsections 3.2 and 3.3. In the current section we introduce an algorithm of solving an oscillatory integral equation having kernel function with stationary phase point(s). The challenging aspect is to handle the existence of the stationary phase point(s) in such type of models. Oscillatory integrals involving stationary phase point(s) create difficulties for the meshless algorithm described in the previous section. Levin’s and modefied Levin’s type of methods fail at the stationary phase point(s) [9]. Such problem is addressed earlier in [1]. Accordingly, the subdivision of the domain interval can be done in two ways. In the first case, to divide the domain sub-interval at the position of the stationary phase point(s), the oscillatory integral having kernel function with stationary point say (ζu , ζv ) ∈ [a, b] × [c, d], looks like: Z bZ a
c
d
Υj dzdy =
Z
b
ζu
Z
d
ζv
Υj dvdu +
Z
b
ζu
Z
ζv
Υj dvdu +
c
Z
a
ζu
Z
d
ζv
Υj dvdu +
Z
a
ζu
Z
ζv
Υj dvdu, (33)
c
where Υj = f (tj , yj , u, v)eι ω g(tj ,yj ,u,v) x(u, v), for each j = 1, ..., N 2 . Using this mode of subdivision, one has to be careful about selecting number of nodal points. 7
In this case, (xj , tj ) 6= (ζu , ζv ). Thus, one should select dense points around stationary phase point(s) to obtain good accuracy. Therefore, (11) can be written in form of (33) as: Ij = Ij1 + Ij2 + Ij3 + Ij4 , where Ij1 Ij2
Z
=
ζu
Z
=
Z
Z
b
=
Υj dvdu,
(35)
Υj dvdu,
(36)
Υj dvdu,
(37)
Υj dvdu.
(38)
ζv
c
Z
ζu
a
Z
d
ζv
ζu
Ij3 = Ij4
Z
b
(34)
d
ζv
Z
ζu
a
ζv
c
In the second case, one can divide the domain sub-interval in such a way that the stationary phase point(s) is enclosed in small sub-interval separating it from the other sub-intervals where the kernel function is free of stationary phase point(s). The u− and v− components of the stationary-point(s) (ζu , ζv ) satisfies ζu ∈ (ζu − ηu , ζu + ηu ) and ζv ∈ (ζv − ηv , ζv + ηv ) where ηu , ηv > 0 are assumed to be relatively small numbers. This technique of sub-division can be preferred to use, as it leads to good accuracy irrespective of uniform or non-uniform nodal distribution. Consequently, the integral part of the integral equation (2) can be written as: Z bZ a
d
Υj dvdu =
c
+ + =
Z
b
Z
d
Υj dvdu +
ζu +ηu ζv +ηv Z ζu +ηu Z d ζu −ηu Z ζu −ηu a Ij1
ζv +ηv Z d
Υj dvdu + Υj dvdu +
Z
Z
b
+ Ij2 + · · · + Ij9 ,
where Ij1
=
ζu −ηu Z ζu −ηu
Ij3 Ij4 Ij5
=
Z
Z
b
Z
b
ζu +ηu
Z
b
ζu +ηu
= =
ζv −ηv Z ζv +ηv
Z
Υj dvdu + Υj dvdu +
ζv −ηv
Z
b
ζu +ηu c Z ζu +ηu ζu −ηu Z ζu −ηu a
ζv −ηv
Z
Υj dvdu
ζv −ηv
Υj dvdu
c
Z
ζv −ηv
Υj dvdu,
c
(39)
ζu +ηu
Ij2 =
Υj dvdu +
ζu +ηu ζv −ηv Z ζu +ηu Z ζv +ηv a
ζv +ηv
ζv +ηv
Z
Z
ζu +ηu
ζu −ηu
Υj dvdu,
(40)
Υj dvdu,
(41)
Υj dvdu,
(42)
Υj dvdu,
(43)
Υj dvdu,
(44)
ζv +ηv
Z
ζv +ηv
ζv −ηv
Z
ζu +ηu
ζu −ηu
d
ζv −ηv
c
Z
d
ζv +ηv
Z
ζv +ηv
ζv −ηv
8
Ij6
=
Z
ζu +ηu
ζu −ηu
Ij7
=
Z
Z
Z
Ij9 =
Z
and
Υj dvdu,
(45)
Υj dvdu,
(46)
Υj dvdu,
(47)
Υj dvdu.
(48)
c
Z
ζu −ηu
a
Ij8 =
ζv −ηv
d
ζv +ηv
ζu −ηu
Z
ζv +ηv
ζv −ηv
a ζu −ηu
a
Z
ζv −ηv
c
The integrals Ij1 − Ij9 shown above are evaluated on the local nodes in their respective subdomains in the following sequel: Dom Ij1 = [ζu + ηu , b] × [ζv + ηv , b], Dom Ij2 = [ζu + ηu , b] × [ζv − ηv , ζv + ηv ], Dom Ij3 = [ζu + ηu , b] × [c, ζv − ηv ], Dom Ij4 = [ζu − ηu , ζu + ηu ] × [ζv + ηv , d], Dom Ij5 = [ζu − ηu , ζu + ηu ] × [ζv − ηv , ζv + ηv ], Dom Ij6 = [ζu − ηu , ζu + ηu ] × [c, ζv − ηv ], Dom Ij7 = [a, ζu −ηu ]×[ζv +ηv , d], DomIj8 = [a, ζu −ηu ]×[ζv −ηv , ζv +ηv ], and DomIj9 = [a, ζu −ηu ]×[c, ζv −ηv ]. The nodes of these sub-domains are generally different from the global nodes occurring in the main interval Dom Ij = [a, b] × [c, d]. We can evaluate all these integrals using the procedure of delaminating quadrature mentioned in subsections 3.2 and 3.3, which transforms a two- dimensional problem into one-dimensional problem. Consequently, one can easily handle the case of stationary phase point(s) in a two-dimensional integral equation by using solution approach of a one-dimensional integral equation. To do this, as a priori requirement, we establish connectivity between the local and global nodes through barycentric interpolation [30, 31]. Let x(vj ) be the given function values for vj , j = 1, 2, ..., N 2 then the function value at inner point can be interpolated as [31]: x(v) = where the weights are defined as: $j = Q
PN 2
$j j=1 v−vj x(vj ) , PN 2 $j j=1 v−vj
1 , j, k = 1, 2, ..., N 2 . (v − v ) j k k6=j
For the uniform nodes, the weights are simplified as: $j = (−1)j
N 2! , − j)!j!
(N 2
whereas for Chebyshev-type nodes, the simplified weights are: 1/2, j = 1 or j = N 2 ; j $j = (−1) δj , δj = 1, otherwise. (49) can be expressed in matrix form as: T x(ζ) = l1 (ζ), l2 (ζ), · · · , lN 2 (ζ) x(v1 ), x(v2 ), · · · , x(vN 2 ) , 9
(49)
where coefficients are lj =
$j v−vj PN 2 $ j=1 v−vj
. r2 local nodes ξk , k = 1, 2, · · · , r2 can be expressed as: l2 (ξ1 ) · · · lN 2 (ξ1 ) x(v1 ) l2 (ξ2 ) · · · lN 2 (ξ2 ) x(v2 ) = M l X. .. .. .. .. . . . .
Consequently, the function values on l1 (ξ1 ) x(ξ1 ) x(ξ2 ) l1 (ξ2 ) = .. .. . . x(ξr2 ) l1 (ξr2 )
l2 (ξr2 ) · · · lN 2 (ξr2 )
x(vN 2 )
For each i = 1, 2, · · · , 9, and ti = 1, 2, · · · , Ni2 ; (uiti , vtii ) be the new nodes in the sub-intervals Dom Iji . Connectivity between the global nodes vk , k = 1, 2, · · · , N 2 and local ones are given in the v− direction as: P11 l1 (v11 ) l2 (v11 ) · · · lN 2 (v11 ) P1 P12 l1 (v21 ) l2 (v21 ) · · · lN 2 (v21 ) P2 P1N 2 ×N 2 = .. = .. .. .. .. = Ol1 PN 2 ×N 2 , .. . . . . . . 1 1 ) l (v 1 ) · · · l l1 (vN 2 N2 2 N 2 (vN 2 )
P1N 2 1
P2N 2 ×N 2 =
P21 P22 .. . P2N 2 2
and so on
P9N 2 ×N 2
=
1
P91 P92 .. . P9N 2 9
l1 (v12 ) l1 (v22 ) .. .
PN 2
1
1
l2 (v12 ) l2 (v22 ) .. .
··· ··· .. .
lN 2 (v12 ) lN 2 (v22 ) .. .
P1 P2 .. .
= 2 2 2 l1 (vN 2 ) l2 (vN 2 ) · · · lN 2 (vN 2 ) PN 2 2
2
2
lN 2 (v19 ) P1 9 l 2 (v2 ) P2 N .. = .. . . 9 9 2 9 l1 (vN 2 ) l2 (vN 2 ) · · · lN (vN 2 ) PN 2
l1 (v19 ) l1 (v29 ) .. . 9
l2 (v19 ) l2 (v29 ) .. .
··· ··· .. .
9
9
= Ol2 PN 2 ×N 2 .
= Ol9 PN 2 ×N 2 ,
where the superscript 1 − 9 indicates the solution on the respective sub-domains; DomIji . The connectivity between global nodes uk , k = 1, 2, · · · , N 2 and local nodes in the u− direction can be expressed as: l1 (u11 ) l2 (u11 ) · · · lN 2 (u11 ) x(u11 ) x(u1 ) x(u12 ) l1 (u12 ) l2 (u12 ) · · · lN 2 (u12 ) x(u2 ) Ψ1 = .. .. .. = Ml1 xr , = .. .. ... . . . . . l1 (u1N 2 ) l2 (u1N 2 ) · · · lN 2 (u1N 2 ) x(u1N 2 ) x(uN 2 ) 1
Ψ2 =
and so on
x(u21 ) x(u22 ) .. .
1
l1 (u21 ) l1 (u22 ) .. .
1
1
l2 (u21 ) l2 (u22 ) .. .
··· ··· .. .
lN 2 (u21 ) lN 2 (u22 ) .. .
x(u1 ) x(u2 ) .. .
= 2 2 2 2 x(uN 2 ) l1 (uN 2 ) l2 (uN 2 ) · · · lN 2 (uN 2 ) x(uN 2 ) 2
2
2
2
10
= Ml2 xr .
Ψ9 =
lN 2 (u91 ) x(u1 ) lN 2 (u92 ) x(u2 ) .. = = Ml9 xr . .. . . 2 (u9 ) l1 (u9N 2 ) l2 (u9N 2 ) · · · lN x(u ) x(u9N 2 ) 2 2 N N x(u91 ) x(u92 ) .. .
l1 (u91 ) l1 (u92 ) .. . 9
l2 (u91 ) l2 (u92 ) .. . 9
··· ··· .. .
9
After establishing connectivity between global and local nodes, the numerical results of the integrals Ij1 − Ij9 can be calculated in the similar way as shown in subsections 3.2 and 3.3 of the same section. The integrals Ij1 , j = 1, 2, ..., N 2 can be written as: Du + ιωΛ1u Ij1 = Q1u j j
where
−1
−1 Ol1 P1 − Q1l Du + ιωΛ1l Ol1 PN 2 , j j
−1 1u −1 l 1l 1l l = Q1u D + ιωΛ O Γ − Q D + ιωΛ O Γ Ml1 xr 2 u 1 u j j 1 j j 1 N
(50)
(51)
ι ω g(t ,u ,b,d) T j j Q1u , 0, · · · , 0, −eι ω g(tj ,uj ,ζu +ηu ,d) , j = e
ι ω g(t ,u ,b,ζ +η ) T v v , j j Q1l 0, · · · , 0, −eι ω g(tj ,uj ,ζu +ηu ,ζv +ηv ) , j = e h i ∂g(tj ,uj ,uN 2 ,d) ∂g(tj ,uj ,u1 ,d) ∂g(tj ,uj ,u2 ,d) Λ1u = diag , and , , · · · , j ∂u ∂u ∂u h i ∂g(tj ,uj ,uN 2 ,ζv +ηv ) ∂g(tj ,uj ,u1 ,ζv +ηv ) ∂g(tj ,uj ,u2 ,ζv +ηv ) Λ1l = diag . , , · · · , j ∂u ∂u ∂u
Similar numerical procedure can be adopted for other integrals; Ij2 − Ij9 to obtain the numerical value of Ij given in (39), over the entire domain as: Ij = + + + + + + + +
−1 Ol1 Γ1 − Q1l Du + ιωΛ1l Ol1 ΓN 2 Ml1 xr j j −1 −1 l Du + ιωΛ2u O2 Γ1 − Q2l Du + ιωΛ2l Ol2 ΓN 2 Ml2 xr j j j −1 −1 l l l 3l 3l Du + ιωΛ3u O Γ − Q D + ιωΛ O Γ M 2 u 1 j 3 N 3 xr j 3 j −1 −1 l 4l 4l l l Du + ιωΛ4u O Γ − Q D + ιωΛ O Γ M 2 u j 4 1 j j 4 N 4 xr −1 −1 l Du + ιωΛ5u O5 Γ1 − Q5l Du + ιωΛ5l Ol5 ΓN 2 Ml5 xr j j j −1 −1 l Du + ιωΛ6u O6 Γ1 − Q6l Du + ιωΛ6l Ol6 ΓN 2 Ml6 xr j j j −1 −1 l Du + ιωΛ7u O7 Γ1 − Q7l Du + ιωΛ7l Ol7 ΓN 2 Ml7 xr j j j −1 −1 l 8l 8l l l Du + ιωΛ7u O Γ − Q D + ιωΛ O Γ M 2 1 u j 8 j j 8 N 8 xr −1 −1 l 9l 9l l l Du + ιωΛ9u O Γ − Q D + ιωΛ O Γ M 2 u j 9 1 j j 9 N 9 xr .
Q1u Du + ιωΛ1u j j
Q2u j Q3u j Q4u j Q5u j Q6u j Q7u j Q8u j Q9u j
−1
11
(52)
Letting Uj =
−1 Ol1 Γ1 − Q1l Du + ιωΛ1l Ol1 ΓN 2 Ml1 j j −1 −1 l 2l 2l l Du + ιωΛ2u O Γ − Q D + ιωΛ O Γ Ml2 2 u j 2 1 j j 2 N −1 −1 l l 3l 3l Ml3 O Γ Du + ιωΛ3u O Γ − Q D + ιωΛ u 3 N2 j 3 1 j j −1 −1 l Du + ιωΛ4u O4 Γ1 − Q4l Du + ιωΛ4l Ol4 ΓN 2 Ml4 j j j −1 −1 l Du + ιωΛ5u O5 Γ1 − Q5l Du + ιωΛ5l Ol5 ΓN 2 Ml5 j j j −1 −1 l Du + ιωΛ6u O6 Γ1 − Q6l Du + ιωΛ6l Ol6 ΓN 2 Ml6 j j j −1 −1 l 7l 7l l Du + ιωΛ7u O Γ − Q D + ιωΛ O Γ Ml7 2 u j 7 1 j j 7 N −1 −1 l l 8l 8l Ml8 O Γ Du + ιωΛ7u O Γ − Q D + ιωΛ u 8 N2 j 8 1 j j −1 −1 l Du + ιωΛ9u O9 Γ1 − Q9l Du + ιωΛ9l Ol9 ΓN 2 Ml9 . j j j
Q1u Du + ιωΛ1u j j
+ Q2u j + Q3u j + Q4u j + Q5u j + Q6u j + Q7u j + Q8u j + Q9u j
−1
(53)
The interpolation matrices Mli , i = 1, . . . , 9 depend on the nature of stationary phase point(s). These matrices Ml1 − Ml9 remain unchanged for fixed stationary phase point(s). By substituting the value of Uj in (52) we have; Ij = Uj xr ,
(54)
where Uj are complex vectors each of order 1 × N 2 and xr is a vector of order N 2 × 1. Putting the value of Ij from (54) into (10) and following the steps given in subsection 3.3 the required solution of the integral equation having stationary-point is obtained.
4
Convergence Analysis
Lemma 4.1. (V ander corput, Stein [32], p.332). Suppose g(y) is real valued and smooth in (a, b) and that g k (y) ≥ 1 ∀ x ∈ (a, b) for fixed value of k then Z b −1 ι ω g(y) e dy ≤ c(k)ω k
(55)
a
holds when 1. k ≥ 2, or
2. k = 1 and g 0 (y) is monotonic. The bound c(k) is independent of g and ω, c(k) = 5.2k−1 − 2. Lemma 4.2. (Stein [32], p.334). Under the assumption on g(y) in Lemma 4.1, it has been concluded that Z b Z b 0 ι ω g(y) ≤ c(k)ω −1 ϕ (b) dy k e ϕ(y)dy |ϕ(b)| + (56) a
a
12
In the procedure of a delaminating quadrature a 2D integral is evaluated by transforming it into a 1D inner and outer integrals. Let a general 1D oscillatory integral mentioned in sections 3.2 and 3.3 is of the form Z b
Ix =
K(y)eι ω g(y) dy,
(57)
a
where K(y) = f (t, y)x(t, y), g(y) = g(t, y). The approximate value of (57) using the Levin’s quadrature theory and the procedures of GDMP and/or LDMP at descretized nodes tk , k = 1, 2, ..., N 2 is given by Z b K(y)eι ω g(y) dy, I n xn = a Z b ∂ = p(y)eι ω g(y) dy, a ∂y h ib = p(y)eι ω g(y) , y=a
= p(b)e
ι ω g(b)
− p(a)eι ω g(a) ,
(58)
where K(y) = K(tk , y) = f (tk , y)x(tk , y), g(y) = g(tk , y) for each fixed k. The error bound of the discretized form xn − I n xn = hn are discussed in the following theorems. Theorem 4.3. Suppose that f (y) and g(y) are suitably smooth functions and g 0 (y) 6= 0∀y ∈ [a, b]. Then (57) and (58) at arbitrary interpolation nodes y0 , y1 , · · · , yM satisfy
3(M + 2) F M +1 (y) ∞ (b − a)M +1 E = |Ix − I n xn | ≤ . (59) ω(M + 1)!
where F (y) = K(y) − L(y), L(y) ≡ p0 (y) + ιωg 0 (y)p(y). More specifically, when the shape parameter tends to zero, then (59) reduces to
3(M + 2)(b − a)M +1 K M +1 (y) ∞ E = |Ix − I n xn | ≤ . ω(M + 1)!
(60)
Proof. Let p(y) =
M X
λj Φ (r).
(61)
p 1 + 2 r2 ,
(62)
j=0
where Φ (r) is MQ-RBF defined as: Φ (r) =
and r2 = (y − yjc )2 , j = 0, 1, ..., M. According to Levin’s quadrature theory we have
L(y) ≡ p0 (y) + ιωg 0 (y)p(y), then L(y) ≡
M X j=0
λj Φ0 (r)
0
+ ιωg (y)
M X j=0
13
λj Φ (r), j = 0, 1, · · · , M,
(63)
(64)
where the interpolation nodes y0 , y1 , · · · , yM ∈ [a, b]. We have Z b ιωg(y) dy , E = |Ix − I n xn | = (K(y) − L(y))e a Z b ιωg(y) = F (y)e dy ,
(65) (66)
a
where F (y) = K(y) − L(y) satisfies the interpolation condition F (yj ) = 0, j = 0, 1, · · · , M . One can find dj ∈ (yj , yj+1 ) 3 F 0 (dj ) = 0, j = 0, 1, · · · , M − 1, then F (y) and F 0 (y) can be written in terms of Lagrange interpolation form as: M F M +1 (ξ1 ) Y F (y) = (y − yjc ), (M + 1)! j=0
and F 0 (y) =
M −1 F M +1 (ξ2 ) Y (y − yjc ), (M )! j=0
for some ξ1 , ξ2 ∈ [a, b]. By Lemma 4.1 and 4.2, for k = 1, c(k) = 3, we have
≤
Z
F 0 (y) dy, a
M +1 M +1
F 3(M + 2)(b − a) (y) ∞ . ω(M + 1)!
E = |Ix − I n xn | ≤ c(k)ω
−1/k
(|F (b)| +
b
(67)
When the shape parameter → 0, the interpolant L(y) converge to Lagrange interpolating polynomial [33, 34] i.e
LM +1 (y) ≡ 0 ⇒ F M +1 (y) ∞ = K M +1 (y) ∞ and above (67), then becomes
3(M + 2)(b − a)M +1 K M +1 (y) ∞ E = |Ix − I n xn | ≤ . ω(M + 1)!
(68)
Definition 1. (Kress [35],p.288). A linear operator I : X → Y from a normed space X into a normed space Y is called compact if for each bounded sequence (xn ) in X the sequence (Ixn ) contains a convergent subsequence in Y , i.e., if each sequence from the set {Ix : x ∈ X, kxk ≤ 1} contains a convergent subsequence. Definition 2. (Kress [36],p.5). A subset A of a normed space X is called complete if every Cauchy sequence of elements in A converges to an element in A. A normed space X is called a Banach space if it is complete.
14
Theorem 4.4. (Kress [35],p.292). Suppose that I : X → X be a compact linear operator on a Banach space X such that I − I is injective. Assume that the sequence I n : X → X of bounded
linear operators is norm convergent i.e I n − I → 0, n → ∞. Then for sufficiently large n the inverse operators (I − I n )−1 : X → X exists and are uniformly bounded. For solution of the equations x − Ix = f and xn − I n xn = fn we have an error estimate kxn − xk ≤ C{kI n x − Ixk + kfn − f k},
(69)
for some constant C. Theorem 4.5. Let I : [a, b] → [a, b] be a compact linear operator such that it satisfies the statement of the Theorem 4.4. The numerical approximation of the operator equation x−Ix = f is given by xn − I n xn = fn solved through meshless methods, either by the GDMP or the LDMP. Then the error bound at M + 1 collocation points satisfies (b − a)M +1 kxn − xk = O . (70) ω Proof. By Theorem 4.4, we have an error estimate kxn − xk ≤ C1 {kI n x − Ixk + kfn − f k},
(71)
for some constant C1 . Since f is the known function therefore both f and fn agree at the given interpolation points. (71) then becomes kxn − xk ≤ C1 {kI n x − Ixk ,
≤ C1 kI n − Ik kxk .
(72) (73)
Using the definition (1), for kxk ≤ 1 we have kxn − xk ≤ C1 kI n − Ik .
(74)
Incorporating the error bound of the integral operator given in Theorem 4.3, (74) simplified to
3(M + 2)(b − a)M +1 K M +1 (y) ∞ kxn − xk ≤ C1 , ω(M + 1)!
(b − a)M +1
K M +1 (y) , = C ∞ ω
1 (M +2) where C = 3C(M +1)! . By the uniform boundedness of K M +1 [11, 37], we obtained (b − a)M +1 kxn − xk = O . ω
This completes the proof.
15
(75) (76)
(77)
5
Numerical Experiments
In this section we have provided examples of oscillatory integral equations having kernel function with and without stationary phase point(s). The benchmark problems having kernels free of stationary phase point(s) are taken from [1]. Performance of the proposed methods is checked on both uniform and non-uniform points distribution. Furthermore, we have also checked efficiency of the proposed methods on non-uniform nodes. The non-uniform nodes are condensed around stationary-point(s) [1] which are useful to alleviate the negative effects of stationary phase point(s) on the accuracy of the method. Other nodal points distribution with uniform boundary points and scattered interior points is shown in Fig. 4. Similarly, Chebeshev-type of nodal points distribution is shown in Fig. 2. MQ-RBF depends on the value of the shape parameter. In general, we have taken this value to be = 0.5 in this paper. The accuracy is measured in terms of maximum relative error norm which is defined as follows: Ljab = |x(tj , yj ) − x ˜(tj , yj )|, Lre =
max |{z}
j=1,...,N 2
j = 1, · · · , N 2 ,
Ljab , |x(tj , yj )|
where the exact solution is represented by x and the computed solution by x ˜. Computations are executed using Core I5. The technique of the LDMP depends upon stencil size, which has been taken as 3 throughout the computation. Figure 1 is for Test Problem 1. Figures 2 and 3 are of Test Problem 2. Figures 4 and 5 are related to Test Problem 3. Figures 6-9 are related to Test Problem 4. Test Problem 1. Consider an oscillatory integral equation [1]; x(t, y) = h(t, y) +
Z
1
−1
Z
1
3 ev
−1
√ eι ω (3 (u + ecos( u + v)
7y 5
) − 5 t v) x(u, v) dvdu,
where h(t, y) = e
√ cos( t + y)
21
+
2e5
ιωy−5ιωt
sin(ω)(4 cos2 (ω) − 1) (− e− 1 + 10 ι ω t + e) . ω (− 1 + 5 ι ω t) √
This example has the exact solution x(t, y) = ecos( t+y) . The integral equation has an oscillatory kernel function free of stationary phase point(s). Fig. 1 (left) shows convergence in terms of relative error Lre of the proposed approach against increasing frequencies. The Lre is decreasing exponentially in the range 10−6 − 10−16 for fixed mesh-points N 2 = 16 and frequency range 104 − 109 . Fig. 1 (right) shows behavior of the Lre where N 2 varies in the range 16 − 324 for fixed frequency ω = 106 . The error norm Lre is not much influenced by increase in the nodal points and is maintained around O(10−11 ). It is witnessed that the performance of proposed methods: the GDMP and the LDMP is very good on a small number of points. Test Problem 2. Consider the integral equation [1]; x(t, y) = h(t, y) +
Z
1
−1
Z
1
−1
et + u ι ω (t u + v) e x(u, v) dvdu, cos(u) 16
−6
−9
10
10 LDMP GDMP
−8
10
−10
Lr e
Lr e
10
−10
10
−12
10
−14
10
−11
10
−16
10 10000
100000
1000000
10000000
100000000
1000000000
16
36
64
100
144
196
256
324
N2
ω
N 2 = 16
ω = 106
Figure 1: Test Problem 1, accuracy in terms of Lre versus ω, (left) and N 2 (right).
where
ι et −ι ω − ι ω t e−1 − e1 + 2ι ω t h(t, y) = cos(t) − ω (1 + ι ω t)
−1 + e2 ι ω
.
Exact solution of the integral equation is x(t, y) = cos(t). This problem is solved using Chebeshev-type of nodal points distribution as shown in Fig. 2. In Fig. 3 (left), accuracy of proposed algorithms the GDMP and the LDMP is shown against the increasing frequencies for fixed N 2 = 100. One can observe from this figure that as long as the frequencies are increasing in the range of 104 − 109 , accuracy of the proposed method increases in the range 10−4 − 10−14 . It is also witnessed from Fig. 3 (right) that the accuracy is retained with increasing nodal points. We have also observed that the proposed algorithms perform well in the case of Chebeshev-type nodes as well. Test Problem 3. Let an integral equation be of the form: Z 1Z 1 5 7v cos v x(t, y) = h(t, y) + eι ω (5 u − 3 − 2 ) x(u, v) dvdu, −1 −1 v + cos u where h(t, y) = y+cos t−2
7e
−5ιω 3
ω − 7e
10ιω 3
ω − 7ω cos(1)e
−31ιω 6
+ 7ω cos(1)e 5ω(49ω 2 − 4)
−ιω 6
− 2ι sin(1)e
−31ιω 6
+ 2ι sin(1)e
The exact solution of the problem is x(t, y) = y + cos t. We have considered this example to justify the results over another type of scattered nodal points as shown in Fig. 4 for N 2 = 100. Accuracy of the proposed algorithms GDMP and LDMP in terms of Lre is in the range of 10−7 −10−16 against increasing frequencies as shown in Fig 5 (left). This figure also confirms that whatsoever the arrangements of nodal points might be, the new proposed algorithms converge and hence are accurate. In Fig. 5 (right), we observed some mild fluctuations in the accuracy against increasing nodes for fixed frequency ω = 106 . 17
−ιω 6
;
1 0.8 0.6 0.4
y−axis
0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1
−0.8
−0.6
−0.4
−0.2
0 x−axis
0.2
0.4
0.6
0.8
1
N 2 = 324 Figure 2: Test Problem 2, model domain of a Chebeshev-type nodal points.
−4
−10
10
10 LDMP GDMP
−6
10
−8
Lr e
Lr e
10
−11
10
−10
10
−12
10
−12
10
−14
10 10000
100000
1000000
10000000
100000000
1000000000
16
36
64
100
144
196
256
324
N2
ω
N 2 = 100
ω = 106
Figure 3: Test Problem 2, accuracy in terms of Lre versus ω, (left) and N 2 (right) incorporating Chebyshev-type nodal points.
18
1 0.8 0.6 0.4
y−axis
0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1
−0.8
−0.6
−0.4
−0.2
0 x−axis
0.2
0.4
0.6
0.8
1
N 2 = 100 Figure 4: Test Problem 3, model domain of a non-uniform data.
−7
−10
10
10 LDMP GDMP
−8
10
−9
10
−10
10
−11
Lr e
Lr e
10
−11
10
−12
10
−13
10
−14
10
−15
10
−12
10
−16
10 10000
100000
1000000
10000000
100000000
16
36
64
100
144
196
256
324
N2
ω
N 2 = 100
ω = 106
Figure 5: Test Problem 3, accuracy in terms of Lre versus ω, (left) and N 2 (right) incorporating non-uniform data.
19
Test Problem 4. Let the Fredholm integral equation be Z 1Z 1 (sin t + e2v + y)eu+v ι ω (u2 + v2 ) x(t, y) = h(t, y) + e x(u, v) dvdu, cos u + e2u+v −1 −1 where ι √ √ √ ι ι ιω − 1 ιω + 1 −ιω π −2 sin t erf −ιω + e ω erf √ − 2y erf −ιω + e ω erf √ h(t, y) = x(t, y)+ erf 2ω −ιω −ιω 2t+y
The exact solution of the problem is x(t, y) = cos(t)+e . In this example the kernel function et+y contains the stationary phase point at 0. To handle such type of problem, we follow the procedure explained in section 3.4. According to the first strategy of domain sub-division, the oscillatory integral having kernel with stationary point 0 can be written as: Z 1Z 1 Z 1Z 1 Z 1Z 0 Z 0Z 1 Z 0Z 0 χ dvdu = χ dvdu + χ dvdu + χ dvdu + χ dvdu, −1
−1
0
0
(sin t+e2v +y)eu+v
0
2
−1
−1
−1
0
−1
2
where χ = cos u+e2u+v eι ω (u + v ) x(u, v). We have tested this mode of domain subdivision both on uniform and non-unform nodal points. Accuracy in terms of Lre increases with the increase of frequency on fixed non-uniform scattered nodal point(s) that are condensed around the stationary point 0 as shown in Fig. 6 (left). In Fig. 6 (right), we observed that for the same number of uniform nodal points, the Lre falls behind the scattered data points for increasing frequencies. In the second case, one can divide the domain interval in such a way that the stationary-point 0 ∈ (−η, η) (η > 0 is a small real number). Consequently, the integral in hand can be written as: Z 1Z 1 Z 1Z 1 Z 1Z η Z 1 Z −η Z ηZ 1 Z ηZ χ dvdu = χ dvdu + χ dvdu + χ dvdu + χ dvdu + −1
−1
+
η Z η
−η
η
Z
η
−η
−1
χ dvdu +
Z
−η −η Z 1
−1
η
χ dvdu +
Z
−1 −η
−1
η
Z
η
−η
−η
χ dvdu +
Z
η −η
−1
Z
−η
−η
χ dvdu.
−1
For this scheme of intervals sub-division, the proposed meshless procedure GDMP and LDMP produces accurate results as shown in Fig. 7. Accuracy in terms of Lre varies between 10−4 to 10−9 against increasing frequencies ω and fixed nodal points N 2 = 36, both for uniform and non-unform nodal points. Performance of the proposed meshless methods, GDMP and LDMP are equally good for uniform and non-uniform nodal points. Furthermore, we have observed from Fig. (9) that both the methods GDMP and LDMP are efficient on coarse grid as well.
6
Conclusion
The proposed solution algorithms are applied on different models of oscillatory integral equations. These models are implemented on uniform and non-uniform nodal points. The models given in test problems 1-3 are examples of Fredholm integral equations having oscillatory kernel function free of stationary phase point(s). The Fredholm integral equations model given in test problem 4 is an example where the oscillatory kernel function is having stationary phase point(s). The 20
η
−η
χ dvdu
−3
−2
10
10 LDMP GDMP
LDMP GDMP
−4
10
−3
10
−5
10
−4
Lr e
Lr e
10 −6
10
−5
10 −7
10
−6
10
−8
10
−9
10 10000
−7
100000
1000000
10000000
100000000
10 10000
1000000000
100000
1000000
10000000
ω
ω
N 2 = 16
N 2 = 16
100000000
1000000000
Figure 6: Test Problem 4, accuracy in terms of Lre versus ω, incorporating non-uniform data condensing around the stationary phase point (left) and uniform data (right) after subdivision of the interval at the position of stationary point.
−3
−3
10
10 LDMP GDMP
LDMP GDMP
−4
−4
10
10
−5
−5
10
Lr e
Lr e
10
−6
10
−7
−7
10
10
−8
−8
10
10
−9
10 10000
−6
10
−9
100000
1000000
10000000
100000000
1000000000
10 10000
100000
1000000
10000000
ω
ω
N 2 = 36
N 2 = 36
100000000
1000000000
Figure 7: Test Problem 4, accuracy in terms of Lre versus ω, incorporating non-uniform data condensing around the stationary point (left) and uniform data (right) after subdivision of the interval into sub-intervals inclosing stationary point.
21
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8 −1
−0.8
−0.6
−0.4
−0.2
N2
0
0.2
0.4
0.6
0.8
1
= 36
Figure 8: Test Problem 4, non-uniform data condensing around the stationary phase point model.
4
10
Uniform nodes−LDMP Scattered nodes−LDMP Uniform nodes−GDMP Scattered nodes−GDMP 3
time
10
2
10
1
10
36
81
144
N2
Figure 9: Test Problem 4, CPU time(s) versus N 2 .
22
proposed methods GDMP and LDMP give better accuracy for larger frequencies on a coarse grid for integral equation free of stationary phase point(s). In case when the kernel function incorporate stationary phase point(s), the domain interval is divided into sub-intervals according to the position of stationary phase point(s). The existence of stationary phase point(s) in the kernel of oscillatory integral equation offer difficulties in its solution. Levin’s type methods [19] fail in such scenario. To handle such type of issues in oscillatory integral equation models, we implemented the proposed methods on two types of domain sub-division. In the first type of domain sub-division, the proposed methods are converging fast when the nodes are non-uniform and condensing around the stationary phase point(s) but diverging after initial convergence when the nodes are uniform. In the second case of domain sub-division, the proposed methods work accurately both for uniform and non-uniform nodal points discretization. Thus, the proposed algorithms GDMP and LDMP are accurate in both the cases, with and without stationary phase point(s).
References [1] Siraj-ul-Islam, I. Aziz, Zaheer-ud-Din, Meshless methods for multivariate highly oscillatory Fredholm integral equations, Eng. Anal. Bound. Elem. 53 (2015) 100–112. [2] G. P. Galdi, K. Pileckas, A. L. Silvestre, On the unsteady poiseuille flow in a pipe., Z. Angew. Math. Phys. 58 (2007) 994–1007. [3] M. A. Bartoshevich, A heat-conduction problem., J. Eng. Phy. Thermophy. 28 (1975) 240– 244. [4] P. Baratella, A Nystrom interpolant for some weakly singular linear Volterra integral equations, Comput. Appl. Math. 231 (2009) 725–734. [5] A. M. Wazwaz, Linear and Non Linear Integral Equations: Methods and Appications, Higher Education Press, Beijing and Springer-Verlog Berlin Heidekberg, Germany, 2011. [6] K. Atkinson, The numerical solution of integral equations of the second kind, Cambridge University Press, Cambridge, 1997. [7] A. J. Jerri, Introduction to Integral Equations with Applications, John Wiley and Sons, INC, 1999. [8] H. Brunner, Y. Ma, Y. Xu, The oscillation of solutions of volterra integral and integrodifferential equations with highly oscillatory kernels, J. Integ. Eq. Appl. 27 (4) (2015) 455– 487. [9] C. Geuzaine, O. Bruno, F. Reitich, On the O(1) solution of multiple-scattering problems, IEEE Trans. Magn. 41 (5) (2005) 1488–1491. [10] O. Bruno, C. Geuzaine, J. Monro, F. Reitich, Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case, R. Soc. Lond. Trans. Ser. 362 (1816) (2004) 629–645. [11] Siraj-ul-Islam, S. Zaman, New quadrature rules for highly oscillatory integrals with stationary points, J. Comput. Appl. Math. 278 (2015) 75–89.
23
[12] H. Guoqiang, W. Jiong, Extrapolation of Nystrom solution for two-dimensional nonlinear Fredholm integral equations, J. Comput. Appl. Math. 134 (2001) 259–268. [13] W. J. Xie, F. R. Lin, A fast numerical solution of two-dimensional Fredholm integral equations of the second kind, Appl. Math. Comput. 59 (2009) 1709–1719. [14] A. Iserles, S. P. Norsett, Efficient quadrature of highly oscillatory integrals using derivatives, Proc. R. Soc. 461 (2005) 1388–1399. [15] D.Huybrechs, S.Vandewalle, A sparse discretisation for integral equation formulations of high frequency scattering problems, Tech. Rep., M.U.Leuven, 2006. [16] D. Huybrechs, S. Vandewalle, A sparse discretization for integral equation formulations of high frequency scattering problems, SIAM J. Sci. Comput. 29 (6) (2007) 2305–2328. [17] J. Li, X. Wang, T. Wang, A universal solution to one-dimensional highly oscillatory integrals, Sci. China. Ser. 51 (10) (2008) 1614–1622. [18] J. Li, X. Wang, T. Wang, C. Shen, Delaminating quadrature method for multi-dimensional highly oscillatory integrals, Appl. Math. Comput. 209 (2) (2009) 327–338. [19] D. Levin, Procedures for computing one and two-dimensional integrals of functions with rapid irregular oscillations, Math. Comput. 38 (1982) 531–538. [20] Siraj-ul-Islam, I. Aziz, W. Khan, Numerical integration of multi-dimensional highly oscillatory, gentle oscillatory and non-oscillatory integrands based on wavelets and radial basis functions., Eng. Anal. Bound. Elem. 36 (2012) 1284–1295. [21] J. Li, X. Wang, S. Xiao, T. Wang, A rapid solution of a kind of 1D Fredholm oscillatory integral equation, J. Comput. Appl. Math. 236 (2012) 2696–2705. [22] F. Ursell, Integral equations with a rapidly oscillating kernel, J. Lond. Math. Soc. 44 (1969) 449–459. [23] Q. Shen, Local rbf-based differential quadrature collocation method for the boundary layer problem, Eng. Anal. Bound. Elem. 34 (2010) 213–228. [24] S. Sarra, A local radial basis function method for advection-diffusion-reaction equations on complexly shaped domains, Appl. Math. Comput. 212 (2012) 9853–9865. [25] G. Liu, Meshfree Methods: Moving Beyond the Finite Element Method, CRC Press, 2003. [26] C. Shu, H. Ding, K. Yeo, Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations, Comp. Meth. Appl. Mech. Eng. 192 (2003) 941–954. ˇ [27] B. Sarler, R. Vertnik, Meshfree explicit local radial basis function collocation method for diffusion problems, Comput. Math. Appl. 51 (2006) 1269–1282. ˇ [28] G. Kosec, B. Sarler, Meshless approach to solving freezing with natural convection, Mater. Sci. Forum. 22 (2010) 205–210. ˇ [29] Siraj-ul-Islam, R. Vertnik, B. Sarler, Radial basis function collocation method for the numerical solution of the two-dimensional transient nonlinear coupled Burger’s equations, App. Math. Model. 36 (3) (2012) 1148–1160. 24
[30] D. Kincaid, W. Cheney, Numerical Analysis, Vol. 2, Math sci comp, Brooks/Cole, Pacific Grove, 2002. [31] J. Berrut, L. Trefethen, Barycentric lagrange interpolation, SIAM Review 46 (3) (2004) 501–517. [32] E. Stein, Harmonic Analysis: Real-rariable Methods, Orthogonality, and Oscillatory Integrals, Princeton University Press, Princeton, NJ. [33] T. Driscoll, B. Fornberg, Interpolation in the limit of increasingly flat radial basis function, Comput. Math. Appl. 43 (2000) 3–5. [34] G. Song, J. Riddle, G. E. Fasshauer, F. J. Hickernell, Multivariate interpolation with increasing flat radial basis functions of finite smoothness, Ad. Comput. Math. 36 (2012) 485–501. [35] R. Kress, Numerical Analysis, Springer-Verlag, New York, Inc., 1998. [36] R. Kress, Linear Integral Equations, Springer-Verlag, Berlin, 1989. [37] S. Xiang, Efficient quadrature for highly oscillatory integrals involving critical points, J. Comput. Appl. Math. 206 (2007) 688–698.
25