Solving biharmonic problems on irregular domains by stably evaluated Gaussian kernel

Abstract The paper extends recently developed idea of stable evaluation of the Gaussian kernel. Owing to this, the Gaussian radial basis function that is sensitive to the shape parameter can be stably evaluated and applied to interpolation problems as well as to solve differential equations, giving highly accurate results. But it can be done only with grids being the Cartesian product of sets of points, what limits the use of this idea to rectangular domains. In the present paper, by the association of an appropriate transformation with the mentioned method, the latter is applied to solve biharmonic problems on quadrilateral irregular domains. As an example, in the present work this approach is applied to solve bending as well as free vibration problem of thin plates. In the paper some strategies for the implementation of the boundary conditions are also presented and examined due to the accuracy. The numerical tests show high accuracy and usefulness of the method.


Introduction
In recent years a significant development of the meshless methods can be observed. Many formulations of the methods have been developed [1][2][3]. The main their advantage is the flexibility in the discretization of the domain -scattered nodes can be used to this end. It simplifies preprocessing procedures and the use of adaptation algorithms, providing powerful tool for analyzing complex-shaped engineering structures. Some of these methods make use of the radial basis functions (RBFs) to form the interpolant (sought solution), since the RBFs are found to be very useful in problems of scattered data approximation [4][5][6].
The RBFs fall into a more general category of functions, called kernels. Most of RBFs contain a constant parameter, which controls their flatness. The value of the shape parameter has significant influence on the stability of the solution process and finally on the accuracy of the results. Some information on the issue can be found in several papers [7][8][9][10].
Unfortunately, the values that theoretically should lead to more accurate results make the system illconditioned, what causes numerical difficulties in the solution process and finally inaccurate results. A significant effort has been done to find an optimal value of the parameter [11][12][13][14], but most of the approaches are a kind of tradeoff between stability and accuracy.
Another group of the approaches aims to find a new stable basis spanned over the space created by the RBFs.
In [15,16] the new basis is obtained by QR or SVD factorization of kernel matrix. In [17] another interesting approach is presented. It differs from the mentioned ones by the fact that the kernel matrix, which is prone to be ill-conditioned, is not formed. In order to obtain the stable basis appropriate factorization of the kernel matrix is done on the base of its approximation following from Mercer's theorem. Thus far, this method has been successfully applied in some approximation problems and in solving differential equations [18][19][20]. The results obtained show that the instability due to the shape parameter is removed and the parameter can be tuned according to the theoretical findings, giving very accurate results.
The method gives promising results, but its variant in [18] can be applied in higher dimensions only on grids that are the Cartesian product of one-dimensional sets of points. It significantly limits the applicability of the method -the scattered nodes cannot be further used as in typical meshless techniques. It makes a serious problem in solving tasks on irregular domains that are often encountered in engineering problems.
To preserve the stable solution and allow the method to be useful for problems defined on irregular domains, in the present paper the method developed in [17] and modified in [18] is appropriately extended. To this end, the method is associated with a domain transformation tech-nique. The latter takes advantage of the blending function interpolation used to transform the quadrilateral regions [21]. So far, this technique has been successfully associated with some pseudospectral methods [22,23].
In the paper, the method is used to solve biharmonic equations, which can describe some engineering problems, e.g. associated with analysis of thin plates. Using collocation methods to solve such type of problems one has to overcome an inconvenience following from the existence of multiple boundary conditions. In the paper two strategies to implement such boundary conditions are proposed and examined.
The paper is organized as follows: in section 2 the basic information on the kernel-based methods as well as on the derivation of the stable basis is provided, in section 3 two strategies to implement boundary conditions are proposed, in section 4 the treatment of irregular domains with the aid of the mentioned transformation is explained, while in section 5 the numerical results with some comments are presented. Finally in section 6 concluding remarks are drawn.

Stable basis for Gaussian kernel
The kernel function is understood as a real-valued function of two variables coming from d-dimensional space, i.e. K : Ω × Ω →, Ω ∈ d . Using such a function one can conveniently construct a data dependent basis {︀ K(·, x 1 ), . . . , K(·, x N ) }︀ that are a crucial point of kernelbased methods. The latter are very useful in solving interpolation problems as well as differential equations defined on domains in higher dimensions.

