An Automatic Methodology for Obtaining Optimum Shape Factors for the Radial Point Interpolation Method

Péricles L. Machado, Rodrigo M.S. de Oliveira, Washington C .B. Souza, Ramon C.F. Araújo, Maria E.L. Tostes Universidade Federal do Pará (UFPA) − Instituto de Tecnologia (ITEC) − Belém, Pará, Brazil pericles.raskolnikoff@gmail.com, rmso@ufpa.br, csakso usa@yahoo.com.br, ramon.araujo@itec.ufpa.br and tostes @ufpa.br Cláudio Gonçalves Universidade do Estado do Amazonas − Manaus, Amazonas, Brazil gonca.claudio@gmail.com


I. INTRODUCTION
One of the most used numerical methods for solving Maxwell's equations in time domain is the finite-difference (FD) technique, on which the finite-difference time-domain method (FDTD) is based [1], [2].Meshless methods, such as the Radial Point Interpolation Method (RPIM), have become an important alternative to solve numerically problems involving partial differential equations [3], [4], [5], [6], due to the fact it provides greater geometric flexibility [7], [8] than FD-based methods.This kind of methdology employs a set of points for representing the analysis region, instead of grids.The field components are locally interpolated by using subgroups of points, called support domains [3] (Fig. 1).As it is well known, Gaussian, multiquadrics, logarithmic and splines functions are used in RPIM method as basis functions that depend on an arbitrary parameters [9].In particular, Gaussian basis depends on a free parameter c, known as the shape factor, which affects the accuracy of the interpolation of functions of interest in a particular support domain.In previous works [3], [8], the parameter c is defined globally usually as c = 0.01.However, this value is not adequate for every support domain, due to loss of accuracy and especially due to matrix inversion impossibilities [9].This issue has been treated in literature by using methods such as the leave-one-outcross-validation algorithm (LOOCV) [9], [10], [11], which calculates a optimum global shape parameter by using statistical analysis.
In order to improve the accuracy of the RPIM method, this paper presents a formulation for computing c in an automatic way, locally for each support domain (Fig. 1), in such way reduced interpolation errors are obtained and matrix inversion difficulties are avoided.This is accomplished by using a high frequency signal, called here calibration function.This way, the proposed method is named Local Shape Factor Calibration Method (LSFCM).
Excellent agreements to analytical solution were observed.

