New collocation path-following approach for the optimal shape parameter using Kernel method

The goal of this work is to develop a numerical method combining Radial Basic Functions (RBF) kernel and a high order algorithm based on Taylor series and homotopy continuation method. The local RBF approximation applied in strong form allows us to overcome the difficulties of numerical integration and to treat problems of large deformations. Furthermore, the high order algorithm enables to transform the nonlinear problem to a set of linear problems. Determining the optimal value of the shape parameter in RBF kernel is still an outstanding research topic. This optimal value depends on density and distribution of points and the considered problem for e.g. boundary value problems, integral equations, delay-differential equations etc. These have been extensively attempts in literature which end up choosing this optimal value by tests and error or some other ad-hoc means. Our contribution in this paper is to suggest a new strategy using radial basis functions kernel with an automatic reasonable choice of the shape parameter in the nonlinear case which depends on the accuracy and stability of the results. The computational experiments tested on some examples in structural analysis are performed and the comparison with respect to the state of art algorithms from the literature is given.


Introduction
Most of the non-linear problems encountered in physical phenomena are resolved using numerical methods, which are based on a subdivision of the domain or a mesh, such as finite difference, finite volume and finite element methods. Creation of these meshes in three dimensions is not an easy work and even impossible virtually for higher dimension. This is where the meshfree or meshless methods will play a crucial role by adjusting a problem involving many space dimensions into one that is virtually one-dimensional. This possible because the meshfree are radially symmetric in nature. The origin of meshless methods dates back to the late seventies, but their development remained very limited until the early nineties. The first meshless method is Smooth Particle Hydrodynamics SPH, which was introduced in 1977 by the British astrophysicist Lucie [23] and by the Australian mathematicians Gingold and Monaghan [15]. This method has been used to solve astrophysical problems, fluid dynamics problems [3,26,27] and problems in solid mechanics [24]. Subsequently, several developments and improvements were proposed [3,5]. Then, the diffuse elements method was proposed by Nayroles et al. [30] which consists to use a Moving Least Square (MLS) approximation with a Galerkin type discretization. Another method, called Element Free Galerkin (EFG) method, was developed in 1994 by Belytschko et al. [4]. A year later, the Reproducing Kernel Particle Method (RKPM) was developed [25]. Another somewhat recent method is called Natural Elements Method (NEM) [10,37]. This method uses shape functions based on geometric constructions such as the Voronoi diagram and the Delaunay triangulation. Currently, work on meshless methods are devoted to improving existing versions and adapting them to specific problems [1,38].
In the literature, there are recent research papers related to collocation and interpolation methods, these methods have successfully and efficiently applied to solve a wide variety of problems which describe various physical phenomena in applied science and technology. As examples for solving generalized Black-Scholes partial differential equation, Kadalbajoo et al. [20] used a cubic B-spline collocation method, Mohammadi [29] a developed a Quintic B-spline collocation approach and Roul et al. [33] described a sixth order numerical method and its convergence. For solving Bratu-type and Lane-Emden problems, Roul and Thula [34] developed a fourth-order B-spline collocation method. Also in the context to solve a class of Lane-Emden singular boundary value problems which describe several phenomena in theoretical physics and astrophysics. In the field of electrohydrodynamic flow of a fluid in a circular cylindrical conduit, Roul [35] presented a fourth-order non-uniform mesh optimal B-spline collocation method for a strongly nonlinear singular boundary value problem. We also used collocation and interpolation methods in our previous work for solving nonlinear elastic and elastoplastic problems [1] and viscoplastic problems [1,28,38].
In this work we are interested to another class of meshless methods developed in 1990 by Edward KANSA [18,19], from the University of California, for solving partial differential equations. This is Radial Basis Function (RBF) method which is based on the interpolation theory using RBFs. It has been used successfully in the following works [17,31,32]. RBF-based method offer numerous advantages, such as no need for a mesh or triangulation, dimension independence, simple implementation in multivariate scattered data approximation [6] and it is becoming a best choice as a method for the numerical solution of partial differential equations [13,14]. For these reasons, we have therefore chosen to use the RBF collocation method in our work.
The objective of this work is to develop a novel approach based on the RBF collocation method and a shape parameter search algorithm that determines the optimal shape parameter with a good precision of the results whatever the points distribution. In this contest, we develop an effective strategy coupling the high order meshless approach, based on the strong form RBF approximation, and an optimal value search algorithm of the shape parameter. The proposed approach is automatic and the choice of the shape parameter depends on the accuracy and stability of the results. For validation, we consider three examples in the context of the structures analysis taking into account the geometric nonlinearity.
The outline for the rest of this paper is as follows: Sect. 2 presents the collocation technique based on RBF interpolation. Section 3 shows a reminder on the considered nonlinear problem. Section 4 presents the collocationpath following approach based on RBF interpolation. Section 5 explains the coupling strategy of the path-following method with the local RBF approximation. Section 6 shows the proposed strategy for funding the optimal shape parameter for nonlinear problems. Some numerical applications are given in Sect. 7, to show the performance of the proposed strategy, and it ends with a conclusion in Sect. 8.

