A gradient reproducing kernel collocation method for boundary value problems

全文

(1)INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng 2013; 93:1381–1402 Published online 31 January 2013 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/nme.4432. A gradient reproducing kernel collocation method for boundary value problems Sheng-Wei Chi1 , Jiun-Shyan Chen2, * ,† , Hsin-Yun Hu3 and Judy P. Yang4 1 Civil 2 Civil. and Materials Engineering Department, University of Illinois at Chicago, IL 60607, USA and Environmental Engineering Department, University of California, Los Angeles, CA 90095, USA 3 Applied Mathematics Department, Tunghai University, Taiwan, R.O.C. 4 Civil Engineering Department, National Chiao Tung University, Taiwan, R.O.C.. SUMMARY The earlier work in the development of direct strong form collocation methods, such as the reproducing kernel collocation method (RKCM), addressed the domain integration issue in the Galerkin type meshfree method, such as the reproducing kernel particle method, but with increased computational complexity because of taking higher order derivatives of the approximation functions and the need for using a large number of collocation points for optimal convergence. In this work, we intend to address the computational complexity in RKCM while achieving optimal convergence by introducing a gradient reproduction kernel approximation. The proposed gradient RKCM reduces the order of differentiation to the first order for solving second-order PDEs with strong form collocation. We also show that, different from the typical strong form collocation method where a significantly large number of collocation points than the number of source points is needed for optimal convergence, the same number of collocation points and source points can be used in gradient RKCM. We also show that the same order of convergence rates in the primary unknown and its first-order derivative is achieved, owing to the imposition of gradient reproducing conditions. The numerical examples are given to verify the analytical prediction. Copyright © 2012 John Wiley & Sons, Ltd. Received 5 May 2012; Revised 21 August 2012; Accepted 17 September 2012 KEY WORDS:. reproducing kernel collocation method; gradient reproducing kernel approximation; weighted collocation method; strong form collocation. 1. INTRODUCTION In the past two decades, significant advancement has been achieved in the development of meshfree methods for solving PDEs based on the Galerkin weak formulation. The approximation functions with compact support such as moving least-squares (MLS) [1–3] and reproducing kernel (RK) [4–6] functions are commonly adopted in Galerkin meshfree methods. With monomial reproducing properties in compactly supported MLS and RK, algebraic convergence rates are obtained [7, 8] and the discrete systems are well-conditioned. Nonetheless, domain integration of the weak equation adds substantial difficulties and complexities to the Galerkin meshfree methods [9–13]. On the other hand, meshfree methods formulated based on the strong form with direct collocation have also been proposed [14–20]. This approach reduces the complexities associated with domain integration and the imposition of boundary conditions. The radial basis functions (RBFs) [21–24] are commonly used in the strong form collocation method [15–17], generally called the radial basis collocation method (RBCM). While the nonlocal RBFs with certain regularity offer exponential. *Correspondence to: Jiun-Shyan Chen, Civil and Environmental Engineering Department, University of California, Los Angeles, CA 90095, USA. † E-mail: [email protected] Copyright © 2012 John Wiley & Sons, Ltd..