II. REVIEW OF THE RPIM METHOD
Consider a function u(x) in space.This function can be interpolated in a support domain Ω centered at x by in which r i (x) = e −c(r/r max ) 2 is the Gaussian radial basis function, r = (x − x i ) 2 + (y − y i ) 2 , c > 0 is the shape factor and r max is the maximum value assumed by r in Ω.Here, x i is the i th node of Ω, i = 1, 2, ..., k and P T (x) is a polynomial function with M terms.In this work, P T (x) is given by [1, x, y] and M = 3.
When (1) is considered for all nodes in Ω, one obtains the matrix equation where, R o = [R T (x 1 ), R T (x 2 ), ..., R T (x k )] T and P o = [P T (x 1 ), P T (x 2 ), ..., P T (x k )] T , with x i ∈ Ω, i = 1...k and U s contains all the values assumed by u(x i ) in Ω.
In order to ensure that a unique solution is obtained for a and b in (2), the condition P T o a = 0 is imposed [3].
With some algebra, it can be shown that and This way, (1) can also be written as in which Ψ(x) is a vector containing samples of shape functions associated to each node i in Ω.It is important to mention that the shape functions in Ω must satisfy the Kronecker delta property [3], [9].
As far as S a and S b are constant matrices (because the nodes' coordinates are fixed), the partial derivative of where l = 1...k, ∂Ψ ∂v = [ ∂Ψ 1 ∂v , ∂Ψ 2 ∂v , ..., ∂Ψ k ∂v ], S a i,l is the element of matrix S a indexed by (i, l), S b j,l is the element of S b indexed by ( j, l) and v = x or v = y.Finally, the partial derivative of u with respect to v can be expressed by In this work, Maxwell's equations are solved in 2-D space, by considering the TMz mode [1].The associated spatial derivatives are approximated by (7) and the field updating equations x, j ∂yΨ j (10) are obtained.In ( 8)-( 10), central finite-differences are used to approximate the time derivatives.
Here, the UPML (Uniaxial-Perfectly Matched Layer) formulation developed by Gedney [12] was implemented for truncating the analysis domain.

III. THE LOCAL SHAPE FACTOR CALIBRATION METHOD (LSFCM)
Previous works consider c ≈ 0 as a global parameter [3], [8] (c = 0.01 is often used).Although small values of c can produce highly accurate results, it is not always possible to compute ( 8)-( 10) due to matrix inversion difficulties [9].This occurs because when a node is placed close enough to another in Ω, as illustrated by Fig. 2, and c is close to zero, the Gaussian basis functions tend to be constant (tend to the unity) between the mentioned nodes and thus the associated matrices tend to be non-invertible.Mathematically, because R o in expanded form is given by [9] in which R i, j o = exp −c(r i, j /r max ) 2 and r i, j = (x i − x j ) 2 + (y i − y j ) 2 , it is observable that for the case of Fig. 2, the distances r 1,n and r 2,n (from nodes 1 and 2 to node n, respectively), with n > 2, are approximately equal.This it is easy to see that the situation depicted by Fig. 2 makes the lines one and two of (11) almost linearly dependent, promoting inversion difficulties for R o .
In these cases, the associated Gaussian basis functions should present higher exponential decays as a function of x.On the other hand, the exponential decays (which are associated to c) should be spatially compatible to the smaller wavelengths present in u( x) for obtaining precise interpolations.Even though one can adjust the shape factors by the trial and error method [9], an automatic procedure is required.This procedure can be obtained if Fig. 3 is carefully observed.It contains plots of the percentage interpolation error (E % ) for a given signal u(x) as a function of c for three different spatial arrangement of nodes: symmetrical (even), slightly irregular and irregular arrangements.As one can observe, as c goes to zero, the interpolation tends to produce very low percentage errors for every case.However, the curve is discontinuous around c = 0 due to the matrix inversion impossibilities in this range.Additionally, it is possible to observe that, for each case, a second root C o exists.The core idea of this work is to use C o as the shape factor for avoiding matrix inversion impossibilities and, in addition, to improve the interpolation precision.In practice, as long as u(x) is not known analytically, the error function E cannot be calculated a priori.
Based on the above discussion, an hypothesis has been formulated: given a function u(x) and its maximum significative frequency f max , the frequency f max can be used as a reference parameter for determining C o for Ω, in such way the lower frequency components of u(x) can also be properly interpolated in Ω.This hypothesys can be verified if, for given a support domain Ω, a calibration function C(x i , y i ), given by is calculated for every x i in Ω.In (12), x i = (x i , y i ) represents the i th node in Ω and K is the associated wavenumber, which is expressed by In ( 13), v 0 is the light speed in vacuum.In order to determine C o , the error is considered in this work.In ( 14), C i (c, x), which is the interpolated version of C(x), is obtained by using ( 1) and (12).From ( 14), the percentage error can be calculated by E % (c) = 100(C i (c, x) −C(x))/C(x).Therefore, we can say that C o is a root of the percentage error, when the calibration function is considered, in such way that C o > 0.
In this work, the modified regula falsi method [13] is used for determining C o from ( 14), in such way that Here, the considered searching range for c is 1 ≤ c ≤ 50.Typically, three to six regula falsi's steps are necessary for obtaining C o which satisfies (15) with E max = 10 −4 .The regula falsi algorithm used in this work is illustrated by Fig. 4.
is the function which root is to be estimated and end if 22: end while 23: x is the final approximation to the root of f .In a few cases, it is not possible to determine if C o exists in the range 1 ≤ c ≤ 50 because E (50).E (1) > 0 (in some cases it does not).For these cases, a minimization algorithm is performed for the function |E (c)| in the referred interval.If C o can not be found in the initial searching range, a new interval is defined for the investigation (e.g.50 ≤ c ≤ 100).Finally, it is of fundamental importance to observe that the spatial derivatives in ( 8)-( 10) are considered separately for determining the values assumed by C o .

IV. NUMERICAL EXPERIMENTS AND DISCUSSION
The RPIM 2-D meshless method was applied for simulating the electromagnetic scattering of a plane wave by a metallic cylinder immersed in free space, such as illustrated by Fig. 5. Here, the TMz mode is employed [8].The parameters of the experiment are: radius of the metallic cylinder a = 100 mm; the electromagnetic wave propagates in vacuum and it is excited by a wideband monocycle pulse (Fig. 6a) with maximum significant frequency f max = 1.5 GHz (Fig. 6b); the dimensions of the analysis region are 3m × 7m.The discrete analysis region is partially shown by Fig. 5b.Finally, Fig. 9 shows the spatial distributions of C o for the present problem.They were obtained for calculating ∂E z /∂x (Fig. 9a) and ∂E z /∂y (Fig. 9b) on ( 9) and ( 8), respectively.As previously discussed, it is possible to see in

V. FINAL REMARKS
The results of the performed experiment confirm the hypothesis in this work: a calibration function can be used for calculating C o , since it contains the higher frequency component of the signal to be propagated.The developed methodology was computationally implemented and it was confirmed numerically that the procedure is suitable for automatically obtaining the shape factors of the Gaussian functions locally (for each support domain).Besides it makes the RPIM method more autonomic and more accurate for applications that involves multi-scale techniques, the new methodology prevents difficulties regarding matrices inversions associated to the use of small values of c.
Although in this paper the developed methodology (LSFCM) has been developed numerically, it could be improved if analytical calculation of C o can be performed.This would suppress the necessity of using root-finding algorithms, such as the regula falsi method used in this paper (the calculation of C o increased the processing time in approximately 25% for the numerical examples in this work).
It should by observed that in many applications, such as those involving the heat equation, the maximum significant frequency is not known straightforwardly.This way, further investigation is necessary in this direction.
It is worth to mention that the proposed methodology can be easily extended to 3D problems.

Figure 2 .
Figure 2. Example of support domain with two nodes positioned close to each other.

Figure 3 .
Figure 3. E % versus c for even, slightly irregular and irregular arrangements of k = 12 interpolating nodes for u( x) = C( x) and graphical definition of C o .

Figure 5 .Figure 6 .
Figure 5. (a) Geometric configuration of the problem and points used for calculating E z and (b) part of the set of points used for representing the analysis region ( E and H are not calculated at the same points in space [8]).

Fig.5a definesFigure 7 .
Fig.5adefines the points where the electric field was calculated in this work, where ℓ x is a distance measured from the right edge of the cylinder surface (parallelly to x).

Figure 9 .
Figure 9. Spatial distribution of log 10 (C o ) for calculating (a) ∂E z /∂x and (b) ∂E z /∂y for the present problem.

Fig. 9
Fig.9 that most values of C o is in the initial searching range for c (from 1 to 50).However, in a few cases (mostly near the cylinder's border), where points are placed closer to each other, higher values of C o were obtained (red points means C o ≈ 100, which is the maximum value assumed by C o in this example).

Table I ANALYTICAL
AND NUMERICAL RESULTS (ℓ x = 20 MILIMETERS)