Kernel-based method
Let us consider a boundary value problem in a general form as where L and B denote linear differential operators imposed on the sought function u in the domain Ω and on the boundary ∂Ω, respectively and f, g are known functions. Alternatively one can introduce corresponding eigenvalue problem as follows where λ is an eigenparameter. After discretizing the domain as well as the boundary using scatter nodes x i , i = 1,. . . , N we can construct the sought solution for the differential problem using the kernel basis as where c denotes the vector of interpolation coefficients to be determined. By introducing Eq. (3) into Eq. (1) and by applying collocation procedure one obtains The structures of system matrix K LB and vector F are as follows Similarly, the eigenvalue equation can be discretized by incorporating Eq. (3) into Eq. (2), yielding the generalized algebraic eigenvalue problem where K is the extended kernel matrix of the form In the above objects, there is the distinction between interior nodes x i j , j = 1, . . . , N i and boundary nodes x b j , j = 1, . . . , N b . Note that the total number of nodes is If only the system matrix of Eq. (4) is not singular, the interpolation coefficients can be computed and the approximate solution in the form of Eq. (3) is determined. As one can see, the procedure of the kernel-based method applied to the solution of differential equations is straightforward.
Moreover, the solution is obtained in analytic form, what can be viewed as another advantage of the method. As the kernel functions, the RBFs are willingly used.
The main drawback of the method is high illconditioning of the system (Eq. (4)). When using RBFs this ill-conditioning can be controlled by properly determined shape parameter but then the theoretically estimated high accuracy of the method is not achievable.
In order to maintain the possibility of achieving spectral accuracy, in the present paper another way to overcome the instability is applied. To this end a stable basis for RBFs that has been derived from Mercer's theorem is used.

Derivation of stable basis
According to Mercer's theorem, every positive definite kernel can be represented by an infinite series composed of eigenvalues λn and eigenfunctions φn of associated Hilbert-Schmidt eigenvalue problem [17]. This representation can be put as follows Since the Gaussian RBF e −ε 2 ‖x−z‖ 2 is a positive definite kernel it can be presented in such a way. Generally, if one operates in d-dimensional spaces, d-dimensional eigenvalues and eigenfunctions for Gaussian kernel follow immediately from the univariate ones via tensor product form of this kernel. In most common one-or two-dimensional cases (Ω ∈ R or Ω ∈ R 2 ), Eq. (6) for Gaussian kernel takes the form or e −ε 2 ‖x−z‖ 2 = ∑︁ The eigenvalues and eigenfunctions from Eq. (7) and (8) assume the following form where Hn is the Hermite polynomial of degree n, and are constants including: ϵ -the shape parameter of the Gaussian function, α -a parameter, which parametrizes the weight function in the Hilbert-Schmidt integral operator and the associated inner product.
To take advantage of Mercer's theorem in numerical computation the series has to be truncated. The number of terms, which is taken into account, in relation to number of nodes leads to three general approaches for derivation of alternate basis. These approaches are described in detail in [24]. In the present work the case that the number of terms equals the number of nodes is considered. It was found [20] that this case give a good trade-off between accuracy and economy of computation.
In this case the Mercer series given by Eq. (6) assumes the form where With this assumption the base vector from Eq. (3) can be put as = ϕ(x) T Λ Φ T and the sought solution from Eq. (3) takes the form In above equations matrix Φ is as follows Following the idea of the kernel-based method all we have to do is to introduce Eq. (14) into Eq. (1) or into Eq. (2), dependently on the problem considered, and collocate the governing equation at the interior nodes and the boundary condition at the boundary nodes, what, in the case of the boundary-value problem, yields In above equations matrix Φ LB has the structure Solving Eq. (15) with respect to interpolation coefficients c and evaluating the sought solution (14) one obtains The main conclusion following from Eq. (16) is that the eigenvalue matrix, which is the main source of illconditioning [17], is algebraically removed. The eigenfunctions contained in vector ϕ(x) play the role of the basis for the assumed solution.
To determine the solution one only needs to take the number of eigenfunctions that equals the number of nodes. With the use of these eigenfunctions, base vector ϕ(x) and matrix Φ LB can be formulated and the solution can be constructed using formula (16). Since the main source of ill-conditioning (matrix Λ) is removed the process of the inversion of Φ LB is much more stable than using Gaussian kernel, even for very small values of the shape parameter. It leads to theoretically proved spectral accuracy. Therefore vector ϕ(x) can be treated as a vector of stable basis for assumed solution. The solution described by Eq. (14) can be alternatively put as where b is the vector containing new interpolation coefficients. It has the form b = Λ Φ T · c. The eigenvalue equation (2), written with the use of expression given by Eq. (17), takes the form where matrix Φ is a modification of matrix Φ introduced in Eqs. (13) and (14), and has the structure Solving such defined eigenvalue problem one does not operate highly ill-conditioned matrix Λ, therefore eigenpairs (λ i , b i ), i =1 . . . N i can be stably and accurately obtained.
In higher dimensions, e.g. x ∈ R 2 , the kernel method from [18] associated with stably evaluated Gaussian kernel work only on grids being the result of the Cartesian product. It is a serious limitation comparing to RBF kernels that work properly on scattered nodes.
Moreover, the use of the collocation procedure requires a special treatment of problems possessing multiple boundary conditions, such as biharmonic problems.
In the next sections some approaches are introduced to extend the possibility of the application of the method in biharmonic problems defined on irregular quadrilateral domains.

