A Coupled RBF Method for the Solution of Elastostatic Problems

Radial basis function (RBF) has been widely used in many scientific computing and engineering applications, for instance, multidimensional scattered data interpolation and solving partial differential equations. However, the accuracy and stability of the RBF methods often strongly depend on the shape parameter. A coupled RBF (CRBF) method was proposed recently and successfully applied to solve the Poisson equation and the heat transfer equation (Appl. Math. Lett., 2019, 97: 93–98). Numerical results show that the CRBF method completely overcomes the troublesome issue of the optimal shape parameter that is a formidable obstacle to global schemes. In this paper, we further extend the CRBF method to solve the elastostatic problems. Discretization schemes are present in detail. With two elastostatic numerical examples, it is found that both numerical solutions of the CRBF method and the condition numbers of the discretized matrices are almost independent of the shape parameter. In addition, even if the traditional RBF methods take the optimal shape parameter, the CRBF method achieves better accuracy.


Introduction
A radial function ϕ(r j ) is a function of the Euclidean norm r j � ‖x − x j ‖ 2 , where x ∈ R n is the center point and x j ∈ R n is a point in the influence domain of x. In the past few decades, the radial function has been used as a special basis function for solving interpolation problems and discretizing partial differential equations (PDEs) by means of collocation techniques [1][2][3], which is referred to as the radial basis function (RBF) method in this paper. e interest in the RBF method has three principal reasons [4]: (i) the approximate value by using RBF can be estimated without using meshes; (ii) RBF gives very accurate results both for interpolation problems and for solving partial differential equations; (iii) there is enough flexibility in the choice of basis functions.
Many RBFs have been proposed and carefully studied in the past few decades [1]. Most of them can be categorized into two groups: infinitely smooth and piecewise smooth (see some typical RBFs listed in Table 1). For the former group, the infinitely smooth RBFs (such as Multiquadric (MQ), Inverse Multiquadric (IMQ), and Gaussian (GA)) lead to spectral convergence, which is a great advantage in actual applications. However, these types of RBFs contain a user-defined positive shape parameter, called c, which controls the stability and accuracy of the RBFs approximation. e optimal shape parameter c may result in an illconditioned linear system [5,6]. us, the value of the shape parameter has to be carefully selected. In contrast, the RBFs associated with the latter group are not infinitely differentiable and seem to be easier to implement than those in the former group due to the free of shape parameter. However, they only lead to an algebraic rate of convergence and thus are rarely used alone in applications. erefore, many researchers prefer to use the former group in many scientific computing and engineering applications.
As stated above, the choice of a suitable value of the shape parameter in the infinitely smooth RBFs is very important and has been an ongoing challenge problem. Up to now, no mathematical theory has been developed to rule the selection of an adequate value but only several approaches are available in the open literature [7][8][9][10][11]. In the early works, many researchers tried to present some empirical formulas based on the number and distribution of data points to select a good value for the shape parameter; see [12][13][14] for examples. ese empirical formulas work well for some special problems but may not be good for other problems. By testing numerous numerical experiments, Rippa [15] pointed out that the value of the optimal shape parameter depends on not only the number and distribution of data points, but also the data vector and the precision of the computation. In the same paper, Rippa also proposed a very interesting algorithm, called the Leave-One-Out Cross Validation (LOOCV) algorithm, to select a good value of the shape parameter. Recently, some improvements of the LOOCV algorithm have been deeply studied; see [16][17][18][19][20] and references therein. Other techniques for finding the optimal shape parameter include the energy-based method [7,9], the sample solution approach [10], the interval method [21], and so on. Another type of technique [22,23], i.e., bypassing the problem of the optimal shape parameter selection by reducing its influence on the stability of the method, is also of great interest. One of the representative works for the latter approach is the coupled RBF (CRBF) method studied recently in [23]. e CRBF method is developed by coupling the classical infinitely smooth RBF with the conical spline. Two numerical examples, i.e., the Poisson equation and the heat transfer equation, have been studied in [23] and showed that the CRBF method sufficiently inherits the high accuracy of the infinitely smooth RBF and the stability property of the piecewise smooth RBF. Most importantly, the CRBF method produces a relatively well-conditioned solution system and maintains almost invariable accuracy as the shape parameter ranges from 1 to 10 9 . erefore, the CRBF method completely overcomes the troublesome issue of the optimal shape parameter that is a formidable obstacle to global schemes for solving these two types of PDEs.
In this paper, we extend the CRBF method, which belongs to a class of meshfree methods, to solve the elastostatic problems to further show its high accuracy and stability. In fact, the application of the meshfree method for the solution of the elastostatic problems is not a new issue. Many meshfree methods have been studied in the past few decades.
ere are the global RBF and local RBF based on strong formulation (also called the meshfree collocation method or Kansa's method) [24][25][26][27][28][29], the global Galerkin weak form (called the radial point interpolation meshfree method) [30][31][32], the local Petrov-Galerkin weak from (called the local radial point interpolation meshfree method) [33], the strong-weak form [34], and so on. For good surveys of meshfree methods, we refer the readers to [35][36][37][38]. Compared with other meshfree methods, the classical RBF method based on strong formulation has the advantages of simple form, easy implementation, and high precision. e objective of this work is to apply the CRBF method studied recently in [23] to solve the elastostatic problems. Discretization scheme of the CRBF method for solving a typical two-dimensional plane stress problem is present in detail. With two numerical examples, our findings show that the CRBF method has better accuracy and convergence rate than the traditional RBF method even if the latter takes the optimal shape parameter. Furthermore, both the accuracy of the CRBF method and the condition numbers of the discretized matrices are also almost independent of the shape parameter.
e rest of this work is organized as follows. In Section 2, we first briefly review the CRBF method and describe the formulation of the CRBF method for the solution of the elastostatic problems. In Section 3, we apply the CRBF method to specifically solve two well-known problems, i.e., the Timoshenko's beam and an infinite plate with a central hole. Both regular and irregular nodes distribution are considered. Condition numbers of the discretized system matrices and accuracy of the CRBF method are studied in detail with a large range of the shape parameter. Comparison of the accuracy of the CRBF method and the classical RBF methods is present. Finally, we end this work with some conclusions and outlooks in Section 4.

The CRBF Method for Elastostatic Problems
In this section, we first briefly review the coupled radial basis function (CRBF) as well as using the CRBF method for interpolation and solving general boundary value problem of partial differential equations. en, we apply the CRBF method to specifically solve the elastostatic problems.

e Construction of CRBFs.
By redefining the expression c 2 + r 2 j in the classical MQ, IMQ, IQ RBFs as (c 2 + r 2 j )/c 2 and coupling the redefined RBFs with the conical spline r 5 j , Zhang [23] proposed a new class of RBFs, which are referred to as the CRBFs in this paper; see specific expressions in Table 2. Note that there still exists the shape parameter, c, tuned by users in the CRBFs. Although there are no theoretical results to guarantee the accuracy and stability of the CRBF method, numerous numerical experiments from the Poisson equation and the heat transfer equation studied in [23] show that the CRBF method has much better accuracy than the classical RBF method and, most importantly, completely overcomes the troublesome issue of the optimal shape parameter. Now, we first turn to use the CRBFs defined in Table 2 to construct the interpolant s(x) for approximating the function f(x) by a linear combination of translates of single CRBF on a scattered data set (x j , f j ) N j�1 . Here, x j where ϕ CRBF (‖x − x j ‖ 2 ) is a coupled radial basis function defined as in Table 2 and λ j N j�1 are the coefficients. For simplicity, the superscript "CRBF" will be omitted hereafter. e unknown coefficients λ j N j�1 can be determined by enforcing the following interpolation conditions: which leads to a system of linear equations where A is the N × N interpolation matrix with entries (A) i,j � ϕ(‖x i − x j ‖ 2 ) and f � (f 1 , . . . , f N ) T is the righthand side vector including the given ordinates at centers. If the interpolation matrix A is nonsingular, by solving the linear system (3), the coefficients λ j are determined and thus, the interpolant function s(x) will be obtained. e coupled radial basis function approximation can be also used as a reliable tool for solving PDEs. Consider the following boundary value problem: where L is a linear differential operator, B is a boundary differential operator, Ω is a bounded and connected domain, and Γ is the boundary of Ω. Assume a set of distinct nodes with N I interior nodes being located inside the domain Ω and N b boundary nodes along its boundary Γ. We are looking for an approximation to the solution s(x) of equation (4) in the vector space spanned by the CRBFs ϕ(‖x − x j ‖ 2 ), j � 1, . . . , N. Suppose that (1) is an approximate solution. For the interior and boundary nodes, we have where By solving the system of linear equations (6), the vector of unknown coefficients λ will be determined and an approximation solution of (4) is obtained.

e Elastostatic Problem and the CRBF Formulation.
For simplicity, we only consider the two-dimensional elastostatic problem here. Note that the CRBF method can be easily extended to solve three-dimensional problems. For the plane stress problem, the elastostatic equation, written in terms of displacements, is where u and v are the horizontal and vertical displacements, σ xx and σ yy are the normal stresses, τ xy � τ yx is the shear stress, f x and f y are the given body forces, u and v are the given displacement constraints on the displacement Mathematical Problems in Engineering boundary Γ u , t x and t y are the given stress constraints on the stress boundary Γ t , E and μ stand for Young's modulus and Poisson's ratio, and l and m represent the cosine of the normal direction outside the slope. e stresses σ xx , σ yy , τ xy and the displacements u, v have the following relationship: By using (10), the stress boundary condition (9) can be expressed in terms of displacements. On the contrary, once the elastostatic equations (7)-(9) are solved for the displacements, the corresponding stresses can be obtained through (10).
Next, we apply the CRBF method studied in the above subsection to solve elastostatic equations (7)- (9). Assume a set of distinct nodes x i with N I interior nodes being located inside the domain Ω, N b 1 nodes along the displacement boundary Γ u , and N b 2 nodes along the stress boundary Γ t . In the CRBF method, the displacement approximate solutions u(x) and v(x) are separately spanned by a set of translated CRBFs, i.e., where a j and b j (j � 1, . . . , N) are the unknown coefficients and ϕ(r j ) with r j � ‖x − x j ‖ 2 is the coupled radial basis function defined as in Table 2. To get approximate solutions u(x) and v(x), the unknown coefficients a j and b j are computed by collocation of (7) at a given set of interior nodes and collocation of the boundary conditions (8)-(9) at boundary nodes. Let λ be a vector of length 2N containing the unknowns Substituting (11) into the elastostatic equation (7) and boundary conditions (8)-(9) leads to the following system of linear equations: where the submatrix A L ∈ R 2N I ×2N and the subvector f L ∈ R 2N I correspond to the interior nodes (i � 1, . . . , N I , j � 1, . . . , N) with If the linear system (13) is solvable, i.e., the coefficient matrix A is nonsingular, the coefficients a j and b j can be obtained. As with the traditional RBF method, it is very difficult to prove the invertibility of CRBF coefficient matrix (13). However, numerous numerical results (see [23] and also the next section) show that the coefficient matrix A obtained in the CRBF method is nonsingular and has a smaller condition number than that of the traditional RBF method. Most importantly, the condition numbers of the CRBF discretized matrices are moderate and stable independent of the shape parameter c. Once the coefficients a j and b j are obtained, we get the approximate solution functions u(x) and v(x) for the horizontal and vertical displacements, respectively. In addition, the stresses can be computed by

Numerical Examples
In this section, two numerical examples are present to show the advantages of the newly developed CRBF method over the traditional RBF method for the solution of the elastostatic problems. In actual computations, we choose two coupled radial basis functions, i.e., the CMQ and CIMQ RBFs. Correspondingly, two traditional radial basis functions, i.e., the MQ and IMQ RBFs, are adopted. In the first example, we consider a cantilever beam carrying an end load. In the second example, we study an infinite plate with a central hole loaded by traction at infinity in the horizontal direction. Both problems are well-known and have analytical solutions. In order to evaluate the performance of the CRBF method as well as the RBF method, we use the following relative errors with respect to the displacement and the stress: are the numerical solutions of the displacement and the stress, respectively, at the ith node, and u e i and σ e i are the corresponding analytical solutions at the ith node.

Cantilever Beam.
For the first example, we consider the cantilever beam of length L, height D, and unit width, which is subjected to traction P at the right free hand; see Figure 1. In particular, a coordinate system with the y-axis centered at the midplane of the beam is used. us, the upper and the lower surfaces of the beam are located at y � ± (D/2). e analytical solution of the displacement is [23] where E is Young's modulus, ] is Poisson's ratio, and I is the moment of inertia. e stresses corresponding to the above displacements are e cantilever beam problem is chosen for numerical example because it is a widely known and well-understood. Many meshfree methods use this problem as a standard test problem; see [26,[29][30][31][32] for examples. However, in many existing works, they did not consider the case that the boundary conditions are necessary to match the exact solutions. erefore, the conclusions based on numerical errors computed using these solutions are questionable. Augarde and Deeks pointed out this fact in [25]; see also [29]. In fact, the displacements given in (19) are not an exact solution of the plane stress equations only if the load is distributed parabolically (the third equation in (20)) and if the essential boundary conditions are applied at x � 0 according to (19). To fairly compare the numerical errors computed by different solution methods, Simonenko et al. suggested using Dirichlet boundary conditions given by (1) in all boundaries of the beam and comparing the numerical solutions with the exact solution for displacements (19) and stresses (20) [29]. In this paper, we adopt the same strategy when computing numerical solutions with the traditional RBF methods and the CRBF methods. Some parameters are taken as L � 12, D � 2, E � 1000, ] � 0.3, P � 1, and I � D 3 /12.

Relationship between Relative Error and Shape
Parameter. We first show that the CRBF methods have high accuracy and the shape parameters have almost no effects on the accuracy of the CRBF methods. To better show these results, we consider both regular and irregular nodes distribution. In Figures 2(a) and 2(b), we plot the regular 73 × It can be found from these figures that the CRBF methods always have better displacement accuracies than the traditional RBF methods. When the optimal shape parameters c of the MQ and IMQ RBFs are obtained by minimizing the relative errors r σ , the traditional RBF methods seem to perform a little better than the CRBF methods. However, the optimal shape parameters c for r σ may not be optimal for r u . erefore, it is a big problem for the traditional RBF methods to find optimal shape parameters. Most importantly, both r u and r σ of the CMQ and CIMQ CRBF methods are almost independent of the shape parameter.
ese numerical results show that the CRBF methods not only have better accuracy but also show better stability than the traditional RBF methods when solving the Cantilever beam problem.

Relationship between Condition Number and Shape
Parameter. Condition number of the discretized matrix is an important index for a solution method. High condition number may result in instability and sometimes make the solution unacceptable. In this subsection, we show numerical results about the condition numbers of the RBF and CRBF discretized matrices. In particular, Figures 4(a) and 4(b) plot the condition numbers of four discretized matrices for regular and irregular nodes distribution with respect to different shape parameters c (varying from 1 to 100), respectively.
From these two figures, we can see that the condition numbers of the MQ and IMQ RBF methods change greatly. When the shape parameter c is very small, the condition numbers of the MQ and IMQ RBF methods are not large. However, in this case, we see from Figure 3 that both displacement and stress relative errors are not good. When the shape parameter c becomes large, the condition numbers of the MQ and IMQ RBF methods become large, too. From Figure 4, we also see that the CRBF methods produce relatively well-conditioned solution systems, whose condition numbers are moderate and almost independent of the shape parameter. ese results show again that the CRBF methods show better stability and robustness than the traditional RBF methods when solving the Cantilever beam problem.

Convergence Analysis, Computational Results, and
Simulation. In this subsection, we study the convergence of the CRBF methods when solving the cantilever beam problem. Here, only four regular nodes distributions, i.e., 13 × 3 (13 nodes in the x direction, 3 nodes in the y direction), 25 × 5, 49 × 9, and 73 × 13, are considered. e computational results are listed in Table 3. e convergence curves of displacement and stress with different nodes distribution are plotted in Figure 5. In Table 3, c opt stands for the optimal shape parameter, Cond(A) represents the condition number of the matrix A, and CPU denotes the elapsed CPU times. Here, the optimal shape parameters c opt are obtained by minimizing the relative errors of displacement r u . In Figure 5, h is the maximum size of nodes arrangement.
From Table 3 and Figure 5, we see that the RBF methods do not show convergence for solving the cantilever beam problem since relative errors r u and r σ do not slow down as number of nodes increases. However, the CRBF methods exhibit very nice convergence properties. For small sizes, the accuracy of the traditional RBF methods is better than that of the CRBF methods. However, for large sizes, the CRBF methods show better computational results.
To further show the computational results, we plot both analytical solution and numerical solution of the CRBF methods with the CMQ RBF in Figures 6 and 7. More specifically, Figures 6(a) and 6(b) plot the displacement for regular 73 × 13 nodes and irregular 949 nodes, respectively. Figure 7(a) plots the normal stress σ xx and Figure 7(b) plots the shear stress τ xy for regular nodes distribution. From these figures, we see that the numerical solutions of the CRBF methods are in good agreement with the analytical solutions.

Plate with a Hole.
For the second test problem, we consider an infinite plate with a central circular hole of radius a subjected to a unidirectional tensile load of σ 0 at infinity in the x direction (see Figure 8(a)). is is a typical plane stress problem and has been used as a standard test case to assess the accuracy of different meshfree methods [29,30,32]. By taking advantage of its symmetry, only a quarter of the model is considered in the analysis (see Figure 8(b)). e second test problem can be discussed in polar coordinates or Cartesian coordinates system. Here, we only consider the Cartesian coordinates system. e analytical solutions of the displacements are where e corresponding analytical solutions of the stresses are In actual computations, we solve the problem with σ 0 � 1, E � 1000, L x � 5, L y � 5, a � 1, v � 0.3, and Dirichlet boundary conditions [29].

Relationship between Relative Error and Shape
Parameter. We first study the influence of shape parameters in the CRBF method on the numerical results. To this end, we take two typical nodes distributions, which are plotted in Figure 9. In particular, Figure 9(a) plots 289 nodes and Figure 9(b) plots 625 nodes. As can be seen from these figures, nodes near round holes are denser than elsewhere as stress concentration will occur around the round hole. In From the numerical results plotted in Figure 10, we see that for the traditional RBF method the relative errors of both displacement and stress vary greatly, while for the CRBF method the relative errors of both displacement and stress remain almost constant when the shape parameter varies from 1 to 100. e optimal parameters of the MQ and IMQ RBFs seem small. e relative errors of displacement of the CRBF method are slightly worse than those of the MQ RBF method and are much better than those of the IMQ RBF method even if the optimal parameters are used. However, in terms of the relative errors of stress, the CRBF methods are always much better than the traditional RBF methods. ese numerical results show again that the CRBF methods not only perform more robustly but also have higher accuracy than the traditional RBF methods for solving the elastostatic problems.

Relationship between Condition Number and Shape
Parameter. In this subsection, we study the condition numbers of the CRBF and RBF discretized matrices arising from the elastostatic equations, in particular the plate with a central hole problem. In Figure 11, we plot the curves of the condition numbers. More specifically, Figure 11(a) shows the relationship between the condition numbers of the discretized matrices and the shape parameters for the case of 289 nodes distribution. Figure 11(b) shows the same items for the case of 625 nodes distribution.   From Figure 11, we see that the condition numbers of the RBF discretized matrices vary greatly with respect to the shape parameter and for most cases are larger than 10 20 , which are very ill-conditioned matrices, while the condition numbers of the CRBF discretized matrices almost remain constant even if the shape parameter c varies from 1 to 100. For the cases of 289 and 625 nodes distribution, the condition numbers of the CRBF discretized matrices are near 10 14 and 10 15 , which are much smaller than those of the RBF discretized matrices.
is, again, shows that the CRBF method is more stable than the traditional RBF method for solving the elastostatic problems.

Convergence Analysis, Computational Results, and
Simulation. In this subsection, we further study the CRBF methods for solving the elastostatic equation, in particular the plate with a central hole problem, and show some computational results. For convergence studies, three different nodes distributions with 289 nodes, 625 nodes, and 1089 nodes are considered here. We plot the 289 nodes and 625 nodes in Figure 9. All these three nodes distributions can be found in [35]. e computational results are present in Table 4. Note that the optimal shape parameters listed in Table 4 are obtained by minimizing the relative errors of displacement. To better show the convergence properties, we plot the convergence curves of the four discussed methods in Figure 12. From these results, we find that both the RBF methods and the CRBF methods exhibit high convergence for solving the second test problem. e optimal shape parameters of the MQ and IMQ RBFs are small. In terms of relative errors, the CRBF methods seem slightly better than the traditional RBF method. However, from the aspect of condition numbers of discretized matrices, the CRBF methods have obvious advantages. Most importantly, there is no need to choose the optimal shape parameters for the CRBF methods. e displacements computed by the CRBF method at nodes are plotted and compared with the exact solutions (see Figure 13(a) for 289 nodes distribution and Figure 13

Conclusions and Outlooks
As is known to all, the traditional RBF method has the advantages of high accuracy and flexibility. When solving the partial differential equations, the RBF method is a truly meshfree method. However, the accuracy and stability of the RBF method strongly depend on the shape parameter. How to choose the optimal parameter for the RBF method is an ongoing challenge problem. A coupled radial basis function (CRBF) method, which was developed by coupling the infinitely smooth RBF with the conical spline, was proposed recently in [23]. Two numerical examples arising from the Poisson equation and the heat transfer equation show that the CRBF method completely overcomes the troublesome issue of the optimal shape parameter. In this paper, we further extend the new methodology of the CRBF method to solve the elastostatic problems. Numerical results show that solutions of both the displacement and the stress (the first derivative of displacement) as well as the condition numbers of the CRBF discretized matrices are almost independent of the shape parameter. Moreover, the CRBF method achieves better accuracy than the traditional RBF methods even with optimal shape parameters. e results of this study further expand the application scope of the CRBF method. Although the CRBF method was successfully applied to solve the elastostatic problems, only global collocation scheme is considered in this paper. e discretized matrices are dense, which is not suitable for solving large scale problems. In addition, only two 2D numerical examples are used to verify the efficiency of the CRBF method; no theoretical result is provided to guarantee the accuracy, convergence, and stability of the CRBF method. Meshfree methods have great advantages in solving 3D problems. erefore, future works should focus on 3D problems, more applications of the CRBF method, studying the local CRBF methods, combining the CRBF with other meshfree methods, and providing some theoretical results.

Data Availability
e data used to support the findings of this research are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest. 14 Mathematical Problems in Engineering