(2) 1382. S.-W. CHI ET AL.. convergence in RBCM [25–27], the linear system of RBCM is typically ill-conditioned [28, 29]. An alternative approach is the employment of smooth approximation with compact support such as the MLS or RK approximation in the strong form collocation method [14, 18, 19, 30, 31]. The reproducing kernel collocation method (RKCM) offers a much better conditioned discrete system than that of RBCM; nevertheless, it converges algebraically [30, 31]. The work in [32] shows that one can construct a localized RBF using a partition of unity function, such as the reproducing kernel enhanced radial basis function, to yield a local approximation while maintaining the exponential convergence in RBCM. This localized RBF, combined with the subdomain collocation method, has been applied to problems with local features, such as problems with heterogeneity [33] or cracks [34] that are difficult to be solved by RBCM. It is noteworthy that higher order derivatives of the approximation functions are needed in the strong form collocation method compared with the Galerkin method. While approximation functions such as RK and MLS can be arbitrarily smooth, taking derivatives of these functions is computationally costly, making RKCM less efficient. In particular, the high complexity in RKCM is caused by taking derivatives of the moment matrix inversion in the multidimensional RK shape functions (see the detailed complexity and error analysis of RKCM in [31] and [30], respectively). Furthermore, for optimal convergence in RBCM and RKCM, using the number of collocation points much larger than the number of source points is needed, and this adds additional computational effort [15, 30]. Motivated by the above mentioned disadvantages in RKCM, a gradient RK approximation is introduced in solving second-order PDEs with strong form collocation, termed the gradient reproducing kernel collocation method (G-RKCM). The idea of gradient RK was first introduced in the Galerkin weak form to achieve a synchronized convergence [35, 36]. The gradient RK approximation in [35, 36] is formulated based on partition of nullity and derivative reproducing conditions, where similar construction has also been introduced in the implicit gradient approximation for localization problems [37]. Different from [35–38] where the gradient RK approximation is used as the enrichment of the standard RK approximation under Galerkin weak formulation, the present approach introduces gradient RK as the ‘assumed strain’ field directly in the strong form. The convergence properties of this G-RKCM approach will be derived, and the complexity of this method in comparison with RKCM will also be analyzed in this paper. The paper is organized as follows. Section 2 reviews the basic equations and the fundamental properties of RK approximation and RKCM. In Section 3, the gradient RK approximation is introduced, and its application to the strong form to construct G-RKCM discrete equations is presented in Section 4. The error analysis of G-RKCM and the choice of collocation points are given in Section 5. The complexities of G-RKCM and RKCM are compared in Section 6. The numerical examples are given in Section 7 to demonstrate the effectiveness of the proposed method. The concluding remarks of the proposed G-RKCM are presented in Section 8. 2. REVIEW OF REPRODUCING KERNEL COLLOCATION METHOD Consider the following boundary value problem: Lu D f in Bh u D h on @h Bg u D g on @g ,. (1). where is the problem domain, @h is the Neumann boundary, @g is the Dirichlet boundary, @ D @h [ @g , L is the differential operator in , and Bh and Bg are the boundary operators on @h and @g , respectively. To solve (1) by strong form collocation, the reproducing kernel approximation of u, denoted by v, is expressed as. u.x/ v.x/ D. Ns X. ‰I .x/aI ,. (2). I D1. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(3) GRADIENT REPRODUCING KERNEL COLLOCATION METHOD. 1383. where Ns is the number of source points, and ‰I .x/ is the reproducing kernel (RK) shape function expressed as ‰I .x/ D C.xI x xI /'a .x xI / ,. (3). where 'a .x xI / is the kernel function, and C .xI x xI / is the correction function C .xI x xI / D. p X. b˛ .x/.x xI /˛ , p > 0. (4). j˛jD0. DW HT .x xI /b.x/. Here we introduce the multi-index notation in d -dimension ˛ D .˛1 , ˛2 , : : : , ˛d /, with the length Pd ˛d ˛d ˛1 ˛2 ˛1 ˛2 ˛ ˛ ˛ of ˛ defined as j˛j D i D1 ˛i , x x1 x2 xd , xI x1I x2I xdI , .x xI / .x1 x1I /˛1 .x2 x2I /˛2 .xd xdI /˛d , and b˛ b˛1 ,˛2 ,:::,˛d . The vectors HT .x xI / and bT .x/ are the corresponding row vectors of ¹.x xI /˛ ºj˛j6p and ¹b˛ .x/ºj˛j6p , respectively. The shape functions are required to satisfy p-th order reproducing conditions given as follows: X I .x/ x˛I D x˛ , j˛j 6 p. (5) I. The coefficients b.x/ are obtained by satisfying (5), and it yields the following RK shape function: ‰I .x/ D HT .0/M1 .x/H .x xI / 'a .x xI /. (6). and M.x/ D. Ns X. H.x xI /HT .x xI /'a .x xI / .. (7). I D1. Introducing RK approximation of u in (2) to the strong form in (1), and evaluating the differential equation and boundary conditions at the collocation points p` 2 , q` 2 @h , and r` 2 @g , we have the following collocation equations: Lv.p` / D f .p` / 8 p` 2 , ` D 1, , Np Bh v .q` / D h .q` / 8 q` 2 @h , ` D 1, , Nq Bg v .r` / D g .r` / 8 r` 2 @g , ` D 1, , Nr .. (8). Collection of the collocation equations yields the following linear system: . . Aa D b,. . . (9). where A D A .L‰/p` , .Bh ‰/q` , .Bg ‰/r` and b D b .f/p` , .h/q` , .g/r` 8 p` 2 , q` 2 @h , and r` 2 @g . Note that the total number of collocation points Np C Nq C Nr is typically much larger than the number of source points Ns for optimal convergence, and hence yields an over-determined system in (9). Remark 2.1 The collocation equations in (8) can be shown to be equivalent to the minimization of the following least-squares functional with quadrature [15], that is, to seek solution ur 2 V D span ¹‰1 , , ‰Ns º, such that E .ur / D inf E .v/ , v2V. where E.v/ D. Z 1 .Lv f/T .Lv f/ d C .Bh v h/T .Bh v h/ d 2 @h Z T 1 Bg v g Bg v g d C 2 @g 1 2. (10). Z. Copyright © 2012 John Wiley & Sons, Ltd.. (11). Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(4) 1384. S.-W. CHI ET AL.. By choosing the quadrature points in (11) the same as the collocation points in (8) in solving (9) by a weighted least-squares method, the equivalence between the solution by minimization of (11) and the solution of (9) can be established; see [15] for details. Remark 2.2 To keep the balance of errors in the domain and boundary terms in the least-squares functional, a weighted least-squares functional has been proposed [15, 30] Z Z 1 ˛h T .Lv f/ .Lv f/ d C .Bh v h/T .Bh v h/ d E.v/ D 2 2 @h Z (12) T ˛g C Bg v g Bg v g d, 2 @g p p where the weights ˛h D 1, ˛g D Ns , with D 1 for Poisson problem and D max¹, º for elasticity for optimal convergence have been proposed. A set of equivalent collocation equations can be obtained Lv .p` / D f .p` / ` D 1, , Np 8 p` 2 , p p ˛h Bh v .q` / D ˛h h .q` / 8 q` 2 @h , ` D 1, , Nq p p ˛g Bg v .r` / D ˛g g .r` / 8 r` 2 @g , ` D 1, , Nr . This RKCM converges in the following norm [30]: ° ± jku ur kj 6 C ku vk2, C k.u v/n k0,@h C ku vk0,@g 6 C kap1 jujpC1, ,. (13). (14). where C is a genetic constant and k is the overlapping number. This result indicates that for RKCM to converge, the RK approximation of degree p > 2 needs to be used. 3. GRADIENT REPRODUCING KERNEL APPROXIMATION Strong form collocation for second-order differential equations requires taking second-order differentiation on the RK shape functions of (6), which is time consuming, especially in calculating higher order derivatives of M1 .x/ at every evaluation point x. Motivated by the reproducing kernel approximation to achieve synchronized convergence in RKPM [35, 36], we consider the approximation of u,ˇ for the strong form of second-order PDEs as follows: u,ˇ wˇ D. Ns X. ‰Iˇ .x/aI ,. (15). I D1. where ‰Iˇ .x/ D C ˇ .xI x xI / 'a .x xI /. (16). P and ˇ D .ˇ1 , ˇ2 , : : : , ˇd /, jˇj D diD1 ˇi 6 j˛j. Here the correction functions in (16) are constructed with monomial bases of degree q C ˇ .xI x xI / D. q X. b˛ˇ .x/ .x xI /˛ , q > 0. j˛jD0 T. (17). ˇ. DW H .x xI / b .x/. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(5) GRADIENT REPRODUCING KERNEL COLLOCATION METHOD. 1385. The coefficients b˛ˇ are obtained from the following gradient reproducing conditions: Ns X. ‰Iˇ x˛I D D ˇ x˛ ,. 0 6 j˛j 6 q,. (18). I D1. where D ˇ @ˇ1 =@ˇ1 x1 @ˇ2 =@ˇ2 x2 : : : @ˇd =@ˇd xd . As shown in [6], Equation (18) is equivalent to Ns X. ‰Iˇ .x xI /˛ D .1/jˇ j D ˇ H .0/ ,. (19). I D1. where D ˇ H.x/ D. ˛Š x˛ˇ .˛ ˇ/Š. (20). and D ˇ H .0/ D ˛Šıˇ ˛ .. (21). Substituting (16) and (17) into (19) gives rise to M.x/bˇ .x/ D .1/jˇ j D ˇ H .0/ ,. (22). where M.x/ is the moment matrix given in (7). Consequently, the gradient RK shape functions are obtained as ‰Iˇ .x/ D .1/jˇ j D ˇ HT .0/ M1 .x/ H .x xI / 'a .x xI / .. (23). Uniqueness in Equation (22) requires that the kernel support to be large enough to ensure the nonsingularity of the moment matrix M.x/. This condition is identical to the requirement of kernel support in the reproducing kernel approximation [4, 5]. For this reason, truncation of kernel support is necessary in the discretization of nonconvex domain. It is noted that M.x/ is the Gram matrix of basis functions H .x xI / with respect to 'a .x xI /. The positivity of the kernel function 'a .x xI / ensures the positive definiteness of M.x/. In this work, 'a .x xI / is chosen to be the quintic B-spline kernel function: 8 2 4 5 11 ˆ 9s2 C 81s 81s , 0 6 s < 13 , ˆ 20 4 4 ˆ ˆ 3 4 5 < 17 15s 63s 2 kx xI k C 8 4 C 135s 243s C 81s , 13 6 s < 23 , 40 4 8 8 sD 'a .s/ D , (24) 2 3 4 5 ˆ 81 81s C 81s 81s C 81s 81s , 2 a ˆ 6 s < 1, ˆ 8 4 4 8 40 3 ˆ : 40 0, s > 1, where s is the normalized nodal distance. If equal order bases are used in the approximation of u and u,ˇ , the term M1 .x/ is identical in all shape functions ‰I and ‰Iˇ . Furthermore, by comparing the shape function for u in (6) with the shape functions for u,ˇ in (23), it appears that HT .0/ in (6) is replaced by .1/jˇ j D ˇ HT .0/ in (23), leading to a significant time saving in computing ‰Iˇ compared with a direct differentiation of ‰I . For the sake of simplicity but without loss of generality, we consider two-dimensional problems in this study. The approximation of u,x and u,y denoted as follows will be used in the next sections, and the simplified derivation of ‰Ix and ‰Iy is given in Appendix A. u,x wx D. Ns X. ‰Ix .x/aI. I D1. u,y wy D. Ns X. (25) ‰Iy .x/aI .. I D1. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(6) 1386. S.-W. CHI ET AL.. Hence, the second-order derivatives of u is obtained by taking direct derivatives of wx and wy , that is, u,xx wx,x D. Ns X. ‰Ix,x .x/aI. I D1. u,yy wy,y D. Ns X. (26) ‰Iy,y .x/aI .. I D1. 4. GRADIENT REPRODUCING KERNEL COLLOCATION METHOD To introduce gradient RK approximation in the discretization of strong form, consider the following boundary value problem: L1 u,x C L2 u,y D f B1h u,x. C B2h u,y. in . D h on @h. Bg u D g. (27). on @g ,. where L1 and L2 are the differential operators in , B1h and B2h are the boundary operators on @h , and Bg is the boundary operator on @g . The explicit forms of the operators and vectors for Poisson and elasticity problems in two dimensions are given in Table I. The approximations of u, u,x and u,y are given as u v D ‰ Ta T u,x wx D ‰ x a T u,y wy D ‰ y a,. (28). ® x ¯Ns ® y ¯Ns s where ‰, ‰ x , ‰ y , and a are the vector forms of ¹‰I ºN I D1 , ‰I I D1 , ‰I I D1 and aI , respectively. We define a least-squares functional associated with the boundary value problem in (27) with approximations u v, u,x wx , u,y wy as Z 1 1 T 1 E v, wx , wy D L wx C L2 wy f L wx C L2 wy f d 2 Z 1 T 1 ˛h C Bh wx C B2h wy h Bh wx C B2h wy h d (29) 2 @h Z T ˛g Bg u g Bg u g d. C 2 @g Table I. Explicit forms of operators for the Poisson and elasticity problems in two dimensions. Operator. Poisson’s problem. Elasticity problem ". L1. @ @x. L2. @ @y. " . B1h. nx. B2h. ny. Bg. 1. . Copyright © 2012 John Wiley & Sons, Ltd.. . C 2/. @ @x. @ @y @ @y @ @x. @ @x @ @x. . C 2/. . C 2/ nx ny ny nx 1 0. 0 1. @ @y. @ @y. ny nx. nx . C 2/ ny . # # . Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(7) 1387. GRADIENT REPRODUCING KERNEL COLLOCATION METHOD. Here the first term accounts for the least-squares residual of the differential equation in the domain, and the second and third terms account for the least-squares residuals of the Neumann and Dirichlet boundary conditions, respectively. Weights ˛h and ˛g are considered in the least-squares residuals for the boundary constraints. Substituting (28) into (29) and considering the stationary condition lead to the variational discrete equation Z. T T L1 ‰ x L1 ‰ x a C L2 ‰ y a f d Z T T T C ıa L2 ‰ y L1 ‰ x a C L2 ‰ y a f d Z T T T C ˛h ıa B1h ‰ x B1h ‰ x a C B2h ‰ y a h d @ Z h T T C ˛h ıaT B2h ‰ y B1h ‰ x a C B2h ‰ y a h d @ Z h C ˛g ıaT Bg ‰ Bg ‰ T a g d. T. ıE D ıa. (30). @g. Performing quadrature rules at the collocation points yields. T. ıE D ıa. Np X. T T T T L1 ‰ x .p` / L1 ‰ x .p` / a C L2 ‰ y .p` / a f .p` /. `D1 2. CL C ıaT ˛h. . ‰. yT. T 1 xT 2 yT .p` / L ‰ .p` / a C L ‰ .p` / a f .p` / w`1. Nq X. T T T T B1h ‰ x .q` / B1h ‰ x .q` / a C B2h ‰ y .q` / a h .q` /. `D1. T T T T CB2h ‰ y .q` / B1h ‰ x .q` / a C B2h ‰ y .q` / a h .q` / w`2. (31). Nr T X T T Bg ‰ .q` / C ıa ˛g Bg ‰ .q` / a g .q` / w`3 D 0, T. `D1. ® ® ¯Np ® ¯Nq ¯Nr where p` , w`1 `D1 , q` , w`2 `D1 , and r` , w`3 `D1 are the pairs of quadrature points and weights in and on @h and @g , respectively. We can rewrite (31) as h T T ıE D ıaT A1 W1 A2 a C A1 a b1 C A2 W1 A1 a C A2 a b1 T T C˛h A3 W2 A3 a C A4 a b2 C ˛h A4 W2 A3 a C A4 a b2 i T C˛g A5 W3 A5 a b3 . D ıaT AT W .Aa b/ D 0,. (32). where 0. 1 A1 C A2 3 p A D @ ˛h A C A4 A, p ˛g A5 Copyright © 2012 John Wiley & Sons, Ltd.. 0. 1 b1 p b D @ ˛h b2 A, p ˛g b3. 0 WD@. 1. W1 W2 3. A. (33). W. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(8) 1388. S.-W. CHI ET AL.. 0. T. L1 ‰ x .p1 / T L1 ‰ x .p2 / .. . 1 xT L ‰ pNp. B B A DB B @ 1. 0 B B A3 D B B @ 0 B B b1 D B @ 0 B W1 D @. T. B1h ‰ x .q1 / T B1h ‰ x .q2 / .. . 1 xT Bh ‰ qNq f .p1 / f .p2 / .. . f pNp. 1. 0. C C C, C A. B B A4 D B B @. 1 C C C, A. 0. w11 ... . 1 wN p. 0. C C C, C A. B B A DB B @. T. L2 ‰ y .p1 / T L2 ‰ y .p2 / .. . 2 yT L ‰ pNp. 2. T. B2h ‰ y .q1 / T B2h ‰ y .q2 / .. . 2 yT Bh ‰ qNq. 1. h .q1 / h .q2 / .. . h qNq 0. C A,. B W2 D @. B B b2 D B @. 1. 1 C C C, A. 1. 0. C C C, C A. B B A5 D B @. 0. g .r1 / g .r2 / .. .. B B b3 D B @. ... . 2 wN q. C C C C A. Bg ‰ T .r1 / Bg ‰ T .r2 / .. . Bg ‰ T .rNr /. 1 C C C A. (35). C C C, A 0 B W3 D @. C A,. (34). 1. g .rNr / 1. w12. 1. 1. w13 ... . 3 wN r. C A (36). From (32), the discrete weighted least-squares equation has the following form: AT W Aa D AT Wb.. (37). Equation (37) is the weighted least-squares approximation of the linear system Aa D b, that is 1 1 0 A1 C A2 b1 3 p p @ ˛h A C A4 A a D @ ˛h b2 A . p p ˛g b3 ˛g A5 ƒ‚ … „ ƒ‚ … „ 0. (38). b. A. The submatrices in matrix A, and the vectors a and b for Poisson and elasticity problems are summarized in Table II. Table II. Submatrices in discrete equations for Poisson and elasticity problems. Submatrix A1IJ A2IJ. Poisson’s problem h i ‰Jx ,x .pI / h. y ‰J ,y. .pI /. i. A3IJ. x. ‰J .qI / nx. A4IJ. y. ‰J .qI / ny. A5IJ. Elasticity problem " ". Copyright © 2012 John Wiley & Sons, Ltd.. y. ". ‰Jx ,y .pI / ‰Jx ,x .pI / y. ‰J ,y .pI /. ‰J ,x .pI /. ‰J ,x .pI /. . C 2/ ‰J ,y .pI /. y. y. ‰J .rI / 0. ‰Jx .qI / ny ‰Jx .qI / nx y. ‰J .qI / ny ‰J .qI / nx. # #. y. . C 2/ ‰Jx .qI / nx ‰Jx .qI / ny y. Œ‰J .rI /. . C 2/ ‰Jx ,x .pI / ‰Jx ,y .pI /. ‰J .qI / nx. #. y. . C 2/‰J .qI / ny . 0 ‰J .rI /. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(9) GRADIENT REPRODUCING KERNEL COLLOCATION METHOD. 1389. 5. CONVERGENCE STUDY We first consider a two-dimensional Poisson boundary value problem (BVP) as a model problem u,xx C u,yy Df in u D g on @g r u n un D h on @h .. (39). As discussed in Section 4, the strong form collocation can be related to the least-squares functional with quadrature. On the basis of the least-squares functional in (29) and considering the BVP in (39), E-norm is defined as follows: °. ±1=2. 2. v, wx , wy D wx,x C wy,y 2 C ˛h kwn k2 C ˛ kvk g 0,@h 0,@g 0, E. (40). where wn D wx nx C wy ny , and vD. Ns X. ‰I aI ,. v 2 V D span ¹‰1 , ‰2 , ‰NS º. ‰Ix aI ,. ¯ ® x wx 2 Wx D span ‰1x , ‰2x , ‰N S. ‰Iy aI. ° ± y wy 2 Wy D span ‰1y , ‰2y , ‰N S. I D1. wx D. Ns X I D1. wy D. Ns X. (41). I D1. Thus, we have °. u v, u,x wx , u,y wy D wx,x C wy,y f 2 C ˛h kwn hk2 0,@h E 0, ±1=2 C˛g kv gk20,@g. p 6 wx,x C wy,y f 0, C ˛h kwn hk0,@h p C ˛g kv gk0,@g. (42). E1 C E2 C E3 . Here E1 is the error from domain, E2 is the error from the Neumann boundary, and E3 is the error from Dirichlet boundary. The individual error norms are estimated as follows:. E1 D wx,x C wy,y f 0,. D wx,x C wy,y u,xx u,yy 0,. 6 kwx,x u,xx k0, C wy,y u,yy 0, (43). 6 kwx u,x k1, C wy u,y 1,. 6 C1 a1 kwx u,x k C C2 a1 wy u,y. 0,. 0,. p ˛h kwn hk0,@h p D ˛h kwn un k0,@h. p p. 6 CN 3 ˛h kwx u,x k1, C CN 4 ˛h wy u,y 1,. p p. 6 C3 a1 ˛h kwx u,x k0, C C4 a1 ˛h wy u,y 0,. E2 D. Copyright © 2012 John Wiley & Sons, Ltd.. (44). Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(10) 1390. S.-W. CHI ET AL.. p ˛g kv gk0,@g p D ˛g kv uk0,@g p 6 CN 5 ˛g kv uk1, p 6 C5 a1 ˛g kv uk0, .. E3 D. (45). We further introduce the following properties of the reproducing kernel approximation of degree p in (5) and the gradient reproducing kernel approximation of degree q in (18): ku vk0, 6 C apC1 jujpC1, ku,x wx k0, 6 C aq jujqC1,. u,y wy 6 C aq juj qC1, . 0,. (46). As can be seen, E1 and E2 are associated with the gradient reproducing kernel approximation (wx and wy ) of the differential equation and the Neumann boundary condition, respectively. It appears that the errors E1 and E2 are in balance without the weight in E2 , thus the weight ˛h is unnecessary. The error term E3 is associated with the reproducing kernel approximation (v), and its balance with the errors E1 and E2 requires the properties in (46). As such, the weights for imposition of boundary conditions in G-RKCM are selected as shown below p p ˛h O.1/, ˛g O aqp1 . (47) Combining the properties in (46) and the weights in (47), we have. . u v, u,x wx , u,y wy 6 aq1 C9 jujqC1, C C10 jujpC1, . E. (48). Assuming the discrete bilinear form associated with the minimization of E-norm in (42) is bounded and coercive, by Lax-Milgram and Cea’s lemmas, there exists an optimal estimate. u v, u,x wx , u,y wy. u uh , u,x uh,x , u,y uh,y 6 CN inf v2V E E wx 2Wx (49) wy 2Wy q1 6a C11 jujqC1, C C12 jujpC1, . Furthermore, considering the balance of errors in the E-norm, and the error properties in (43)–(45), we have. . (50). u uh O aq1 , u,x uh,x C u,y uh,y O aq1 1,. u uh. 0,. 1,. O .aq / ,. u,x uh,x. 0,. 1,. C u,y uh,y. 0,. O .aq /.. (51). For elasticity, similar procedures are followed to obtain. E1 6 C1 a1 kwx u,x k0, C C2 a1 wy u,y 0,. p p. E2 6 C3 a1 ˛h kwx u,x k0, C C4 a1 ˛h wy u,y 0, p E3 6 C5 a1 ˛g kv uk0, ,. (52). where D max ¹, º. For balance of errors between E1 , E2 , and E3 , the following weights are selected: p p (53) ˛g O aqp1 . ˛h O.1/, Similar convergence properties to the Poisson problem as given in (48)–(51) can be obtained for elasticity problems. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(11) 1391. GRADIENT REPRODUCING KERNEL COLLOCATION METHOD. Remark 5.1 The results in (49) indicate that the convergence of this method is only dependent on the polynomial degree q in the approximation of u,x and u,y , and is independent of the polynomial degree p in the approximation of u. Furthermore, q > 2 is needed for convergence. Remark 5.2 The collocation points in the strong form collocation method play a similar role as the quadrature points in the least-squares method as discussed in Section 3. For strong form collocation method based on approximation for u, such as the RKCM [30–32], it requires second-order differentiation of the approximation functions. Typically, higher order differentiation in the approximation function requires higher order quadrature rule for sufficient accuracy in the solution process. Taking RKCM for a Poisson problem for example, we have ˇ ˇ ˇZ ˇ Z^ ˇ ˇ ˇ vd vdˇ 6 C hrC1 N rC3 kvk2 . (54) c s 1, ˇ ˇ ˇ ˇ . . R^ In the above, denotes numerical integration, hc D 1=Nc , Nc and Ns are the numbers of collocation points and source points in one dimension, respectively, and r is the parameter related to the accuracy of numerical integration method, for example, r D 1 for Trapezoidal rule. Here, v D r rv involves second-order differentiation of the approximation in v. For the proposed G-RKCM, v is replaced by r Œwx , wy , which requires only first-order differentiation of wx and wy , and we have ˇ ˇ ˇZ ˇ Z^ ˇ ˇ ˇ r Œwx , wy d r Œwx , wy dˇ 6 C hrC1 N rC1 kwn k2 , (55) 1, c s ˇ ˇ ˇ ˇ . . where wn D wx nx C wy ny . NsrC3 D Nc.rC1/ NsrC3 o.1/ for integration error to be under For RKCM, it requires hrC1 c control, and thus necessitates the use of more collocation points Nc than source points Ns in the collocation method, and that leads to an over-determined system in its collocation equations. For the proposed G-RKCM, we need hrC1 NsrC1 D Nc.rC1/ NsrC1 o.1/, and thus allows the use of c Nc D Ns for sufficient accuracy as will be shown in the numerical examples. 6. COMPLEXITY ANALYSIS In this section, we analyze the complexity of RKCM and the proposed G-RKCM. For complexity comparison of RKCM and G-RKCM, consider the solution of the following Poisson problem: u D f uDg un D h. in on @g on @h ,. (56). where D rr and un D r un. We consider the following two formulations in the approximations: RKCM W u v D. Ns X. ‰I .x/aI , u,˛ v,˛ D. I D1. Ns X. ‰I ,˛ .x/aI , u,˛˛ v,˛˛ D. I D1. Ns X. ‰I ,˛˛ .x/aI. I D1. (57) GRKCM W u v D. Ns X. ‰I .x/aI , u,˛ w˛ D. I D1. Ns X I D1. ‰I˛ .x/aI , u,˛˛ w˛,˛ D. Ns X. ‰I˛,˛ .x/aI ,. I D1. (58) Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(12) 1392. S.-W. CHI ET AL.. where ˛ D 1, 2, ‰I .x/ is the RK shape function of degree p, and ‰I˛ is the gradient RK shape function with degree q. Consider a set of collocation points ° ± Np Nq r C D ¹p` º`D1 , ¹q` º`D1 , ¹r` ºN (59) , p` 2 , q` 2 @g , r` 2 @h . `D1 Introducing RK approximation in (57) into the strong form (56), and enforcing the residual to be zero at the collocation points to yield . Ns X. ‰I .p` / aI D f .p` / ,. p` 2 ,. ` D 1, , Np. (60a). I D1 Ns p p X ˛g ‰I .q` / aI D ˛g g .q` / ,. q` 2 @g ,. ` D 1, , Nq. (60b). I D1 N. s p X p .r ‰I .r` / n .r` // aI D ˛h h .r` / , ˛h. r` 2 @h ,. ` D 1, , Nr. (60c). I D1. Note that for RKCM, the second-order derivative on the shape function ‰I D ‰I ,xx C ‰I ,yy is needed in (60a), while for G-RKCM, this term is replaced by ‰I D ‰Ix,x C ‰Iy,y . Similarly. in (60c), r‰I D Œ‰I ,x , ‰I ,y for RKCM, while r‰I D ‰Ix , ‰Iy for G-RKCM. It is therefore imperative to analyze the operating counts of ‰I , ‰I ,˛ , and ‰I˛ , ˛ D 1, 2 in the following. We denote multiplication and division operations by M/D, and the addition and subtraction operations by A/S. For RKCM, operation counts following [30] are: ² M/D: S 3 C .2k C 1/ S 2 C S C 1 ‰I (61) A/S W S 3 C .k 2/ S 2 C S 1 ² ‰I ,˛ ² ‰I ,˛˛. M/D: 3S 3 C .8k C 4/ S 2 C 3S C 2 A/S: 3S 3 C .4k 5/ S 2 C S C 1. M/D W A/S W. 6S 3 C .20k C 12/ S 2 C 6S C 4 , 6S 3 C .10k 11/ S 2 C S C 12. (62). (63). where S D .p C d /Š=.pŠd Š/, p is the reproducing degree of RK approximation, d is the space dimension, and k is the kernel support overlapping number. For G-RKCM, the operating count for ‰I is the same as (61), and the operating counts for ‰I˛ and ‰I˛,˛ are ² M/D: SN 3 C .2k C 1/ SN 2 C SN C 1 ‰I˛ (64) A/S: SN 3 C .k 2/ SN 2 C SN 1 ‰I˛,˛. ². M/D W A/S:. 3SN 3 C .8k C 4/ SN 2 C 3SN C 2 , 3SN 3 C .4k 5/ SN 2 C SN C 1. (65). where SN D .q C d /Š=.qŠd Š/, and q is the reproducing degree of gradient RK approximation. Note that the complexity of ‰I˛ is the same as that for ‰I , and the complexity of ‰I˛,˛ is the same as that for ‰I ,˛ , with p in S replaced by q in SN . The computational complexities of these shape functions in two dimensions for p D q D 2 are shown in Tables III and IV, respectively. The kernel support overlapping number is based on normalized kernel support of k D 4S in two dimensions. The collocation equations in (60) lead to a linear system Aa D b. Copyright © 2012 John Wiley & Sons, Ltd.. (66) Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(13) 1393. GRADIENT REPRODUCING KERNEL COLLOCATION METHOD. Table III. Complexity comparison of shape function calculation in RKCM and G-RKCM in two dimensions. pD2. RKCM M/D A/S M/D A/S M/D A/S. ‰I ‰I ,˛ ‰I ,˛˛. 1987 1013 7724 3931 19,048 9558. pDqD2. G-RKCM M/D A/S M/D A/S M/D A/S. ‰I ‰I˛ ‰I˛,˛. 1987 1013 1987 1013 7724 3991. Table IV. Complexity comparison of shape function calculation in RKCM and G-RKCM in three dimensions. pD2. RKCM ‰I ‰I ,˛ ‰I ,˛˛. M/D A/S M/D A/S M/D A/S. 17,111 8809 67,642 34,511 167,684 84,922. pDqD2. G-RKCM ‰I ‰I˛ ‰I˛,˛. M/D A/S M/D A/S M/D A/S. 17,111 8809 17,111 8809 67,642 34,511. In (66), the matrix A is with dimension Nc Ns , where Nc D Np C Nq C Nr is the total number of collocation points, Ns is the number of source points, Nc > Ns for RKCM, while Nc D Ns for G-RKCM. Thus, the solution time for solving the linear system (66) also favors G-RKCM in addition to its simplicity in shape function calculations as discussed above. Furthermore, the computation time in constructing the linear system in G-RKCM is also considerably less than that in RKCM. For example, let NN p , NN q , and NN r be the counter parts of Np , Nq , and Nr in G-RKCM, and NN c D NN p C NN q C NN r D Ns . It can be shown that the construction times for the linear system of (66) are RKCM W Ns 38, 096Np C 1987Nq C 7724Nr. (67). G-RKCM W Ns 15, 548NN p C 1987NN q C 1987NN r .. (68). By considering the fact that Np C Nq C Nr > NN p C NN q C NN r as discussed above, the CPU advantage in G-RKCM is trivial. 7. NUMERICAL EXAMPLES In the following numerical examples, both RK shape functions and gradient RK shape functions are constructed with the quintic B-spline kernel function. For comparison, the solutions of the proposed G-RKCM method are compared with analytical solutions and RKCM solutions. In the solution of p p BVP, the boundary weights of ˛h D 1, ˛g D Ns are used for RKCM following [15], while p p ˛h D 1, ˛g D aqp1 are used for G-RKCM, with D 1 for the Poisson problem and D max ¹, º for elasticity as discussed in Section 5. In the convergence plots of the numerical examples, the numbers shown in the legends represent the rate of convergence of a given norm. 7.1. Approximation of a sine function The RK shape functions and gradient RK shape functions are employed to approximate sin. x/ sin. y/ and the associated derivative in the domain Œ0, 1 Œ0, 1, respectively. The L2 error norms of the function approximation by RK shape function with p D 1 and p D 2 are shown Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(14) 1394. S.-W. CHI ET AL.. Figure 1. Convergence of L2 norms in approximating a sine function and its derivative.. in Figure 1(a) while the L2 error norms of the approximation of sine function derivative with q D 1 and q D 2 are shown in Figure 1(b). The same number of collocation points as that of source points is used in this study. The convergence rates are in agreement with the theoretical values. 7.2. Two-dimensional Poisson problem Consider a two-dimensional Poisson problem as follows: u .x, y/ D x 2 C y 2 e xy in D .0, 1/ .0, 1/ u .x, y/ D e xy on @.. (69). The numbers of source points and collocation points employed for RKCM in the convergence study are ¹10 10, 15 15, 20 20, 25 25, 30 30º and ¹19 19, 29 29, 39 39, 49 49, 59 59º, respectively, and the number of collocation points are the same as the number of source points ¹10 10, 15 15, 20 20, 25 25, 30 30º for G-RKCM. Figure 2 compares L2 norms of u and u,˛ obtained by the proposed G-RKCM with various degrees of bases, and RKCM with p D 2. As predicted by the theory in Section 4, G-RKCM requires at least second-order bases in the gradient RK approximation for convergence, similar to the convergence requirement for RKCM [30]. The results also show that the rate of convergence in G-RKCM is determined by the degree of bases in the gradient RK approximation (q), although higher degree of bases in the RK approximation (p) improves the solution accuracy in u. It is also shown that the L2 error norms of u and u,˛ have essentially the same convergence rates and are consistent with the error analysis results given in Section 4. The CPU comparison for RKCM and G-RKCM shown in Figure 3 demonstrates the effectiveness of the proposed G-RKCM.. Figure 2. Convergence of L2 norms of u and u,˛ in two-dimensional Poisson problem. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(15) GRADIENT REPRODUCING KERNEL COLLOCATION METHOD. 1395. Figure 3. CPU comparison of RKCM and G-RKCM.. 7.3. Infinite long cylinder under internal pressure An infinite long elastic cylinder subjected to an internal pressure is depicted in Figure 4 , where a plane strain condition in the out of plane direction is assumed. Because of symmetry, only a quarter of the domain is modeled by G-RKCM as shown in Figure 5(a). The corresponding boundary value problem is ij ,j D 0 in . (70). on 1 hi D P ni u2 D 0, h1 D 0 on 2 hi D 0 on 3 u1 D 0, h2 D 0 on 4 ,. (71). Figure 4. An infinite long cylinder subjected to an internal pressure.. Figure 5. (a) Quarter model and (b) distribution of source points and collocation points for RKCM. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(16) 1396. S.-W. CHI ET AL.. where ij D Cij kl u.k,l/ and hi D ij nj . The analytical solutions to this problem are given by P a2 r b2 .1 C

(17) / .1 2

(18) / C 2 .1 C

(19) / ur .r/ D E .b 2 a2 / r 2 2

(20) b Pa rr .r/ D 2 1 2 b a2 r

(21) 2 b2 Pa 1C 2 , .r/ D 2 b a2 r. (72). where P is the internal pressure, b and a are the outer and inner radii of the cylinder, respectively. The distribution of source points and collocation points for RKCM is shown in Figure 5(b). Five levels of discretization with source points ¹66, 222, 469, 808, 1238º are employed in the convergence study. The number of collocation points is approximately four times the source points for the RKCM whereas the collocation points are the same as the source points for G-RKCM. As shown in Figure 6, disregarding the degree of basis p, the G-RKCM with quadratic basis q D 2 achieves the similar rate of convergence as the RKCM with quadratic basis while it yields better accuracy than RKCM in this problem. The errors in the G-RKCM with q D 2 along the radial direction are also compared with those in the RKCM in Figure 7. In general, the stress results obtained by G-RKCM are less oscillatory in comparison with those by RKCM. The results also show that for G-RKCM, the solution is predominated by the order of basis functions (q) in the gradient RK shape functions, and is nearly independent of the order of basis functions (p) in the RK shape functions. The condition number of matrix A in G-RKCM (Equation (38)) is compared with those of corresponding matrices in RKCM and RBCM in Figure 8, where in RBCM, the commonly used RBF of the following form is adopted with a shape parameter c D 3 in this problem: 12 gI .x/ D .x xI /2 C c 2 .. (73). As can be seen in Figure 8, the condition number of the linear system in G-RKCM is the smallest among others and it grows with the lowest rate when the number of points increases. 7.4. Beam under shear load Consider a plane-strain elastic cantilever beam subjected to a tip shear traction P shown in Figure 9. The corresponding boundary value problem and boundary conditions are given as ij ,j D 0, 0 < x < L, D=2 < y < D=2. (74). Figure 6. Convergence of L2 norms of u and u,˛ in the cylinder problem. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(22) GRADIENT REPRODUCING KERNEL COLLOCATION METHOD. 1397. Figure 7. Displacement and stresses along radial direction of the cylinder.. Figure 8. Condition numbers of discrete equations in G-RKCM, RKCM, and RBCM.. Figure 9. Cantilever problem statement. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(23) 1398. S.-W. CHI ET AL.. .1/ at x D 0, y D 0, u1 D u2 D 0 .2/ at x D 0, y D ˙D=2, u1 D 0, h2 D 0 .3/ on x D L, D=2 6 y 6 D=2, h1 D 0, h2 D. 6P D3. .4/ on x D 0, Dt =2 < y < 0, 0 < y < D=2, h1 D .5/ on 0 < x < L, y D ˙D=2, h1 D h2 D 0.. . . D2 y2 4 12PL y, h2 D3. D. 6P D 3. . D2 4. y. The analytical solutions to the problem are h i 2 N y 2 D4 u1 .x, y/ D 6PENyI .6L 3x/ x C .2 C

(24) / h i , P D2 x 2 2 .3L .L .4 u2 .x, y/ D 6E C 3

(25) y N x/ x x/ C C 5

(26) / N NI 4. 2. . (75). (76). where I D D 3 t =12 , EN D E= 1

(27) 2 , and

(28) N D

(29) = .1

(30) /. Six levels of discretization are performed in the convergence study with source points ¹175, 25 7, 33 9, 41 11, 49 13, 57 15º, and collocation points in both G-RKCM and RKCM are the same as the source points in this problem. The L2 norms of u and u,˛ obtained by the proposed G-RKCM with various degrees of bases are compared with those obtained by RKCM (p D 2) in Figure 10. Again, almost independent of the degree of basis p, the G-RKCM with quadratic basis q D 2 achieves the similar rate of convergence as the RKCM with quadratic basis. The comparison of shear stress solutions along x D L=2 obtained by G-RKCM with q D 2 and RKCM with p D 2 is shown in Figure 11, where Ns D 25 7 is used. The results of shear stress obtained by G-RKCM are less oscillatory compared with that obtained by RKCM.. Figure 10. Convergence of L2 norms of u and u,˛ in cantilever problem.. Figure 11. Comparison of shear stress along x D L=2 in cantilever problem. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(31) GRADIENT REPRODUCING KERNEL COLLOCATION METHOD. 1399. 8. CONCLUSION Although RKCM using direct RK approximation of strong form has shown an enhanced conditioning and sparsity in its discrete system compared with RBCM using radial basis approximation of strong form, and it also resolved the domain integration issues in the weak form based Galerkin meshfree method, the method suffers from the high level of complexity involved in computing the second-order derivatives of RK shape functions and the need of using the number of collocation points much larger than the number of source points for optimal convergence. To resolve these issues, in this work we propose a G-RKCM by formulating the derivatives of RK shape functions directly based on the partition of nullity and discrete derivative reproducing conditions to eliminate the need of taking second derivatives of the Gram matrix involved in RKCM for solving second-order PDEs. We also showed that in the proposed G-RKCM the number of collocation points needs not be greater than that of source points required in RKCM. The error analysis showed that the rate of convergence in G-RKCM is determined by the polynomial degree in the gradient RK approximation, and is independent of the polynomial degree in the RK approximation. Furthermore, G-RKCM yields the same convergence rates in L2 norms of u and u,˛ . The complexity analysis provided precious operating counts of both RKCM and G-RKCM and clearly demonstrated the significant computational efficiency of G-RKCM over RKCM. The numerical results confirmed the analytical predictions, and showed that the proposed G-RKCM yields similar convergence property as the RKCM in both L2 norms of u and u,˛ , yet it is roughly 10 times computationally more efficient than RKCM. APPENDIX A Consider the approximation of u,x and u,y in two dimensions as follows: u,x wx D. Ns X. ‰Ix .x/aI. I D1. u,y wy D. Ns X. (A.1) ‰Iy .x/aI ,. I D1. where ‰Ix .x/ D C 1 .xI x xI / 'a .x xI / ‰Iy .x/ D C 2 .xI x xI / 'a .x xI /. .. (A.2). For demonstration purpose, consider a case with linear bases q D 1 in two dimensions: i i i C i .xI x xI / D b00 .x/ C b10 .x/ .x xI / C b01 .x/ .y yI / DW HT .x xI / bi .x/,. i D 1, 2, (A.3). where the coefficients b˛i 1 ˛2 .x/ are determined by satisfying the partition of nullity and first-order derivative reproducing conditions shown below: Ns X. ‰Ix .x/ D 0,. I D1. Ns X. Ns X. ‰Ix .x/xI D 1,. I D1. ‰Iy .x/ D 0,. I D1. Copyright © 2012 John Wiley & Sons, Ltd.. Ns X I D1. Ns X. ‰Ix .x/yI D 0. (A.4). ‰Iy .x/yI D 1.. (A.5). I D1. ‰Iy .x/xI D 0,. Ns X I D1. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(32) 1400. S.-W. CHI ET AL.. From (A.4), multiplying (A.4a) by x and subtracting (A.4b) leads to Ns X. ‰Ix .x/ .x xI / D 1.. (A.6). I D1. Similarly, multiplying (A.4a) by y and subtracting (A.4c) yields Ns X. ‰Ix .x/ .y yI / D 0.. (A.7). I D1. Applying the same procedures to (A.5), we have Ns X. ‰Iy .x/ .x xI / D 0. (A.8). ‰Iy .x/ .y yI / D 1.. (A.9). I D1. Ns X I D1. The first-order derivative reproducing conditions in (A.6) (A.9) can be equivalently written as Ns X. ‰Ix .x/ D 0,. I D1. Ns X I D1. Ns X. ‰Ix .x/ .x xI / D 1,. I D1. ‰Iy .x/ D 0,. Ns X. Ns X. ‰Ix .x/ .y yI / D 0. (A.10). ‰Iy .x/ .y yI / D 1.. (A.11). I D1. ‰Iy .x/ .x xI / D 0,. I D1. Ns X I D1. From which we can express the first-order derivative reproducing conditions (A.10) and (A.11) as Ns X. ‰Ix .x/H .x xI / D H,x .0/. (A.12). ‰Iy .x/H .x xI / D H,y .0/.. (A.13). I D1. Ns X I D1. Substituting (A.2) and (A.3) into (A.12) and (A.13) gives rise to M.x/b1 .x/ D H,x .0/. (A.14). M.x/b2 .x/ D H,y .0/,. (A.15). where M.x/ is the moment matrix given in (7). Consequently, the gradient RK shape functions are obtained as 1 ‰Ix .x/ D HT ,x .0/M .x/H .x xI / 'a .x xI / 1 ‰Iy .x/ D HT ,y .0/M .x/H .x xI / 'a .x xI /. Copyright © 2012 John Wiley & Sons, Ltd.. (A.16). Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(33) GRADIENT REPRODUCING KERNEL COLLOCATION METHOD. 1401. REFERENCES 1. Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent development. Computer Methods in Applied Mechanics and Engineering 1996; 139:3–49. 2. Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. International Journal for Numerical Methods in Engineering 1994; 37:229–256. 3. Lancuster P, Salkauskas K. Surface generated by moving least squares methods. Mathematics of Computation 1981; 37:141–158. 4. Chen JS, Pan C, Wu CT, Liu WK. Reproducing kernel particle methods for large deformation analysis of nonlinear structures. Computer Methods in Applied Mechanics and Engineering 1996; 139:195–227. 5. Liu WK, Jun S, Zhang YF. Reproducing kernel particle methods. International Journal for Numerical Methods in Fluids 1995; 20:1081–1106. 6. Liu WK, Li S, Belytschko T. Moving least-square reproducing kernel methods (i) Methodology and convergence. Computer Methods in Applied Mechanics and Engineering 1997; 143:113–154. 7. Babuska I, Melenk JM. The partition of unity method. International Journal for Numerical Methods in Engineering 1997; 40:727–758. 8. Han W, Meng X. Error analysis of the reproducing kernel particle method. Computer Methods in Applied Mechanics and Engineering 2001; 190:6157–6181. 9. Babuska I, Banerjee U, Osborn JE, Zhang Q. Effect of numerical integration on meshless methods. Computer Methods in Applied Mechanics and Engineering 2009; 198:2886–2897. 10. Beissel S, Belytschko T. Nodal integration of the element-free Galerkin method. Computer Methods in Applied Mechanics and Engineering 1996; 139:49–74. 11. Bonet J, Kulasegaram S. Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations. International Journal for Numerical Methods in Engineering 2000; 47:1189–1214. 12. Chen JS, Wu CT, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin mesh-free methods. International Journal for Numerical Methods in Engineering 2001; 50:435–466. 13. Puso MA, Chen JS, Zywicz E, Elmer W. Meshfree and finite element nodal integration methods. International Journal for Numerical Methods in Engineering 2008; 74:416–446. 14. Aluru NR. A point collocation method based on reproducing kernel approximation. International Journal for Numerical Methods in Engineering 2000; 47:1083–1121. 15. Hu HY, Chen JS, Hu W. Weighted radial basis collocation method for boundary value problems. International Journal for Numerical Methods in Engineering 2007; 69:2736–2757. 16. Kansa EJ. Multiqudrics—a scattered data approximation scheme with applications to computational fluiddynamics. I. Surface approximations and partial derivatives. Computers & Mathematics with Applications 1992; 19:127–145. 17. Kansa EJ. Multiqudrics—a scattered data approximation scheme with applications to computational fluiddynamics. II. Solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers & Mathematics with Applications 1992; 19:147–161. 18. Kim DW, Kim Y. Point collocation methods using the fast moving least-square reproducing kernel approximation. International Journal for Numerical Methods in Engineering 2003; 56:1445–1464. 19. Onate E, Idelsohn S, Zienkiewicz OC, Taylor RL. A finite point method in computational mechanics. Application to convective transport and fluid flow. International Journal for Numerical Methods in Engineering 1996; 39:3839–3866. 20. Zhang X, Liu XH, Song KZ, Lu MW. Least-squares collocation meshless method. International Journal for Numerical Methods in Engineering 2001; 51:1089–1100. 21. Franke R. Scattered data interpolation: tests of some methods. Mathematics of Computation 1982; 98:181–200. 22. Hardy RL. Multiquadric equations of topography and other irregular surfaces. Journal of Geophysical Research (1971); 176:1905–1915. 23. Hardy RL. Theory and applications of the multiquadric-biharmonic method: 20 years of discovery. Computers & Mathematics with Applications 1990; 19:163–208. 24. Madych WR. Miscellaneous error bounds for multiquadric and related interpolatory. Computers and Mathematics with Applications 1992; 24:121–138. 25. Franke R, Schaback R. Solving partial differential equations by collocation using radial functions. Applied Mathematics and Computation 1998; 93:73–82. 26. Madych WR, Nelson SA. Bounds on multivariate polynomials and exponential error estimates for multiquadric interpolation. Journal of Approximation Theory 1992; 70:94–114. 27. Wendland H. Meshless Galerkin methods using radial basis functions. Mathematics of Computation 1999; 68:1521–1531. 28. Hon YC, Schaback R. On unsymmetric collocation by radial basis functions. Applied Mathematics and Computation 2001; 119:177–186. 29. Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations. Computers & Mathematics with Applications 2000; 4:123–137. 30. Hu HY, Chen JS, Hu W. Error analysis of collocation method based on reproducing kernel approximation. Numerical Methods for Partial Differential Equations 2011; 27 (3):554–580. 31. Hu HY, Lai CK. A study on convergence and complexity of reproducing kernel collocation method. Interaction Multiscale Mech 2009; 2:295–319. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(34) 1402. S.-W. CHI ET AL.. 32. Chen JS, Hu W, Hu HY. Reproducing kernel enhanced local radial basis collocation method. International Journal for Numerical Methods in Engineering 2008; 75:600–627. 33. Chen JS, Wang L, Hu HY, Chi SW. Subdomain radial basis collocation method for heterogeneous media. International Journal for Numerical Methods in Engineering 2009; 80:163–190. 34. Wang L, Chen JS, Hu HY. Subdomain radial basis collocation method for fracture mechanics. International Journal for Numerical Methods in Engineering 2010; 83(7):851–876. 35. Li S, Liu WK. Reproducing kernel hierarchical partition of unity Part I: Formulation and theory. International Journal for Numerical Methods in Engineering 1999; 45:251–288. 36. Li S, Liu WK. Reproducing kernel hierarchical partition of unity Part II: Applications. International Journal for Numerical Methods in Engineering 1999; 45:289–317. 37. Chen JS, Zhang X, Belytschko T. An implicit gradient model by a reproducing kernel strain regularization in strain localization problems. Computer Methods in Applied Mechanics and Engineering 2004; 193:2827–2844. 38. Li S, Liu WK. Synchronized reproducing kernel interpolant via multiple wavelet expansion. Computational Mechanics 1998; 21:28–47.. Copyright © 2012 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2013; 93:1381–1402 DOI: 10.1002/nme.

(35)