The use of the stable basis in solving biharmonic problems
Let us introduce a biharmonic problem in the form Since ∆ 2 is the biharmonic operator, two boundary conditions on the boundary are required to make the problem well-posed. Similarly one can consider corresponding eigenvalue problem as where λ is the eigenparameter.
If Ω is a rectangular area, the kernel-based method with the stable basis can be used to obtain analytic, approximate solution. In order to implement two boundary conditions at each boundary nodes two strategies are proposed.

Implementation of boundary conditions -strategy 1
In this strategy the distinction is made between not only interior and boundary nodes but also between nodes adjacent to the boundary, denoted in Figure 1 by crosses. At the interior nodes (x i k -unfilled circle) the governing equation is applied while at the boundary nodes (x b kfilled circle) the first boundary condition is implemented. The nodes denoted by crosses are unused during the collocation procedure. Instead of this the second boundary condition is written at these boundary nodes which are adjacent to the unused points. This gives the same number In this strategy the solution is described by Eq. (16), where vector of base functions ϕ(x) is unchanged, the only difference is in some entries of system matrix Φ LB . The latter has now following form where Note, that nodes x b2 k are a subset of all boundary nodes. The subset contains these boundary nodes that are the projection of the adjacent nodes on the proper boundary, as it is shown in Figure 1. The approach leads to the system matrix of size N × N.
The associated eigenvalue problem is also defined by matrix from Eq. (21) as well as by matrix Φ introduced in Eq. (18).

Implementation of boundary conditions -strategy 2
In this strategy the distinction is made only between interior and boundary nodes as Figure 2 suggests. At the interior nodes the governing equation is applied, while at the boundary nodes two boundary conditions are implemented. To maintain the square system matrix Φ LB , the number of eigenfunctions has to be increased. Instead of N eigenfunctions that correspond to total number of nodes, now their number is increased to N+N b . The vector of basis functions assumes the form The structures and entries of the objects included in Eq. (16) and in the associated eigenvalue problem are similar as previously (Eq. (21)) with the difference that matrix Φ B2 is created with the use of all boundary nodes x b k . The approach leads to system matrix of size (N + N b ) × (N + N b ).

Treatment of irregular domain
The kernel-based method with stable basis can be applied on grids being a result of the Cartesian product. It implies that the domain of the problem has a rectangular shape. When the boundaries are curved the direct application of the method is impossible. To extent the usefulness of the method the advantage of the domain transformation can be taken.
In this approach, graphically presented in Figure 3, one needs a transformation between regular domain in (ξ , η) plane and irregular one in (x, y) plane, what can be gen- To this end, a method that uses so-called blending functions [21] can be applied. In the method the transformation between the regions is done by following interpolation function where s = [x, y] T , s i are parametric curves that represent the curvilinear boundaries and s i are Cartesian coordinates of the corner points of the quadrilateral region in (x, y) plane. If only it is possible to find parametric representations of the curves the method exactly maps computational domain onto physical one.
By the use of this transformation all spatial derivatives with respect to physical variables, contained in the equation, have to be written as derivatives with respect to ξ and η.
With the chain rule of differentiation the transformation of the first derivatives is obtained as where J denotes the Jacobian matrix. Then, the invers transformation is as follows where the invers of the Jacobian matrix has the form In Eq. (26) |J| denotes determinant of the Jacobian matrix. Similarly higher order derivatives with respect to the physical coordinates, in terms of the computational ones, can be obtained. In this way one is able to derive the transformation formulas for all partial derivatives contained in the biharmonic equation. This equation, i.e. Eq.  (21) u, ξξ + D (22) u, ξη + D (23) u,ηη + D (11) u, ξ + D (12) u,η = f Similarly the eigenvalue problem, corresponding to Eq. (20), can be written.
The coefficients contained in Eq. (27) depend on the derivatives of physical variables with respect to computational ones and their explicit form can be found in [22].
Since Eq. (27) physically can describe the deflection of thin plate under load represented by function f, the associated boundary conditions can follow from the type of constrains on plate edges.
Typical homogenous boundary conditions are associated with the clamped edge (vanishing normal rotation to the edge) or with simply-supported edge (vanishing normal moment to the edge). They have the following form in (x, y) coordinates In above equations ϑ is the angle between the normal to the boundary and x axes, while υ is Poisson's ratio. These boundary conditions transformed onto computational region in (ξ , η) coordinates assume the form The explicit form of all coefficients contained in Eqs. (30) and (31) is shown in [22]. Once the transformation between the domains is established one can take advantage of the above equations. Following the procedure presented in section 2, the sought solution, now written as a linear combination of eigenfunctions of Gaussian kernel formed in computational coordinates where ζ = [ξ , η] T should be introduced into Eq. (27) or into the corresponding eigenvalue problem as well as into appropriate boundary conditions (Eq. (30) or (31)). After the collocation procedure, where the chosen strategy for the implementation of boundary conditions is involved, interpolation coefficients b can be found, giving the analytic approximate solution.

