SciELO - Scientific Electronic Library Online

vol.10 issue2A new multilevel smoothing method for wavelet-based algebraic multigrid poisson problem solver author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand




Related links


Journal of Microwaves, Optoelectronics and Electromagnetic Applications

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,
IIUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil,
IIIUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil,
IVUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil,
VUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil,
VIUniversidade do Estado do Amazonas - Manaus, Amazonas, Brazil,




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 [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-out cross-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.



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 [3].

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 [3], [9].

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 [1]. 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 [12] was implemented for truncating the analysis domain.



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 Ro in expanded form is given by [9]



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 [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 (%) 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 [13] 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 = 10–4. 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.



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 ƒ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 [8].

The problem was also solved analytically by using the solution presented by [14], 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).



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.



[1] A. Taflove and S. C. Hagness, Computational Electrodynamics, The Finite-Difference Time-Domain Method, 3rd ed. Artech House, 2005.         [ Links ]

[2] K. Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Trans. Antennas and Propagation, vol. 14, pp. 302–307, 1966.         [ Links ]

[3] J. G. Wang and G. R. Liu, "A point interpolation meshless method based on radial basis functions," Int. J. Numer. Method, vol. 54, pp. 1623–1648, 2002.         [ Links ]

[4] 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 ]

[5] 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 ]

[6] 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 ]

[7] G. R. Liu and Y. T. Gu, An Introduction to Meshfree Methods and Their Programming. Springer, 2005.         [ Links ]

[8] 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. 61–65.         [ Links ]

[9] G. E. Fasshauer, Meshfree Approximation Methods with MatLab. World Scientific Publishing Company, 2007.         [ Links ]

[10] G. E. Fasshauer and J. G. Zhang, "On choosing "optimal" shape parameters for rbf approximation," Numerical Algorithms, vol. 45, no. 1-4, pp. 345–368, 2007.         [ Links ]

[11] 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. 244–247.         [ Links ]

[12] S. D. Gedney, "An anisotropic perfectly matched layer absorbing media for the truncation of fdtd latices," IEEE Trans. Antennas and Propagation, vol. 44, pp. 1630–1639, 1996.         [ Links ]

[13] S. C. Chapra and R. P. Canale, Numerical Methods for Engineers, 6th ed. McGraw-Hill, 2010.         [ Links ]

[14] 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

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License