- Email: [email protected]
Fast method of approximate particular solutions using Chebyshev interpolation
Engineering Analysis with Boundary Elements 64 (2016) 290–294
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Fast method of approximate particular solutions using Chebyshev interpolation A.R. Lamichhane a, D.L. Young b, C.S. Chen a,c,n a
Department of Mathematics, University of Southern Mississippi, USA Department of Civil Engineering, National Taiwan University, Taiwan c College of Mathematics, Taiyuan University of Technology, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 30 October 2015 Received in revised form 20 December 2015 Accepted 21 December 2015 Available online 14 January 2016
The fast method of approximate particular solutions (FMAPS) is based on the global version of the method of approximate particular solutions (MAPS). In this method, given partial differential equations are discretized by the usual MAPS and the determination of the unknown coefficients is accelerated using a fast technique. Numerical results confirm the efficiency of the proposed technique for the PDEs with a large number of computational points in both two and three dimensions. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Radial basis functions Meshless method Method of approximate particular solutions Chebyshev interpolation Fast summation method Fast method of approximate particular solutions
1. Introduction During the past two decades, the radial basis functions (RBFs) have been widely applied for solving partial differential equations (PDEs). Kansa [13] first proposed the so-called RBF collocation method in 1990 for computationally solving fluid dynamic problems. One of the attractions of the Kansa method is its simplicity for solving problems in high dimensional and irregular domains. Due to the popularity of the Kansa method, several other RBF collocation methods have been proposed. Among them, the method of approximate particular solutions (MAPS) [5,6] is another effective method using the particular solution of the given RBFs as the basis function. The formulation of the MAPS produces a full and dense matrix system which is often very ill-conditioned. Traditionally, this matrix system is solved by using direct or iterative methods. Direct methods such as Gaussian elimination require OðN3 Þ operations for an N N system of equations. For iterative methods, we may obtain the approximate solution in k steps with each step needing a matrix vector multiplication OðN 2 Þ. In recent years, the MAPS has been widely applied to solve physical and engineering problems like Navier–Stokes equations [3], Stokes flow problems [4], Linear elasticity equations [2], Timen
Corresponding author. E-mail address: [email protected] (C.S. Chen).
http://dx.doi.org/10.1016/j.enganabound.2015.12.015 0955-7997/& 2016 Elsevier Ltd. All rights reserved.
fractional diffusion equations [10], Inverse problem of nonhomogeneous convection-diffusion equation [12], Diffusion equation with non-classical boundary [1] and Nonhomogeneous cauchy problems of elliptic PDEs [15]. Simulations of these kinds of problems involve a large number of interpolation points. The high computational cost using traditional solvers has become an issue. In this paper we pay special attention on how to develop a fast algorithm to alleviate the issue of high cost for solving large-scale problems using the MAPS. Consequently, we propose to couple the MAPS with a fast summation method to reduce the computational time by multiplying a matrix and a vector in each step inside the iterative method. This fast summation method is based on the Chebyshev interpolation [9]. As we shall see in the section of numerical results, we have successfully solved a 2D problem using 694,541 collocation points with only 77.15 s of computer running time and 343,000 collocation points with 105.15 s for a 3D problem. Moreover, we do not compromise the accuracy for our proposed fast computation. The structure of the paper is as follows. In Section 2, we give a brief review of the MAPS and the closed-form particular solution of the Gaussian for Laplacian in 2D and 3D. In Section 3, we review the algorithm of fast summation method (FSM). In Section 4, we propose the fast method of approximate particular solutions (FMAPS) by coupling the FSM and the MAPS as a fast algorithm for solving PDEs which require a large number of collocation points. A
A.R. Lamichhane et al. / Engineering Analysis with Boundary Elements 64 (2016) 290–294
specific algorithm has been given. In Section 4.1, to demonstrate the efficiency of the proposed method, two numerical examples in 2D and 3D are given. Finally, a short conclusion is drawn in Section 4.1.1.
2. Method of approximate particular solutions (MAPS)
ΔuðxÞ ¼ f ðxÞ; x A Ω;
ð1Þ
uðxÞ ¼ gðxÞ; x A ∂Ω;
ð2Þ
where Δ is the Laplacian operator, Ω and ∂Ω are the interior and boundary of the computational domain, respectively. Suppose f xi gni¼ 1 are the interpolation points containing ni interior points in Ω and nb boundary points on ∂Ω; i.e., n ¼ ni þnb . Let ϕ be a given radial basis function. By the MAPS [6], we assume the solution to (1) and (2) can be approximated by ^ uðxÞ uðxÞ ¼
αj Φð J x xj J Þ;
ð3Þ
where J J is the Euclidean norm, fαj g are the undetermined coefficients, and
ΔΦ ¼ ϕ:
ð4Þ
By the collocation method, from (1) and (2), we have
αj ϕ J xi xj J ¼ f ðxi Þ;
1 r i rni ;
ð5Þ
ni þ 1 r i r n:
ð6Þ
j¼1 n X
αj Φ J xi xj J ¼ gðxi Þ;
Z
1
EiðxÞ ¼ x
eu du; u
ð11Þ
and γ C 0:5772156649015328 is the Euler–Mascheroni constant [11]. Note that EiðxÞ is the special function known as the exponential integral function [11]. Theorem 2. Let ϕðrÞ ¼ expð cr 2 Þ, and ΔΦðrÞ ¼ ϕðrÞ in 3D. Then, 8 pffiffiffiffi pffiffiffi π > > < 3=2 erf ð crÞ; r a 0; 4c r ð12Þ ΦðrÞ ¼ 1 > > : ; r ¼ 0; 2c where erfðxÞ is the special function defined as follows [11]: Z x 2 2 erfðxÞ ¼ pffiffiffiffi e u du:
π
j¼1
n X
Theorem 1. Let ϕðrÞ ¼ expð cr 2 Þ, c 4 0, and ΔΦðrÞ ¼ ϕðrÞ in 2D. Then, 8 1 1 > 2 > < Eiðcr Þ þ log ðrÞ; r a0; 4c 2c ð10Þ ΦðrÞ ¼ 1 > > ðγ þ log ðcÞÞ; r ¼ 0; : 4c where
For simplicity, let us consider the following Poisson equation with Dirichlet boundary condition
n X
291
ð13Þ
0
In this paper, we adopt the Gaussian, ϕ, as the RBF basis function and the corresponding particular solutions, Φ, as the basis functions for the approximation of the partial differential equation. The particular solutions of the Gaussian in (10) and (12) contain the special functions, EiðxÞ and erf ðxÞ, which is costly in terms of numerical evaluation. The efficiency can be significantly improved using compiled MATLAB MEX functions. In the next section, we briefly introduce a fast summation method for the matrix vector multiplication.
j¼1
From (5) and (6), we can formulate a linear system of equations Aα ¼ F
Consider the evaluation of the sum of the form
where " #
sðxÞ ¼
ϕ A¼ Φ
ni þ1 r k r n;
; ij
Φ ¼ Φ J xk xj J
; kj
1 ri r ni ; 1 r j rn;
α ¼ ½α1 α2 ⋯ αn T ;
F ¼ ½f ðx1 Þ ⋯ f ðxni Þ gðxni þ 1 Þ ⋯ gðxn ÞT . For a more general form of the PDEs such as the following convection-diffusion-reaction equation in 2D, ∂u ∂u þ gðxÞ þ hðxÞu ¼ lðxÞ; x ¼ ðx; yÞ A Ω; ∂x ∂y
ð8Þ
bj κ j j x xj j j ;
ð14Þ
αj ðkϕð J x x J Þ þ f ðxÞΦx J x xj J þ gðxÞΦy J x xj J
P M ðξ; ηÞ ¼
1 1 2 MX þ T ðξÞT i ðηÞ; M M i¼1 i
ð15Þ
MD
j¼1
þ hðxÞΦ J x xj J Þ ¼ lðxÞ:
where κ is either the RBFs or the particular solution of the corresponding RBFs. We can evaluate the sum (14) in an efficient way by using the Chebyshev interpolation technique as described in [9]. From [9], let
where ξ; η A ½ 1; 1; T i is the first kind Chebyshev polynomial of order i. Let H be a hypercube in D dimension which contains all the
we have [5] N X
n X j¼1
ϕ ¼ ϕ J xi xj J
kΔu þ f ðxÞ
3. Fast summation method (FSM)
ð7Þ
ð9Þ
One of the key procedures in the implementation of the MAPS is to obtain the closed-form expression for the particular solutions of the corresponding RBFs. The derivation of the particular solutions for the well-known RBFs has already been known [7,16,18– 20]. Recently, the closed-form particular solutions of the Gaussian have been derived in [14], which are stated as follows:
given collocation points fxi gni¼ 1 . Consider fξl gl ¼ 1 , fηl gM l ¼ 1 be two identical sets of Chebyshev nodes in ½ 1; 1D . Then by using linear transformations, we can map xi ; xj into ξi ; ηj , and ξl ; ηm in ½ 1; 1D into xl ; xm in H, respectively. Instead of directly evaluating (14), we approximate it by using the following functional approximation XX κ ð j j x yj j Þ ¼ k j j xl ym j j Q M ðξ l ; ξÞQ M ðηm ; ηÞ; ð16Þ l
D
m
where ξ; η are the linear transformations of x; y, respectively, and
292
A.R. Lamichhane et al. / Engineering Analysis with Boundary Elements 64 (2016) 290–294
in 2D, Q M ðξ; ηÞ ¼ P M ðξ1 ; η1 ÞP M ðξ2 ; η2 Þ;
ξ ¼ ðξ1 ; ξ2 Þ; η ¼ ðη1 ; η2 Þ;
and in 3D, Q M ðξ; ηÞ ¼ P M ðξ1 ; η1 ÞP M ðξ2 ; η2 ÞP M ðξ3 ; η3 Þ; η ¼ ðη1 ; η2 ; η3 Þ:
ξ ¼ ðξ1 ; ξ2 ; ξ3 Þ;
Using (16), (14) can be separated into the product of the different sums as follows: n X
sðxi Þ ¼
n X bj κ j j xi xj j j bj
j¼1
j¼1
M M X X D
4.1. Algorithm
!
D
κ ð j j xl xm j j ÞQ M ðξl ; ξi ÞQ M ðηm ; ηj Þ
l¼1m¼1
¼
MD X l¼1
Q M ðξ l ; ξ i Þ
MD X
κ ð j j xl xm j j Þ
m¼1
n X
bj Q M ðηm ; ηj Þ:
ð17Þ
Wm ¼
bj Q M ðηm ; ηj Þ; m ¼ 1; …; M ; D
and the computational cost for these Chebyshev weights is OðM D nÞ. Then, we compute the sðxÞ at the Chebyshev nodes xl as M X
Find a hypercube H in D dimension which contains fxi gni¼ 1 by using linear transformations, find ξi ; ηj in ½ 1; 1D corre-
sponding to xi ; xj , respectively, find xl ; xm in H corresponding to ξl ; ηm , respectively, compute the following matrices: R ¼ ðRÞm;j ¼ Q M ðηm ; ηj Þ;
IR ¼ ðI R Þl;i ¼ Q M ðξl ; ξi Þ;
m ¼ 1; 2; …; M D ; j ¼ 1; 2; …; n; l ¼ 1; 2; …; M D ; i ¼ 1; 2; …; ni ;
BR ¼ ðBR Þl;i ¼ Q M ðξl ; ξi Þ;
D
W m κ ð j j xl xm j j Þ; l ¼ 1; …; M D ;
l ¼ 1; 2; …; M D ; i ¼ ni þ1; ni þ 2; …; n;
IK ¼ ðI K Þl;m ¼ ϕðj j xl xm j j Þ;
m¼1
with computational cost OðM 2D Þ. Eventually, we compute the sðxÞ at the collocation points xi as M X
4.1.2. Step 1: pre-computational step
j¼1
sðxl Þ ¼
4.1.1. Input Collocation points fxi gni¼ 1 , Chebyshev nodes ξl ; ηm , dimension of the domain D.
j¼1
If we choose those M D Chebyshev nodes such that M D o o n, then the summation can be computed efficiently. To evaluate (14) approximately, let us first evaluate the last term of (17) which gives us the weights at the Chebyshev nodes xm n X
this, we can use any iterative method which includes the vector multiplication with matrix A. In this paper, we have used GMRES iterative method to find the unknown coefficients. As we know that, in the GMRES, at each iteration, we have to multiply a matrix and an updated vector. This matrix and vector multiplication is obtained faster by the technique that we have described earlier. This kind of approach which is called FMAPS, actually solves the PDEs by using MAPS without explicitly computing matrix A which requires a minimum storage. The computational procedure for the FMAPS is described as follows.
BK ¼ ðBK Þl;m ¼ Φðj j xl xm j j Þ;
l ¼ 1; 2; …; M D ; m ¼ 1; 2; …; M D ; l ¼ 1; 2; …; M D ; m ¼ 1; 2; …; M D :
D
sðxi Þ ¼
sðxl ÞQ M ðξl ; ξi Þ; i ¼ 1; …; n;
scales at the computational cost of OðM nÞ. Overall, this algorithm pffiffiffi like Oðð2nM D Þn þ ðM 2D ÞÞ. If we choose M o o ðn=ð1 þ 2ÞÞ1=D , then the above mentioned algorithm is faster than Oðn2 Þ. In the next section, we give a description to accelerate the MAPS using this fast summation technique.
4.1.3. Step 2: iterative step Iterative methods converge in k steps with each step needing a matrix vector multiplication Oðn2 Þ. At each iterative step k, we follow the iterative algorithm as it is except the matrix vector multiplication of the matrix A and the k-dimensional updated column vector B, which is computed by using the FSM procedure in the vectorized form as follows: First, compute
4. Fast method of approximate particular solutions (FMAPS)
R B ¼ M A RM
The matrix A of the linear system Aα ¼ F, arising from the MAPS discretization in (5) and (6), can be viewed as the formulation of two block matrices; i.e., " # PA A¼ ; QA
and then compute
l¼1 D
PA ¼ ϕðj j xi xj j j Þ 1 r i r n ;1 r j r n ; Q A ¼ Φðj j xi xj j j Þ n i
i
þ 1 r i r n;1 r j r n
:
Hence, if we want to multiply A by any n-dimensional column vector B ¼ ½b1 ; b2 ; …; bn T , then this can be done by " # PA B ; ð18Þ AB¼ QA B and each row of those individual products PA B, Q A B can be considered like the summation defined in (14) which can be computed in an efficient way by using the FSM described in the previous section. If we use some iterative methods to solve the linear system instead of using the direct solvers like Gaussian elimination, then these fast products help us to solve the linear system faster. For
1
;
IM ¼ I0R ðIK MÞ; BM ¼ B0R ðBK MÞ; which approximates (18) as PA B IM i.e.,
where
D
and
Q A B BM ;
"
# IM : AB BM The undetermined coefficient α in (7) can be obtained by the iterative method. 4.1.4. Step 3: evaluation step Once we determine the unknown coefficient vector α from step 2, we use the given nt test points and RBF centers to evaluate the desired solution at the test points. This is also done by multiplying the corresponding evaluation matrix with the coefficient vector α using the similar FSM technique. As we can see, in this procedure we are not computing matrix A explicitly. The only matrix evaluation we have to do is for the
A.R. Lamichhane et al. / Engineering Analysis with Boundary Elements 64 (2016) 290–294
matrices described in the pre-computational step. Now, in the next section, we give some numerical examples, which validate our proposed method.
5. Numerical results Numerical experiments have been done by using MATLABⓒ on a desktop computer which has a 64-bit operating system and 16 GB memory with Intel(R) Core(TM) i7-4770K CPU @ 3.50 GHz processor. The Root Mean Square Error (RMSE) is used to measure the accuracy of the numerical results. The RMSE is defined as follows: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u X u 1 nt RMSE ¼ t ð19Þ ðu^ uj Þ2 ; nt j ¼ 1 j
293
Table 2 RMSE and CPU time for a large number of collocation points in the square domain using the FMAPS. n
RMSE
CPU time
2502
1.124E-04
1.65
5002
8.197E-05
11.04
7502
1.709E-04
35.00
10002
2.964E-04
56.90
1 0.5
where nt ; u^ j , and uj are the number of test points, approximate solution, and the exact solution, respectively. The special function exponential integral in the particular solution of the Gaussian [14] is evaluated using the rational approximation techniques [8] with the help of the boost library [17].
0 −0.5
Example 1. Consider the convection-diffusion-reaction equation in 2D as defined in (8) with the Dirichlet boundary condition u ¼ mðx; yÞ. We choose k ¼ 1, f ðx; yÞ ¼ y cos ðyÞ, gðx; yÞ ¼ sinhðxÞ, hðx; yÞ ¼ x2 þ y2 , and lðx; yÞ, mðx; yÞ are obtained by the analytic solution:
−1 −1
uðx; yÞ ¼ sin ðπ xÞcoshðyÞ cos ðπ xÞsinhðyÞ:
−0.5
0
0.5
1
Fig. 1. The profile of gear-shape domain
The above differential equation has been considered in Reference [5]. For the computational domain, we consider the standard unit square Ω ¼ ½0; 12 . To perform the test, we choose various numbers of collocation points and 1000 randomly distributed test points inside the domain to evaluate the RMSE errors. For the FSM, we employ 12 12 Chebyshev nodes and the shape parameter for the Gaussian RBF is chosen as 1. To show the efficiency of the proposed method, we test the MAPS using standard Gaussian elimination (GE), iterative method (GMRES), and the proposed FSM with GMRES. In Table 1, we show the proposed method is far more efficient than the direct Gaussian elimination and iterative method. As far as the accuracy is concerned, they all produce a similar accuracy. When the number of collocation points is increased, the computational cost using the two traditional linear solvers becomes too expensive. To illustrate the efficiency of the proposed method, we choose a huge number of collocation points up to one million as shown in Table 2. As we can see in the table, for the case of one million collocation points, the required CPU time is only 56.90 s which is considered to be extremely fast. Note that the tolerance of GMRES is set as 1E-05 in the above computation. Next, we consider the gear-shape domain as shown in Fig. 1. The parametric equation of the gear-shape curve can be written as
Ω ¼ ðx; yÞ : x ¼ rðtÞ cos t; y ¼ rðtÞ sin t; 0 r t r 2π Table 1 RMSE and CPU time using various numbers of the collocation points and solvers in the square domain. n
RMSE
GE
GMRES
FSM þ GMRES
402
6.163E-05
0.32
0.34
0.05
602
6.637E-05
1.79
1.50
0.13
802
3.740E-05
6.28
4.70
0.20
1002
2.751E-04
18.25
11.36
0.35
1202
1.813E-04
42.86
23.29
0.40
Table 3 RMSE and CPU time using various numbers of the collocation points for the gearshape domain. ðni; nbÞ
RMSE
CPU (s)
(5195, 1000) (14,552, 3000) (25,953, 5000) (58,623, 10,000) (104,405, 20,000) (163,297, 25,000) (654,541, 40,000)
8.860E-05 8.382E-05 8.851E-05 9.371E-05 8.907E-05 9.833E-05 1.591E-04
0.52 1.23 1.97 3.29 6.01 17.81 77.15
where 1 rðtÞ ¼ 1 þ 10 tanhð10 sin ð12tÞÞ:
In this case, we choose 969 randomly distributed test points inside the domain for the evaluation of the RMSE. Table 3 shows the RMSE and the CPU time for various combinations of interior and boundary points using the FMAPS. In this case, the tolerance of GMRES is again set as 1E-05. Table 3 validates the efficiency of the FMAPS. Despite a large number of collocation points are used in this example, the numerical results are two orders less accurate than the results obtained in Reference [5]. Note that the large number of collocation does not necessary produce more accurate results due to the ill-conditioning of the resultant matrix and round-off errors. Furthermore, to increase the efficiency, the iterate method GMRES has been adopted in our algorithm. To increase the accuracy, we need to reduce the tolerance of stopping criteria of GMRES which will inevitably slow the computation. There is a trade-off between accuracy and efficiency in the proposed fast algorithm.
294
A.R. Lamichhane et al. / Engineering Analysis with Boundary Elements 64 (2016) 290–294
Table 4 RMSE and CPU time for different sizes of the computational points in the unit cube. FSM þGMRES
n
RMSE
GE
GMRES
133
1.804E-05
0.31
0.64
1.09
153
3.156E-05
0.76
1.20
1.44
173
6.025E-05
1.77
1.57
1.47
193
3.898E-05
4.18
3.08
1.66
213
1.217E-04
8.40
4.05
1.52
233
1.494E-04
17.42
6.88
1.97
253
3.118E-04
33.46
11.08
2.41
273
2.204E-04
64.39
17.64
3.04
Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.enganabound. 2015.12.015.
n
RMSE
303
3.514E-05
7.04
403
5.253E-05
14.92
503
3.717E-05
29.96
603
5.219E-05
62.035
703
6.037E-05
105.15
CPU (s)
Example 2. Consider a three dimensional Poisson equation with a Dirichlet boundary condition in a unit cube ½0; 13 ; u ¼ gðx; y; zÞ;
The first author acknowledges the financial support from the Department of Civil Engineering, National Taiwan University to carry out this research work.
Appendix A. Supplementary data
Table 5 RMSE and CPU time for various numbers of the collocation points in the unit cube by FMAPS.
Δu ¼ f ðx; y; zÞ;
Acknowledgement
ðx; y; zÞ A Ω; ðx; y; zÞ A ∂Ω;
where Ω and ∂Ω are the interior and boundary of the unit cube, f and g are defined according to the exact solution uðx; y; zÞ ¼ sin ðπ xÞ sin ðπ yÞ sin ðπ zÞ: Table 4 shows the accuracy and CPU time using the three techniques including the FMAPS for the Poisson equation with Dirichlet boundary condition. The tolerance for the GMRES is set as in the two dimensional case and 123 Chebyshev nodes in the 3D have been used for the FSM. As shown in Table 4, the FMAPS is far more efficient than the other two traditional approaches. It is noteworthy that in Table 5, the FMAPS needs only 105:15 s to solve the given problem with 703 ¼ 343; 000 collocation points.
6. Conclusion In this paper we propose to couple the FSM with GMRES to speed up the computational efficiency of the MAPS. FMAPS has been proven to be a very effective alternative for solving the PDEs when a large number of computational points is required for both two and three dimensional problems. The numerical results show the proposed method is highly efficient. The difficulty of global RBF collocation methods such as the MAPS for solving large-scale elliptic PDEs has been alleviated. It would be interesting to compare the performance of the proposed approach which is a global method to localized RBF collocation methods for solving largescale problems. To our knowledge, the proposed algorithm is the first MPS fast method which is a RBF-based method for solving PDEs. Furthermore, in this paper we are focusing on the improvement of efficiency. To further improve the accuracy of FMAPS will be the subject of our future research.
References [1] Abbasbandy S, Ghehsareh HR, Alhuthali MS, Alsulami HH. Comparison of meshless local weak and strong forms based on particular solutions for a nonclassical 2-d diffusion model. Eng Anal Bound Elem 2014;39:121–8. [2] Bustamante CA, Power H, Florez W, Hang CY. The global approximate particular solution meshless method for two-dimensional linear elasticity problems. Int J Comput Math 2013;90(5):978–93. [3] Bustamante CA, Power H, Florez WF. A global meshless collocation particular solution method for solving the two-dimensional Navier–Stokes system of equations. Comput Math Appl 2013;65:1939–55. [4] Bustamante CA, Power H, Sua YH, Florez W. A global meshless collocation particular solution method (integrated radial basis function) for twodimensional Stokes flow problems. Appl Math Model 2013;37:4538–47. [5] Chen CS, Fan CM, Wen PH. The method of particular solutions for solving elliptic problems with variable coefficients. Int J Comput Methods 2011;8:545–59. [6] Chen CS, Fan CM, Wen PH. The method of particular solutions for solving certain partial differential equations. Numer Methods Partial Differ Equ 2012;28:506–22. [7] Cheng AH-D. Particular solutions of Laplacian, Helmholtz-type, and polyharmonic operators involving higher order radial basis functions. Eng Anal Bound Elem 2000;24:531–8. [8] Cody WJ, Thacher Jr. HC. Rational Chebyshev approximations for the exponential integral E1 ðxÞ. Math Comput 1968;22:641–9. [9] Fong W, Darve E. The black-box fast multipole method. J Comput Phys 2009;228:8712–25. [10] Fu Z, Chen W, Ling L. Method of approximate particular solutions for constantand variable-order fractional diffusion models. Eng Anal Bound Elem 2015;57:37–46. [11] Arfken G. Mathematical Methods for Physicists. Orlando, FL: Academic Press; 1985. [12] Jiang T, Li M, Chen CS. The method of particular solutions for solving inverse problems of a nonhomogeneous convection-diffusion equation with variable coefficients. Numer Heat Transf Part A 2012;61(5):338–52. [13] Kansa EJ. Multiquadrics – a scattered data approximation scheme with applications to computational fluid dynamics – II. Comput Math Appl 1990;19 (8/9):147–61. [14] Lamichhane AR, Chen. CS. The closed-form particular solutions for Laplace and biharmonic operators using a Gaussian function. Appl Math Lett 2015;46:50– 6. [15] Li M, Chen W, Tsai CH. A regularization method for the approximate particular solution of nonhomogeneous cauchy problems of elliptic partial differential equations with variable coefficients. Eng Anal Bound Elem 2012;36:274–80. [16] Muleshkov AS, Golberg MA, Chen CS. Particular solutions of Helmholtz-type operators using higher order polyharmonic splines. Comput Mech 1999;23:411–9. [17] BOOST c þ þ library, 〈http://www.boost.org〉. [18] Tsai CC. Automatic particular solutions of arbitrary high-order splines associated with polyharmonic and poly-Helmholtz equations. Eng Anal Bound Elem 2011;35(7):925–34. [19] Tsai CC. Generalized polyharmonic multiquadrics. Eng Anal Bound Elem 2015;50(0):239–48. [20] Tsai CC, Cheng AH-D, Chen CS. Particular solutions of splines and monomials for polyharmonic and products of Helmholtz operators. Eng Anal Bound Elem 2009;33:514–21.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Fast method of approximate particular solutions using Chebyshev interpolation A.R. Lamichhane a, D.L. Young b, C.S. Chen a,c,n a
Department of Mathematics, University of Southern Mississippi, USA Department of Civil Engineering, National Taiwan University, Taiwan c College of Mathematics, Taiyuan University of Technology, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 30 October 2015 Received in revised form 20 December 2015 Accepted 21 December 2015 Available online 14 January 2016
The fast method of approximate particular solutions (FMAPS) is based on the global version of the method of approximate particular solutions (MAPS). In this method, given partial differential equations are discretized by the usual MAPS and the determination of the unknown coefficients is accelerated using a fast technique. Numerical results confirm the efficiency of the proposed technique for the PDEs with a large number of computational points in both two and three dimensions. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Radial basis functions Meshless method Method of approximate particular solutions Chebyshev interpolation Fast summation method Fast method of approximate particular solutions
1. Introduction During the past two decades, the radial basis functions (RBFs) have been widely applied for solving partial differential equations (PDEs). Kansa [13] first proposed the so-called RBF collocation method in 1990 for computationally solving fluid dynamic problems. One of the attractions of the Kansa method is its simplicity for solving problems in high dimensional and irregular domains. Due to the popularity of the Kansa method, several other RBF collocation methods have been proposed. Among them, the method of approximate particular solutions (MAPS) [5,6] is another effective method using the particular solution of the given RBFs as the basis function. The formulation of the MAPS produces a full and dense matrix system which is often very ill-conditioned. Traditionally, this matrix system is solved by using direct or iterative methods. Direct methods such as Gaussian elimination require OðN3 Þ operations for an N N system of equations. For iterative methods, we may obtain the approximate solution in k steps with each step needing a matrix vector multiplication OðN 2 Þ. In recent years, the MAPS has been widely applied to solve physical and engineering problems like Navier–Stokes equations [3], Stokes flow problems [4], Linear elasticity equations [2], Timen
Corresponding author. E-mail address: [email protected] (C.S. Chen).
http://dx.doi.org/10.1016/j.enganabound.2015.12.015 0955-7997/& 2016 Elsevier Ltd. All rights reserved.
fractional diffusion equations [10], Inverse problem of nonhomogeneous convection-diffusion equation [12], Diffusion equation with non-classical boundary [1] and Nonhomogeneous cauchy problems of elliptic PDEs [15]. Simulations of these kinds of problems involve a large number of interpolation points. The high computational cost using traditional solvers has become an issue. In this paper we pay special attention on how to develop a fast algorithm to alleviate the issue of high cost for solving large-scale problems using the MAPS. Consequently, we propose to couple the MAPS with a fast summation method to reduce the computational time by multiplying a matrix and a vector in each step inside the iterative method. This fast summation method is based on the Chebyshev interpolation [9]. As we shall see in the section of numerical results, we have successfully solved a 2D problem using 694,541 collocation points with only 77.15 s of computer running time and 343,000 collocation points with 105.15 s for a 3D problem. Moreover, we do not compromise the accuracy for our proposed fast computation. The structure of the paper is as follows. In Section 2, we give a brief review of the MAPS and the closed-form particular solution of the Gaussian for Laplacian in 2D and 3D. In Section 3, we review the algorithm of fast summation method (FSM). In Section 4, we propose the fast method of approximate particular solutions (FMAPS) by coupling the FSM and the MAPS as a fast algorithm for solving PDEs which require a large number of collocation points. A
A.R. Lamichhane et al. / Engineering Analysis with Boundary Elements 64 (2016) 290–294
specific algorithm has been given. In Section 4.1, to demonstrate the efficiency of the proposed method, two numerical examples in 2D and 3D are given. Finally, a short conclusion is drawn in Section 4.1.1.
2. Method of approximate particular solutions (MAPS)
ΔuðxÞ ¼ f ðxÞ; x A Ω;
ð1Þ
uðxÞ ¼ gðxÞ; x A ∂Ω;
ð2Þ
where Δ is the Laplacian operator, Ω and ∂Ω are the interior and boundary of the computational domain, respectively. Suppose f xi gni¼ 1 are the interpolation points containing ni interior points in Ω and nb boundary points on ∂Ω; i.e., n ¼ ni þnb . Let ϕ be a given radial basis function. By the MAPS [6], we assume the solution to (1) and (2) can be approximated by ^ uðxÞ uðxÞ ¼
αj Φð J x xj J Þ;
ð3Þ
where J J is the Euclidean norm, fαj g are the undetermined coefficients, and
ΔΦ ¼ ϕ:
ð4Þ
By the collocation method, from (1) and (2), we have
αj ϕ J xi xj J ¼ f ðxi Þ;
1 r i rni ;
ð5Þ
ni þ 1 r i r n:
ð6Þ
j¼1 n X
αj Φ J xi xj J ¼ gðxi Þ;
Z
1
EiðxÞ ¼ x
eu du; u
ð11Þ
and γ C 0:5772156649015328 is the Euler–Mascheroni constant [11]. Note that EiðxÞ is the special function known as the exponential integral function [11]. Theorem 2. Let ϕðrÞ ¼ expð cr 2 Þ, and ΔΦðrÞ ¼ ϕðrÞ in 3D. Then, 8 pffiffiffiffi pffiffiffi π > > < 3=2 erf ð crÞ; r a 0; 4c r ð12Þ ΦðrÞ ¼ 1 > > : ; r ¼ 0; 2c where erfðxÞ is the special function defined as follows [11]: Z x 2 2 erfðxÞ ¼ pffiffiffiffi e u du:
π
j¼1
n X
Theorem 1. Let ϕðrÞ ¼ expð cr 2 Þ, c 4 0, and ΔΦðrÞ ¼ ϕðrÞ in 2D. Then, 8 1 1 > 2 > < Eiðcr Þ þ log ðrÞ; r a0; 4c 2c ð10Þ ΦðrÞ ¼ 1 > > ðγ þ log ðcÞÞ; r ¼ 0; : 4c where
For simplicity, let us consider the following Poisson equation with Dirichlet boundary condition
n X
291
ð13Þ
0
In this paper, we adopt the Gaussian, ϕ, as the RBF basis function and the corresponding particular solutions, Φ, as the basis functions for the approximation of the partial differential equation. The particular solutions of the Gaussian in (10) and (12) contain the special functions, EiðxÞ and erf ðxÞ, which is costly in terms of numerical evaluation. The efficiency can be significantly improved using compiled MATLAB MEX functions. In the next section, we briefly introduce a fast summation method for the matrix vector multiplication.
j¼1
From (5) and (6), we can formulate a linear system of equations Aα ¼ F
Consider the evaluation of the sum of the form
where " #
sðxÞ ¼
ϕ A¼ Φ
ni þ1 r k r n;
; ij
Φ ¼ Φ J xk xj J
; kj
1 ri r ni ; 1 r j rn;
α ¼ ½α1 α2 ⋯ αn T ;
F ¼ ½f ðx1 Þ ⋯ f ðxni Þ gðxni þ 1 Þ ⋯ gðxn ÞT . For a more general form of the PDEs such as the following convection-diffusion-reaction equation in 2D, ∂u ∂u þ gðxÞ þ hðxÞu ¼ lðxÞ; x ¼ ðx; yÞ A Ω; ∂x ∂y
ð8Þ
bj κ j j x xj j j ;
ð14Þ
αj ðkϕð J x x J Þ þ f ðxÞΦx J x xj J þ gðxÞΦy J x xj J
P M ðξ; ηÞ ¼
1 1 2 MX þ T ðξÞT i ðηÞ; M M i¼1 i
ð15Þ
MD
j¼1
þ hðxÞΦ J x xj J Þ ¼ lðxÞ:
where κ is either the RBFs or the particular solution of the corresponding RBFs. We can evaluate the sum (14) in an efficient way by using the Chebyshev interpolation technique as described in [9]. From [9], let
where ξ; η A ½ 1; 1; T i is the first kind Chebyshev polynomial of order i. Let H be a hypercube in D dimension which contains all the
we have [5] N X
n X j¼1
ϕ ¼ ϕ J xi xj J
kΔu þ f ðxÞ
3. Fast summation method (FSM)
ð7Þ
ð9Þ
One of the key procedures in the implementation of the MAPS is to obtain the closed-form expression for the particular solutions of the corresponding RBFs. The derivation of the particular solutions for the well-known RBFs has already been known [7,16,18– 20]. Recently, the closed-form particular solutions of the Gaussian have been derived in [14], which are stated as follows:
given collocation points fxi gni¼ 1 . Consider fξl gl ¼ 1 , fηl gM l ¼ 1 be two identical sets of Chebyshev nodes in ½ 1; 1D . Then by using linear transformations, we can map xi ; xj into ξi ; ηj , and ξl ; ηm in ½ 1; 1D into xl ; xm in H, respectively. Instead of directly evaluating (14), we approximate it by using the following functional approximation XX κ ð j j x yj j Þ ¼ k j j xl ym j j Q M ðξ l ; ξÞQ M ðηm ; ηÞ; ð16Þ l
D
m
where ξ; η are the linear transformations of x; y, respectively, and
292
A.R. Lamichhane et al. / Engineering Analysis with Boundary Elements 64 (2016) 290–294
in 2D, Q M ðξ; ηÞ ¼ P M ðξ1 ; η1 ÞP M ðξ2 ; η2 Þ;
ξ ¼ ðξ1 ; ξ2 Þ; η ¼ ðη1 ; η2 Þ;
and in 3D, Q M ðξ; ηÞ ¼ P M ðξ1 ; η1 ÞP M ðξ2 ; η2 ÞP M ðξ3 ; η3 Þ; η ¼ ðη1 ; η2 ; η3 Þ:
ξ ¼ ðξ1 ; ξ2 ; ξ3 Þ;
Using (16), (14) can be separated into the product of the different sums as follows: n X
sðxi Þ ¼
n X bj κ j j xi xj j j bj
j¼1
j¼1
M M X X D
4.1. Algorithm
!
D
κ ð j j xl xm j j ÞQ M ðξl ; ξi ÞQ M ðηm ; ηj Þ
l¼1m¼1
¼
MD X l¼1
Q M ðξ l ; ξ i Þ
MD X
κ ð j j xl xm j j Þ
m¼1
n X
bj Q M ðηm ; ηj Þ:
ð17Þ
Wm ¼
bj Q M ðηm ; ηj Þ; m ¼ 1; …; M ; D
and the computational cost for these Chebyshev weights is OðM D nÞ. Then, we compute the sðxÞ at the Chebyshev nodes xl as M X
Find a hypercube H in D dimension which contains fxi gni¼ 1 by using linear transformations, find ξi ; ηj in ½ 1; 1D corre-
sponding to xi ; xj , respectively, find xl ; xm in H corresponding to ξl ; ηm , respectively, compute the following matrices: R ¼ ðRÞm;j ¼ Q M ðηm ; ηj Þ;
IR ¼ ðI R Þl;i ¼ Q M ðξl ; ξi Þ;
m ¼ 1; 2; …; M D ; j ¼ 1; 2; …; n; l ¼ 1; 2; …; M D ; i ¼ 1; 2; …; ni ;
BR ¼ ðBR Þl;i ¼ Q M ðξl ; ξi Þ;
D
W m κ ð j j xl xm j j Þ; l ¼ 1; …; M D ;
l ¼ 1; 2; …; M D ; i ¼ ni þ1; ni þ 2; …; n;
IK ¼ ðI K Þl;m ¼ ϕðj j xl xm j j Þ;
m¼1
with computational cost OðM 2D Þ. Eventually, we compute the sðxÞ at the collocation points xi as M X
4.1.2. Step 1: pre-computational step
j¼1
sðxl Þ ¼
4.1.1. Input Collocation points fxi gni¼ 1 , Chebyshev nodes ξl ; ηm , dimension of the domain D.
j¼1
If we choose those M D Chebyshev nodes such that M D o o n, then the summation can be computed efficiently. To evaluate (14) approximately, let us first evaluate the last term of (17) which gives us the weights at the Chebyshev nodes xm n X
this, we can use any iterative method which includes the vector multiplication with matrix A. In this paper, we have used GMRES iterative method to find the unknown coefficients. As we know that, in the GMRES, at each iteration, we have to multiply a matrix and an updated vector. This matrix and vector multiplication is obtained faster by the technique that we have described earlier. This kind of approach which is called FMAPS, actually solves the PDEs by using MAPS without explicitly computing matrix A which requires a minimum storage. The computational procedure for the FMAPS is described as follows.
BK ¼ ðBK Þl;m ¼ Φðj j xl xm j j Þ;
l ¼ 1; 2; …; M D ; m ¼ 1; 2; …; M D ; l ¼ 1; 2; …; M D ; m ¼ 1; 2; …; M D :
D
sðxi Þ ¼
sðxl ÞQ M ðξl ; ξi Þ; i ¼ 1; …; n;
scales at the computational cost of OðM nÞ. Overall, this algorithm pffiffiffi like Oðð2nM D Þn þ ðM 2D ÞÞ. If we choose M o o ðn=ð1 þ 2ÞÞ1=D , then the above mentioned algorithm is faster than Oðn2 Þ. In the next section, we give a description to accelerate the MAPS using this fast summation technique.
4.1.3. Step 2: iterative step Iterative methods converge in k steps with each step needing a matrix vector multiplication Oðn2 Þ. At each iterative step k, we follow the iterative algorithm as it is except the matrix vector multiplication of the matrix A and the k-dimensional updated column vector B, which is computed by using the FSM procedure in the vectorized form as follows: First, compute
4. Fast method of approximate particular solutions (FMAPS)
R B ¼ M A RM
The matrix A of the linear system Aα ¼ F, arising from the MAPS discretization in (5) and (6), can be viewed as the formulation of two block matrices; i.e., " # PA A¼ ; QA
and then compute
l¼1 D
PA ¼ ϕðj j xi xj j j Þ 1 r i r n ;1 r j r n ; Q A ¼ Φðj j xi xj j j Þ n i
i
þ 1 r i r n;1 r j r n
:
Hence, if we want to multiply A by any n-dimensional column vector B ¼ ½b1 ; b2 ; …; bn T , then this can be done by " # PA B ; ð18Þ AB¼ QA B and each row of those individual products PA B, Q A B can be considered like the summation defined in (14) which can be computed in an efficient way by using the FSM described in the previous section. If we use some iterative methods to solve the linear system instead of using the direct solvers like Gaussian elimination, then these fast products help us to solve the linear system faster. For
1
;
IM ¼ I0R ðIK MÞ; BM ¼ B0R ðBK MÞ; which approximates (18) as PA B IM i.e.,
where
D
and
Q A B BM ;
"
# IM : AB BM The undetermined coefficient α in (7) can be obtained by the iterative method. 4.1.4. Step 3: evaluation step Once we determine the unknown coefficient vector α from step 2, we use the given nt test points and RBF centers to evaluate the desired solution at the test points. This is also done by multiplying the corresponding evaluation matrix with the coefficient vector α using the similar FSM technique. As we can see, in this procedure we are not computing matrix A explicitly. The only matrix evaluation we have to do is for the
A.R. Lamichhane et al. / Engineering Analysis with Boundary Elements 64 (2016) 290–294
matrices described in the pre-computational step. Now, in the next section, we give some numerical examples, which validate our proposed method.
5. Numerical results Numerical experiments have been done by using MATLABⓒ on a desktop computer which has a 64-bit operating system and 16 GB memory with Intel(R) Core(TM) i7-4770K CPU @ 3.50 GHz processor. The Root Mean Square Error (RMSE) is used to measure the accuracy of the numerical results. The RMSE is defined as follows: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u X u 1 nt RMSE ¼ t ð19Þ ðu^ uj Þ2 ; nt j ¼ 1 j
293
Table 2 RMSE and CPU time for a large number of collocation points in the square domain using the FMAPS. n
RMSE
CPU time
2502
1.124E-04
1.65
5002
8.197E-05
11.04
7502
1.709E-04
35.00
10002
2.964E-04
56.90
1 0.5
where nt ; u^ j , and uj are the number of test points, approximate solution, and the exact solution, respectively. The special function exponential integral in the particular solution of the Gaussian [14] is evaluated using the rational approximation techniques [8] with the help of the boost library [17].
0 −0.5
Example 1. Consider the convection-diffusion-reaction equation in 2D as defined in (8) with the Dirichlet boundary condition u ¼ mðx; yÞ. We choose k ¼ 1, f ðx; yÞ ¼ y cos ðyÞ, gðx; yÞ ¼ sinhðxÞ, hðx; yÞ ¼ x2 þ y2 , and lðx; yÞ, mðx; yÞ are obtained by the analytic solution:
−1 −1
uðx; yÞ ¼ sin ðπ xÞcoshðyÞ cos ðπ xÞsinhðyÞ:
−0.5
0
0.5
1
Fig. 1. The profile of gear-shape domain
The above differential equation has been considered in Reference [5]. For the computational domain, we consider the standard unit square Ω ¼ ½0; 12 . To perform the test, we choose various numbers of collocation points and 1000 randomly distributed test points inside the domain to evaluate the RMSE errors. For the FSM, we employ 12 12 Chebyshev nodes and the shape parameter for the Gaussian RBF is chosen as 1. To show the efficiency of the proposed method, we test the MAPS using standard Gaussian elimination (GE), iterative method (GMRES), and the proposed FSM with GMRES. In Table 1, we show the proposed method is far more efficient than the direct Gaussian elimination and iterative method. As far as the accuracy is concerned, they all produce a similar accuracy. When the number of collocation points is increased, the computational cost using the two traditional linear solvers becomes too expensive. To illustrate the efficiency of the proposed method, we choose a huge number of collocation points up to one million as shown in Table 2. As we can see in the table, for the case of one million collocation points, the required CPU time is only 56.90 s which is considered to be extremely fast. Note that the tolerance of GMRES is set as 1E-05 in the above computation. Next, we consider the gear-shape domain as shown in Fig. 1. The parametric equation of the gear-shape curve can be written as
Ω ¼ ðx; yÞ : x ¼ rðtÞ cos t; y ¼ rðtÞ sin t; 0 r t r 2π Table 1 RMSE and CPU time using various numbers of the collocation points and solvers in the square domain. n
RMSE
GE
GMRES
FSM þ GMRES
402
6.163E-05
0.32
0.34
0.05
602
6.637E-05
1.79
1.50
0.13
802
3.740E-05
6.28
4.70
0.20
1002
2.751E-04
18.25
11.36
0.35
1202
1.813E-04
42.86
23.29
0.40
Table 3 RMSE and CPU time using various numbers of the collocation points for the gearshape domain. ðni; nbÞ
RMSE
CPU (s)
(5195, 1000) (14,552, 3000) (25,953, 5000) (58,623, 10,000) (104,405, 20,000) (163,297, 25,000) (654,541, 40,000)
8.860E-05 8.382E-05 8.851E-05 9.371E-05 8.907E-05 9.833E-05 1.591E-04
0.52 1.23 1.97 3.29 6.01 17.81 77.15
where 1 rðtÞ ¼ 1 þ 10 tanhð10 sin ð12tÞÞ:
In this case, we choose 969 randomly distributed test points inside the domain for the evaluation of the RMSE. Table 3 shows the RMSE and the CPU time for various combinations of interior and boundary points using the FMAPS. In this case, the tolerance of GMRES is again set as 1E-05. Table 3 validates the efficiency of the FMAPS. Despite a large number of collocation points are used in this example, the numerical results are two orders less accurate than the results obtained in Reference [5]. Note that the large number of collocation does not necessary produce more accurate results due to the ill-conditioning of the resultant matrix and round-off errors. Furthermore, to increase the efficiency, the iterate method GMRES has been adopted in our algorithm. To increase the accuracy, we need to reduce the tolerance of stopping criteria of GMRES which will inevitably slow the computation. There is a trade-off between accuracy and efficiency in the proposed fast algorithm.
294
A.R. Lamichhane et al. / Engineering Analysis with Boundary Elements 64 (2016) 290–294
Table 4 RMSE and CPU time for different sizes of the computational points in the unit cube. FSM þGMRES
n
RMSE
GE
GMRES
133
1.804E-05
0.31
0.64
1.09
153
3.156E-05
0.76
1.20
1.44
173
6.025E-05
1.77
1.57
1.47
193
3.898E-05
4.18
3.08
1.66
213
1.217E-04
8.40
4.05
1.52
233
1.494E-04
17.42
6.88
1.97
253
3.118E-04
33.46
11.08
2.41
273
2.204E-04
64.39
17.64
3.04
Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.enganabound. 2015.12.015.
n
RMSE
303
3.514E-05
7.04
403
5.253E-05
14.92
503
3.717E-05
29.96
603
5.219E-05
62.035
703
6.037E-05
105.15
CPU (s)
Example 2. Consider a three dimensional Poisson equation with a Dirichlet boundary condition in a unit cube ½0; 13 ; u ¼ gðx; y; zÞ;
The first author acknowledges the financial support from the Department of Civil Engineering, National Taiwan University to carry out this research work.
Appendix A. Supplementary data
Table 5 RMSE and CPU time for various numbers of the collocation points in the unit cube by FMAPS.
Δu ¼ f ðx; y; zÞ;
Acknowledgement
ðx; y; zÞ A Ω; ðx; y; zÞ A ∂Ω;
where Ω and ∂Ω are the interior and boundary of the unit cube, f and g are defined according to the exact solution uðx; y; zÞ ¼ sin ðπ xÞ sin ðπ yÞ sin ðπ zÞ: Table 4 shows the accuracy and CPU time using the three techniques including the FMAPS for the Poisson equation with Dirichlet boundary condition. The tolerance for the GMRES is set as in the two dimensional case and 123 Chebyshev nodes in the 3D have been used for the FSM. As shown in Table 4, the FMAPS is far more efficient than the other two traditional approaches. It is noteworthy that in Table 5, the FMAPS needs only 105:15 s to solve the given problem with 703 ¼ 343; 000 collocation points.
6. Conclusion In this paper we propose to couple the FSM with GMRES to speed up the computational efficiency of the MAPS. FMAPS has been proven to be a very effective alternative for solving the PDEs when a large number of computational points is required for both two and three dimensional problems. The numerical results show the proposed method is highly efficient. The difficulty of global RBF collocation methods such as the MAPS for solving large-scale elliptic PDEs has been alleviated. It would be interesting to compare the performance of the proposed approach which is a global method to localized RBF collocation methods for solving largescale problems. To our knowledge, the proposed algorithm is the first MPS fast method which is a RBF-based method for solving PDEs. Furthermore, in this paper we are focusing on the improvement of efficiency. To further improve the accuracy of FMAPS will be the subject of our future research.
References [1] Abbasbandy S, Ghehsareh HR, Alhuthali MS, Alsulami HH. Comparison of meshless local weak and strong forms based on particular solutions for a nonclassical 2-d diffusion model. Eng Anal Bound Elem 2014;39:121–8. [2] Bustamante CA, Power H, Florez W, Hang CY. The global approximate particular solution meshless method for two-dimensional linear elasticity problems. Int J Comput Math 2013;90(5):978–93. [3] Bustamante CA, Power H, Florez WF. A global meshless collocation particular solution method for solving the two-dimensional Navier–Stokes system of equations. Comput Math Appl 2013;65:1939–55. [4] Bustamante CA, Power H, Sua YH, Florez W. A global meshless collocation particular solution method (integrated radial basis function) for twodimensional Stokes flow problems. Appl Math Model 2013;37:4538–47. [5] Chen CS, Fan CM, Wen PH. The method of particular solutions for solving elliptic problems with variable coefficients. Int J Comput Methods 2011;8:545–59. [6] Chen CS, Fan CM, Wen PH. The method of particular solutions for solving certain partial differential equations. Numer Methods Partial Differ Equ 2012;28:506–22. [7] Cheng AH-D. Particular solutions of Laplacian, Helmholtz-type, and polyharmonic operators involving higher order radial basis functions. Eng Anal Bound Elem 2000;24:531–8. [8] Cody WJ, Thacher Jr. HC. Rational Chebyshev approximations for the exponential integral E1 ðxÞ. Math Comput 1968;22:641–9. [9] Fong W, Darve E. The black-box fast multipole method. J Comput Phys 2009;228:8712–25. [10] Fu Z, Chen W, Ling L. Method of approximate particular solutions for constantand variable-order fractional diffusion models. Eng Anal Bound Elem 2015;57:37–46. [11] Arfken G. Mathematical Methods for Physicists. Orlando, FL: Academic Press; 1985. [12] Jiang T, Li M, Chen CS. The method of particular solutions for solving inverse problems of a nonhomogeneous convection-diffusion equation with variable coefficients. Numer Heat Transf Part A 2012;61(5):338–52. [13] Kansa EJ. Multiquadrics – a scattered data approximation scheme with applications to computational fluid dynamics – II. Comput Math Appl 1990;19 (8/9):147–61. [14] Lamichhane AR, Chen. CS. The closed-form particular solutions for Laplace and biharmonic operators using a Gaussian function. Appl Math Lett 2015;46:50– 6. [15] Li M, Chen W, Tsai CH. A regularization method for the approximate particular solution of nonhomogeneous cauchy problems of elliptic partial differential equations with variable coefficients. Eng Anal Bound Elem 2012;36:274–80. [16] Muleshkov AS, Golberg MA, Chen CS. Particular solutions of Helmholtz-type operators using higher order polyharmonic splines. Comput Mech 1999;23:411–9. [17] BOOST c þ þ library, 〈http://www.boost.org〉. [18] Tsai CC. Automatic particular solutions of arbitrary high-order splines associated with polyharmonic and poly-Helmholtz equations. Eng Anal Bound Elem 2011;35(7):925–34. [19] Tsai CC. Generalized polyharmonic multiquadrics. Eng Anal Bound Elem 2015;50(0):239–48. [20] Tsai CC, Cheng AH-D, Chen CS. Particular solutions of splines and monomials for polyharmonic and products of Helmholtz operators. Eng Anal Bound Elem 2009;33:514–21.