Collocation technique based on RBF method
In this section, the local version of the classical meshless radial basis function collocation (Kansa) method is introduced. For a recent chronological and historical development scheme of RBF methods for solving partial differential equations, we refer the reader to the paper [2] . Consider a vector x in ℝ d and the Euclidian norm ‖ ⋅ ‖ 2 on ℝ d , the radial basis functions have the form (‖x − x j ‖ 2 ) which assumed to be strictly positive definite. Using a set of arbitrarily scattered nodes localized in a local support domain ⊂ ℝ d , and assigning to each point x a set of nearest neighbors points x j included in the support domain, the RBF approximation can be written as: We focus on RBFs which are mostly identified on the basis of smoothness i.e. infinitely differentiable and contain free parameters c i called the shape parameters. We will be mostly interested in the multi-quadrics and the exponential functions Other popular used RBFs or compactly supported functions will not play a role in this paper. The coefficients = [ 1 , ⋯ , N ] T are unknowns yet to be determined by enforcing the interpolation constraints as follows: which leads to N linear equations which can be expressed in the matrix form as where R ij = (‖x i − x j ‖), i, j = 1 ⋯ N is called interpolation matrix and defined symmetric and positive for all the neighbors, and In order to ensure the existence of R −1 and a well-conditioned R, we must choice a good shape parameter c. Unfortunately, there is no theoretical best value of c, and it has to be determined by numerical experiments to stabilize the solution which is very difficult and expensive in terms of calculation time. However, it is very difficult to find the best optimal value of c. In this paper, we are interested to develop a novel approach that determines the optimal parameter c with a good precision of the results whatever the points distribution.
The approximation of the function s(x) can be written in the following form:

Nonlinear elastic problem statement in strong form
In this section, we present the equations of equilibrium in strong formulation for a nonlinear elastic structure (geometrical nonlinearity) in large deformations [1]. It is assumed that the displacements and forces imposed are proportional to a single scalar parameter called the load parameter, then the problem to be solved is written as follows: where T is the first tensor of Piola-Kirchhoff, S represent the second tensor of Piola-Kirchhoff, D is the fourth order elastic behavior tensor in the case of plane stresses, is the domain occupied by the structure where U and F are the boundary associated to an imposed displacement U d and applied loading F respectively with N is the normal applied toward the outside of the boundary F .
The problem (7) can be written in the following matrix form: w h e r e :