Numerical results
To validate the method the biharmonic problems have been solved on two irregular domains presented in Figure 4. As a measure of the quality of the results the mean relative error (MRE) has been assumed. For the boundaryvalue problem the error is obtained by the formula while for the eigenvalue one In Eqs. (33)  are the reference values computed by the differential quadrature method [25], which is a kind of pure numerical, pseudospectral technique. To compute the reference results this technique was also associated with domain transformation.
In Figure 5 two first vibration modes are presented. In order to show the stability of the method due to the shape parameter, in Figure 6 the influence of this parameter on the accuracy as well as on the condition number of the system is shown, when solving BVP. For comparison, similar results obtained by Gaussian kernel are also presented in the figure.

Triangular domain with corner cut-out
In this case the transformation between physical and computational domain is done by the function + 0.75 · (ξ + 1) y = 0.25 · sin (︀ 0.125 · (η + 1) · π )︀ (1 − ξ ) + 0.75 · (ξ + η + ξη + 1) The results obtained are presented in Tables 5 -8 dependently on the type of the boundary conditions and the strategy of their implementation. In Figure 7 two first vibration modes are presented. The influence of the shape parameter on the accuracy and condition number is shown in Figure 8.

Comments to the results
It should be mentioned that all the results contained in the tables have been obtained using Chebyshev-Gauss-Lobatto grid (points concentrated near boundaries) and with the shape parameter ε = 0.01. General analysis of the results obtained indicates that the method provides solutions that are very close to the reference ones.
Comparing MRE between appropriate tables it is easy to notice that the implementation of boundary conditions with the use of strategy 1 leads to more accurate results, despite of the fact that this strategy employs smaller number of interior nodes to discretize the governing equation than strategy 2. The reason may lie in conditioning of the system. Strategy 2 uses more base functions then strategy 1 for the same number of total points. It means that the Hermite polynomials of higher degrees are employed in the discretization process, what makes conditioning of the system worse.
This fact also means that increasing the number of nodes does not always lead to more accurate results, regardless of the strategy used.
Analyzing Figures 6 and 8 one can conclude that the derived bases is stable due to the shape parameter -the influence of this parameter on the conditioning of the system is not as significant as for the conventional Gaussian kernel. Therefore very small values of the parameter still allow to carry out stable computation and give very accurate results, according to the theoretical findings on RBF kernels.

Conclusion
Although the modified kernel method with stable basis is an improvement of the well-known RBF collocation methods it limits their applicability. The RBF kernels easily allow to take advantage of the scattered nodes, while the mentioned method does not. The method can be applied on the Cartesian product grids. It means that the method cannot be used directly to solve the problems with curvilinear boundaries.
In the present paper an extension of the applicability of the kernel method with stable basis has been shown. The method, which can be directly used for problems on regular domains has been applied on irregular ones. By using the blending functions, the transformation between irregular and computational domain has been done. Therefore the computation can be carried out on the regular grid.
Solving higher order equations such as biharmonic ones, the problem with implementation of multiple boundary conditions appears, when using collocation methods. In the paper this problem has been solved by the use of two strategies.
The results obtained indicate that the method provides very accurate results, maintaining its main feature -stability due to the shape parameter. It allows, according to the theoretical findings, to take very small value of this parameter and expect very accurate results.
On the other hand more collocation points require to use basis functions of higher degrees, what worsen the conditioning of the system and the results may not be improved in this case.
Therefore the reasonable choice is to use a coarse grid and control the accuracy by the shape parameter.
It is worth to notice that the use of the method requires to derive, operate and evaluate complicated mathematical expressions, what can viewed as a drawback. But with the use of modern computer algebra systems this task can be significantly simplified.