On-line version ISSN 2179-1074
J. Microw. Optoelectron. Electromagn. Appl. vol.10 no.2 São Caetano do Sul Dec. 2011
An automatic methodology for obtaining optimum shape factors for the radial point interpolation method
Péricles L. MachadoI; Rodrigo M.S. de OliveiraII; Washington C.B. SouzaIII; Ramon C.F. AraújoIV; Maria E.L. TostesV; Cláudio GonçalvesVI
IUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil, firstname.lastname@example.org
IIUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil, email@example.com
IIIUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil, firstname.lastname@example.org
IVUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil, email@example.com
VUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil, firstname.lastname@example.org
VIUniversidade do Estado do Amazonas - Manaus, Amazonas, Brazil, email@example.com
In this letter, a methodology is proposed for automatically (and locally) obtaining the shape factor c for the Gaussian basis functions, for each support domain, in order to increase numerical precision and mainly to avoid matrix inversion impossibilities. The concept of calibration function is introduced, which is used for obtaining c. The methodology developed was applied for a 2-D numerical experiment, which results are compared to analytical solution. This comparison revels that the results associated to the developed methodology are very close to the analytical solution for the entire bandwidth of the excitation pulse. The proposed methodology is called in this work Local Shape Factor Calibration Method (LSFCM).
Index Terms: improved numerical precision, matrix inversion difficulties, optimum shape factor calculation, radial point interpolation method (RPIM).
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 , . Meshless methods, such as the Radial Point Interpolation Method (RPIM), have become an important alternative to solve numerically problems involving partial differential equations , , , , due to the fact it provides greater geometric flexibility ,  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  (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 . 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 , , 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 . This issue has been treated in literature by using methods such as the leave-one-out cross-validation algorithm (LOOCV) , , , 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() in space. This function can be interpolated in a support domain Ω centered at by
in which is the Gaussian radial basis function,, c > 0 is the shape factor and rmax is the maximum value assumed by r in Ω. Here, i is the ith node of Ω, i = 1, 2, ..., k and PT () is a polynomial function with M terms. In this work, PT () is given by [1, x, y] and M = 3.
When (1) is considered for all nodes in Ω, one obtains the matrix equation
where, Ro =[RT (1), RT (2), ..., RT (k)]T and Po =[PT (1), PT (2), ..., PT (k)]T , with i ∈ Ω, i = 1 ... k and Us contains all the values assumed by u(i) in Ω.
In order to ensure that a unique solution is obtained for a and b in (2), the condition a = 0 is imposed .
With some algebra, it can be shown that
This way, (1) can also be written as
in which Ψ() 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 , .
As far as Sa and Sb are constant matrices (because the nodes' coordinates are fixed), the partial derivative of Ψl() with respect to v is given by
where l = 1 ...k, , is the element of matrix Sa indexed by (i, l), is the element of Sb 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 . The associated spatial derivatives are approximated by (7) and the field updating equations
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  was implemented for truncating the analysis domain.
III. THE LOCAL SHAPE FACTOR CALIBRATION METHOD (LSFCM)
Previous works consider c ≈ 0 as a global parameter ,  (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 . 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 Ro in expanded form is given by 
in which = exp (c (ri, j / rmax)2) and ri, j = , it is observable that for the case of Fig. 2, the distances r1,n and r2,n (from nodes 1 and 2 to node n, respectively), with n > 2, are approximately equal. This way, it is evident that . Additionally, when exp(0) = 1. Observing that , 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 Ro.
In these cases, the associated Gaussian basis functions should present higher exponential decays as a function of . On the other hand, the exponential decays (which are associated to c) should be spatially compatible to the smaller wavelengths present in u() for obtaining precise interpolations. Even though one can adjust the shape factors by the trial and error method , an automatic procedure is required.
This procedure can be obtained if Fig.3 is carefully observed. It contains plots of the percentage interpolation error (%) for a given signal u() 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 Co exists. The core idea of this work is to use Co as the shape factor for avoiding matrix inversion impossibilities and, in addition, to improve the interpolation precision. In practice, as long as u() 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() and its maximum significative frequency max, the frequency max can be used as a reference parameter for determining Co for Ω, in such way the lower frequency components of u() can also be properly interpolated in Ω.
This hypothesys can be verified if, for given a support domain Ω, a calibration function C(xi, yi), given by
is calculated for every i in Ω . In (12), i =(xi, yi) represents the ith node in Ω and K is the associated wavenumber, which is expressed by
In (13), v0 is the light speed in vacuum. In order to determine Co, the error
is considered in this work. In (14), Ci(c, ), which is the interpolated version of C(), is obtained by using (1) and (12). From (14), the percentage error can be calculated by %(c) = 100(Ci(c, ) C())/C(). Therefore, we can say that Co is a root of the percentage error, when the calibration function is considered, in such way that Co > 0.
In this work, the modified regula falsi method  is used for determining Co 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 Co which satisfies (15) with max = 104. The regula falsi algorithm used in this work is illustrated by Fig.4.
In a few cases, it is not possible to determine if Co exists in the range 1 < c < 50 because (50). (1) > 0 (in some cases it does not). For these cases, a minimization algorithm is performed for the function |(c)| in the referred interval. If Co 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 Co.
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 . 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 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.
For the performed experiments, initially global values of c were used for testing purposes (c = 0.1, c = 7.4 and c = 8.5) with k = 12. The average spacing among points is Δa = (Fig.5b), where λ = v0/ max. Then, c was locally calculated (specifically for each support domain) by applying the methodology presented in this paper, and a new simulation was executed. For this case, the parameters Δa and k were kept unchanged (Δa = and k = 12). The precision and stability criteria of the RPIM algorithm follows .
The problem was also solved analytically by using the solution presented by , and additional numerical data were generated by using the FDTD method. The Fourier transform was applyed to the transient signals in order to make the comparisons to the analytical solution feasible.
Fig.7 shows a graphical comparison among the analytical and numerical solutions for electric field at ℓx = 20 mm, with Δa = , for local and global shape parameters. For FDTD, due to the staircase effect, it was necessary to discretize space with Δ = in order to get results closer to that generated with RPIM (Δa = ). Fig.8 shows similar results for ℓx = 38 mm.
In Figs.7 and 8, it is possible to see that the use of local values Co produces the closest curves to the analytical solution for the entire band of frequencies. When the RPIM method employs c = 7.4, for example, it is possible to see errors % of 15.62% for the higher frequency components. However, with the local shape factors Co, the RPIM algorithm produces a maximum error of 1.24% for Δa = (see Tables I and II for numerical data regarding ℓx = 20 mm).
For the space discretization illustrated by Fig.5b, it is possible to use c = 0.1 for the entire domain, with no difficulties for inverting the matrices. It is possible to see from Figs.7 and 8 that the maximum error produced by the RPIM algorithm in this situation is close to 2%. However, it is worth to emphasize that it is not always possible to use such small global value for c, specially if highly heterogeneous distributions of points are needed to represent the analysis region. For the arrangement of points of Fig.5b, it was not possible to set c = 0.01.
Fig.5a defines 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).
Finally, Fig.9 shows the spatial distributions of Co for the present problem. They were obtained for calculating ∂Ez/∂x (Fig.9a) and ∂Ez/∂y (Fig.9 b) on (9) and (8), respectively. As previously discussed, it is possible to see in Fig.9 that most values of Co 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 Co were obtained (red points means Co ≈ 100, which is the maximum value assumed by Co in this example).
V. FINAL REMARKS
The results of the performed experiment confirm the hypothesis in this work: a calibration function can be used for calculating Co, 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 Co 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 Co 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.
The authors are thankful to UFPA, to UEA and FAPEAM for infrastructure and financial support. We also acknowledge reviewers for all the valuable recommendations for improving the quality this paper.
 A. Taflove and S. C. Hagness, Computational Electrodynamics, The Finite-Difference Time-Domain Method, 3rd ed. Artech House, 2005. [ Links ]
 K. Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Trans. Antennas and Propagation, vol. 14, pp. 302307, 1966. [ Links ]
 J. G. Wang and G. R. Liu, "A point interpolation meshless method based on radial basis functions," Int. J. Numer. Method, vol. 54, pp. 16231648, 2002. [ Links ]
 Y. Yu and Z. Chen, "A 3-d radial point interpolation method for meshless time-domain modeling," Microwave Theory and Techniques, IEEE Transactions on, vol. 57, no. 8, pp. 2015 2020, aug. 2009. [ Links ]
 X. Chen, Z. Chen, Y. Yu, and D. Su, "An unconditionally stable radial point interpolation meshless method with laguerre polynomials," Antennas and Propagation, IEEE Transactions on, vol. 59, no. 10, pp. 3756 3763, oct. 2011. [ Links ]
 T. Kaufmann, C. Engstrom, C. Fumeaux, and R. Vahldieck, "Eigenvalue analysis and longtime stability of resonant structures for the meshless radial point interpolation method in time domain," Microwave Theory and Techniques, IEEE Transactions on, vol. 58, no. 12, pp. 3399 3408, dec. 2010. [ Links ]
 G. R. Liu and Y. T. Gu, An Introduction to Meshfree Methods and Their Programming. Springer, 2005. [ Links ]
 T. Kaufmann, C. Fumeaux, and R. Vahldieck, "The meshless radial point interpolation method for time-domain electromagnetics," in IEEE MTT-S International Microwave Symposium, 2008, pp. 6165. [ Links ]
 G. E. Fasshauer, Meshfree Approximation Methods with MatLab. World Scientific Publishing Company, 2007. [ Links ]
 G. E. Fasshauer and J. G. Zhang, "On choosing "optimal" shape parameters for rbf approximation," Numerical Algorithms, vol. 45, no. 1-4, pp. 345368, 2007. [ Links ]
 T. Kaufmann, C. Engstrom, and C. Fumeaux, "Residual-based adaptive refinement for meshless eigenvalue solvers," in Electromagnetics in Advanced Applications (ICEAA), 2010 International Conference on, 2010, pp. 244247. [ Links ]
 S. D. Gedney, "An anisotropic perfectly matched layer absorbing media for the truncation of fdtd latices," IEEE Trans. Antennas and Propagation, vol. 44, pp. 16301639, 1996. [ Links ]
 S. C. Chapra and R. P. Canale, Numerical Methods for Engineers, 6th ed. McGraw-Hill, 2010. [ Links ]
 C. A. Balanis, Advanced Engineering Eletromagnetics. John Wiley and Sons, New York, 1989. [ Links ]
received 15 July 2011; for review 26 July 2011; accepted 12 Dec. 2011