Numerical path-following continuation
Our aim in this section is to review briefly numerical methods for solving nonlinear continuation (or homotopy) problems in the form: where G ∶ ℝ N+1 → ℝ N is a smooth function or operator, u ∶ ℝ N → ℝ is a real-valued function which can represents the "solution" (i.e., flow field, displacements, etc.) for various values of a real scalar physical parameter ∈ ℝ (i.e., Reynold's number, load, etc.). In geometric terms, the Eq. (9) defines a one-dimensional implicit curve in ℝ n+1 under appropriate regularity and consistency properties of the mapping G. The task then is to find the solution for some -intervals, that is, a path of solutions, [u( ), ]. Numerical continuation is the procedure of solving such systems of nonlinear equations. The key idea of continuation is to follow the underlying curve of solutions. The trivial approach for numerical implementation is called a parameter continuation [16,22] which draws a solution path by starting at a point = init by repeatedly increasing and solving G(u, ) for u until the desired value of is reached. While parameter continuation is intuitive, simple and seems to be a reasonable method for solving G(u, ) = 0 , it fails at points (u, ) that violate the assumptions of the implicit function theorem where the Jacobian G u is singular.
For this purpose, we will focus in this paper on a class of continuation based on parameterising the solution branches by a pseudo arc-length, say [u(s), (s)] to avoid the problems of simple parameter continuation. In fact, these pseudo arc-length continuation methods [11,12,21], are that most singular points on the solution branches can be handled without much difficulty. We get ride of these difficulties by not parameterizing the solution u by . Instead, we parameterize the bifurcation solutions using an pseudo arc-length parameter s, and specify how far along the current solution branch we want to proceed. To be more precise, we let s denote the arclength parameter along a solution arc (u(s), (s)) . We compute the tangent [u(s 0 )],̇(s 0 ) (where the dots denote differentiation with respect to s) of a known solution at s = s 0 . So The smooth homotopy path u(s), (s) must satisfy the differential equation: obtained from differentiating G(u, ) = 0 with respect to s. In addition, the norm is the Euclidean norm and s is arclength, we impose the arclength condition: The kernel of the Jacobian has exactly two vectors of unit norm which correspond to the two possible directions of traversing the curve. In general, one will wish to traverse the solution curve in a consistent direction. In order to specify the orientation of traversing, an additional equation for closing the system must be added to define the path parameter a to G(u, ) = 0 so that the number of equations equals the number of unknowns. We obtained a nonsingular Jacobian for the reformulated problem. The proof of this assertion is based on the implicit function theorem and can be found in many references, see for instance [11,12,21]. Hence, we work with the augmented equations: The normalization equation N = 0 is an linearization of (11) and require the new solution u(s) and (s) to satisfy This equation tells us that it forces the new solution to lie on a hyperplane perpendicular to the tangent vector to the solution curve at s 0 at a distance (s − s o ) (see Fig. 1).

Coupling of the path-following method with the local RBF approximation
The mixed unknown vector of problem (10) is searched in our modeling in the form of a Taylor series development from a known solution point ( T 0 , S 0 , 0 , U 0 and 0 ) as follows: Let the path parameter "a" denote the arclength along a solution arc (u(a), (a)) which depends smoothly on this parameter and x = (u T , ) T . As mentioned in the previous section we need to close our system by adding an equation and we have a chance of obtaining a non singular Jacobian for the our problem statement. Several choices are possible and the most used [1,7] is a pseudoarc-length parametrization. It corresponds to the projection of the pair (U − U 0 , − 0 ) , ie, the displacement and the load parameter, on the tangent direction U 1 , 1 which allows us to have an additional condition as follows: To develop a pseudo-arclength continuation method, we assume that x depends smoothly on a. Then one can differentiate G(x) = 0 in addition "a" is arclength and the norm is the Euclidean norm. Therefore, we get the following equations depending of the order p: By injecting the developments [see Eq. (14)] into the nonlinear problem (8) and identifying the coefficients having the same powers of the parameter a, we obtain a sequence of linear problems given by: Problem at order 1: The matrix Ŝ 0 contains the stress components of the starting solution of each order defined as follows: In the problems (17) and (18), the principle unknown U p is approximated using the local RBF method which is given by: where [ (X )] is the matrix of the RBF shape functions and U p is a vector which contains the new unknowns at order p.These terms are defined at each point X by the following expressions: By injecting the approximations of U p in p , then p in S p and then S p in T p , and after an assembly technique, we obtain a compact problems as follows: Order p for 2 ≤ p ≤ N order : where [K T ] is the stiffness tangent matrix evaluated at the starting point ( T 0 , S 0 , 0 , U 0 ). Using the continuation procedure, we consider T (a max ) , S(a max ) , (a max ) , U(a max ) as the new starting point of the new solution branch. The convergence radius a max is defined as in the reference [1,7]:

Estimation of the optimal shape parameter
For the RBF kernel methods, there is no rule governing the choice of the shape parameter c, this coefficient can be selected from numerical tests of each RBF function to stabilize the solution which is very difficult and expensive in terms of calculation time. In the present work, we develop an algorithm to determinate quickly and automatically the best value of c. The idea of the strategy for determining the optimal shape parameter c by combining the RBF method, high order algorithm Taylor expansion and numerical continuation approach adapted to the algorithms developed in the references [8,39]. In this strategy, we consider the new definition of the shape parameter c = ds where is a coefficient and ds is the average influence domain. Note that ds = 1 N ∑ N i=1 d i where d i is the distance between the ith point and its nearest neighbors.
The principle for determining the optimal shape parameter is presented in Fig. 2. The the relative error ( error( k n ) ) of the displacement at order 1 decreases as the shape parameter increases at the beginning. Continuously increasing , the relative error becomes larger and the proposed algorithm gives incorrect result. The goal is to minimize the relative error of the displacement at order 1 of the high order mesh-free algorithm.
To be more precise, we propose a simplified presentation of our adapted algorithm illustrated in Fig. 3 where the detail of the part 1 to part 3 of the algorithm is depicted separately in Figs. 4, 5 and 6. The master algorithm shows then the proposed strategy for estimating the optimal shape parameter. This first order error estimator allows us to ensure a well-conditioned tangent matrix [K T ] at higher order algorithm and by iteration which is the same used in the higher orders. For this reason, in the minimization of the relative error, we are

Numerical analysis
In this application, we study with the proposed high order algorithm the effect of the shape parameter on the solution of the non linear elastic problem. The accuracy of this numerical solution is controlled by the FEM method as a reference solution.

Study on the approximation accuracy according to the shape parameter values of a bi-dimensional structure in tension
In this study, we consider the geometrical nonlinearity of a bi-dimensional elastic structure in tension using a 2D plate of geometry (100 mm × 50 mm) . This plate is fixed at x = 0 and submitted to an imposed load F at x = L ; with F = 1 (see Fig. 7). The parameters of the high order algorithm are the truncation order p = 15 and the tolerance parameter = 10 −8 .
As the exact solution of this example is not available, so the accuracy of the numerical solution will be computed using double mesh principle. In For any fixed value of number of points Np, the maximum point-wise error G Np will be calculated by: is the numerical solution at the same point x i which is obtained by increasing the original mesh with 2Np number of points. The order of convergence is defined as : Table 1 depicts the maximum point-wise error and order of convergence for = 10 4 and several values of number of points Np. From the results, it can be observed that as the grid size h decreases, the maximum absolute errors decrease, which shows the convergence to the computed solution.
We performed a simulation using the proposed strong form RBF approximation, distributed points number 231 and MQ function with q = 0.5 . Figure 8 shows the deformed configuration for = 179,690.
The evolution of load F versus displacement u is plotted in Fig. 9 at point P (see Fig. 8) for the proposed algorithm with q = 0.5 and FEM based on Newton-Raphson method. After several tests using different values of the shape parameter , the best result is obtained by = 10. Figure 10 shows the evolution of load F , versus displacement u at point P, obtained by the proposed algorithm with automatic choice of the shape parameter , by the strong form MLS approximation [1] and by FEM.
To show the performance of our algorithm, we are looking for precision in the interval [38.2 mm, 40.2 mm] . For the best choice of the shape parameter = 10 , the relative error is approximately equal to 3% (see Fig. 11a). Fig. 6 Part 3 of the resolution strategy for estimating the optimal shape parameter (see Fig. 3)

Fig. 7
Bi-dimensional elastic plate subjected to a tensile F  On the other hand, the relative error does not exceed 1% (see Fig. 11b) if the proposed algorithm is coupled with the strategy for determining the optimal choice of the shape parameter . So, we can get good results with a automatic choice of .
To give an idea on the functioning of the search algorithm for determination of the optimal shape parameter, we plot the evolution of relative error variation versus shape parameter in Fig. 12a. This figure represents the first estimate (k = 0) of the optimal shape parameter which gives a value between 35 and 36.
So, in the second estimate (k = 1) , we are looking for precision in the interval [35,36] or the accuracy of optimal to the first decimal place as shown Fig. 13 Figure 14 presents the evolution of the optimal shape parameter optimal with respect to number of steps for the new algorithm which shows that optimal varies versus the number of steps. Therefore the parameter can vary slightly with respect to the increase in loading or deformation.
The conditioning of the system matrix and the accuracy of the RBF method are affected by the shape parameter. The shape parameter must not be too small to ensure that the system matrix to be well-conditioned. However, to obtain a good accuracy for the RBF method a small shape parameter is required, but this choice results in the system matrix being ill-conditioned. Therefore, the best precision and the best conditioning cannot be achieved at the same time, so we must search a compromise. This is known as the Uncertainty Principle [36]. The stability result of the proposed method can be determined so that the resulting system matrix has a condition number in the range [10 13 , 10 15 ] [9]. This condition number is given by: where the R is the momentum matrix assigned to the MQ RBF interpolation. Figure 15 presents the condition number in each point of the domain for the MQ RBF interpolation. We notice that the condition number is in the range [10 13 , 10 15 ] . On the other hand, the evolution of the condition number (see Fig. 15) depends of the optimal shape parameter as shown  (a) Zoom of Fig.8 (b) Zoom of Fig.9(right) in Fig. 14. So the stability result of the proposed method can be verified by seeking a compromise between precision and conditioning. Now we performed a simulation using the strategy for determining optimal choice of the shape parameter and the Inverse Multi-Quadric (IMQ) function as radial basis function ( q = −0.5 ), we give the displacement at point P using the present algorithm with automatic choice of shape parameter in Fig. 16. Figure 17 represents the evolution of the optimal shape parameter optimal versus the number of steps. We can notice the same remark that the increasing the load or the deformation can result a small increase in the value of the shape parameter .
For the radial basis function given by Eq. (2), for a given value of q, the new algorithm with automatic choice of the shape parameter seeks the good precision. For example we can test the RBF function with q = 0.25 , the proposed present algorithm can determinate automatically the adapted shape parameter using the proposed strategy given in Fig. 3. Figure 18 shows the evolution of load with respect to displacements U with automatic choice of compared to the FEM at point P. Figure 19 represents the evolution of optimal shape parameter versus number of steps.
We can test other function as exponential [see Eq. (3)]. Figure 20 shows the evolution of load with respect to displacements u with automatic choice of and it is compared to FEM at point P. Figure 21 represents the evolution of optimal shape parameter optimal versus number of steps. We note that it is constant in the case of exponential function.

Example of a bi-dimensional structure in bending
We consider a thin elastic plate clamped at the left end and subjected at the right end to a bending load. The used mechanical and geometrical characteristics of the plate are: Young's modulus E = 210 GPa , Poisson's ratio = 0.3 , length L = 100 mm , width H = 10 mm , bending load F = 1 MPa . Regarding the high order algorithm, we use the same parameters of the previous example with a distributed points number N = 729 . Figure 22 shows the initial and deformed configurations of plate in bending. Second estimate of the optimal shape parameter α optimal with n=4 (α n-3 1 , error1) (α n-2 1 , error2) (α n-1 1 , error3) α optimal =α n-3 0 + =35.3 (α n 1 , error4)   Fig. 19 Evolution of optimal shape parameter optimal with respect to number of steps using a MQ function with q = 0.25 Figure 23 represents the evolution of displacements u and v of point P versus loading parameter using MQ function where q = 0.5 . We can see that the calculation by the high order algorithm based on FEM stops because the convergence radius of this high order algorithm tends to zero which shows the advantage of the high order algorithm based on RBF to simulate this kind of problems. Figure 24 represent the evolution of the optimal shape parameter optimal versus the number of steps. As in the case of a bi-dimensional structure in tension, the parameter can vary with respect to the increase in loading or deformation.
Using exponential function as RBF function, Fig. 25a, b represent the evolution of displacements of point P versus the loading parameter and the evolution of the optimal shape parameter optimal versus the number of steps, respectively. As the example of traction, we also note that it is constant in the case of exponential function.
As shown in this paper, the selection of the shape parameter has a key influence on the numerical accuracy and efficiency in the meshless modelling based on RBF interpolation in strong form. The optimal value of the shape parameter is highly related to the adopted nodal scheme for discretization and the problem itself, and it is often determined by trial calculation traditionally. In what follows, the objective is to compare the computation time for finding good results between the proposed approach coupled with the algorithm for automatic selection of shape parameter and without using this algorithm. In this comparison, we are looking for a good result by testing the values of the shape parameter c ∈ [1 ds, 40 ds] with a precision of one digit after the decimal point. Note that using the proposed approach with a fixed c, 400 tests with c = [ds, 1.1 ds, 1.2 ds, 1.3 ds....40 ds] must be carried out and compare between them to choose the good value of c, and then good results. The CPU time required to reach the load * using the two algorithms is reported in Table 2. We can observe that the corresponding time for the proposed algorithm with automatic selection of c optimal is much less than the time needed for the proposed algorithm with a fixed c, which confirms the effectiveness of the proposed approach.

Conclusion
The main goal of this paper is to establish a new numerical collocation path-following approach using RBF kernel method and a high order algorithm based on Taylor series and homotopy continuation method. The second objective and an outstanding research was to suggest an optimal choice of the shape parameter in the RBF interpolation in view of applications for the simulation of nonlinear problems. For this purpose, we considered as paradigm nonlinear elasticity problems with large deformation which includes large strain, displacement, and rotation. Our choice has at least two reasons: the simulation of large structural deformations with the finite element method poses several challenges and there is an abundant literature related to this problem where many techniques can employed to try to avoid the element distortion or split the total amount of force or displacement to be applied into smaller parts or change the original problem either by remeshing with more or less elements. The problem still does not converge or a long time consuming process even if take into consideration all suggestions. The proposed RBF approach has been used successfully in conjunction with continuation procedures on the model problem using as arclength the load parameter and Taylor series development of known solution points namely the first and second tensor of Piola-Kirchhoff, the displacement and the normal applied toward the outside of the boundary. After injecting different approximations, we have then obtain a compact problem up to order p. Each nonlinear system can be then solved by computing the iterative solution of linear systems of equations determined by Jacobian matrices associated with the nonlinear system of equations. The pseudo-arclength parametrization continuation technique is very powerful and efficient procedure applied in the vicinity of turning points to overcome the singularity of the Jacobian matrix to improve convergence properties when an adequate starting value for an iterative method is not available.
Since many radial basis function kernel methods contain free shape parameter which plays an important role for the accuracy of the method, we develop an algorithm to determinate quickly and automatically the best value of this parameter. Numerical results show that the proposed adaptive algorithm with automatic selection of the shape parameter is effective and provides reasonable more accuracy compared to the same algorithm with a fixed   and given shape parameter. The limitation of this adaptive algorithm is the fact there is no guarantee to find the global minimum, but it can get stuck in a local minimum. The values of the optimal shape parameters depend on the number of collocations, the shape of computation domains and the distributions of nodes. The possible perspective in this direction is to realize a comparative study with different algorithms from literature for the optimal choice of this free parameter in the RBF techniques depending on the accuracy and the precision. Our guess is that this study will lead to confirm that our algorithm is very efficient as a compromise between the accuracy and stability of RBF approximation before the interpolation matrix becomes overly ill-conditioned. In addition, for structure modeling in large deformation, this new adaptive algorithm gives a good result with a good precision compared to FEM. The comparison with FEM shows that the relative error does not exceed 1% . The proposed strategy of this adaptive algorithm is currently work in progress to other application fields such as the mixing process observed in Friction Stir Welding.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Human and animal rights This article does not contain any studies with human participants or animals performed by the author.

Informed consent
No individual participants were included in the study so there are no subjects for which informed consent requirements arise.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.