An automatic methodology for obtaining optimum shape factors for the radial point interpolation method

An automatic methodology for obtaining optimum shape factors for the radial point interpolation method

Péricles L. Machado^I; Rodrigo M.S. de Oliveira^II; Washington C.B. Souza^III; Ramon C.F. Araújo^IV; Maria E.L. Tostes^V; Cláudio Gonçalves^VI

^IUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil, [e-mail]

^IIUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil, [e-mail]

^IIIUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil, [e-mail]

^IVUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil, [e-mail]

^VUniversidade Federal do Pará (UFPA) - Instituto de Tecnologia (ITEC) - Belém, Pará, Brazil, [e-mail]

^VIUniversidade do Estado do Amazonas - Manaus, Amazonas, Brazil, [e-mail]

ABSTRACT

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).

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-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.

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 r_maxis the maximum value assumed by r in Ω. Here, _i is the i^th node of Ω, i = 1, 2, ..., k and P^T () is a polynomial function with M terms. In this work, P^T () 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 (₁), R^T (₂), ..., R^T (_k)]^T and P_o=[P^T (₁), P^T (₂), ..., P^T (_k)]^T , with _i∈ Ω, i = 1 ... k and U_scontains 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

and

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 S_aand S_bare 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 S_aindexed by (i, 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

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

in which = exp (–c (r_{i, j}/ r_max)²) and r_{i, j}= , 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 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 R_o.

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 C_oexists. The core idea of this work is to use C_oas 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 C_ofor Ω, 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(x_i, y_i), given by

is calculated for every

_iin Ω . In (12), _i=(x_i, y_i) represents the i^th node in Ω and K is the associated wavenumber, which is expressed by

In (13), v₀ is the light speed in vacuum. In order to determine C_o, the error

is considered in this work. In (14), C_i(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(C_i(c, ) – C())/C(). Therefore, we can say that C_ois 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 _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 C_oexists 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 C_ocan 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 ƒ_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 λ = v₀/ ƒ_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 Δ_aand 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 C_oproduces 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 C_o, 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 C_ofor the present problem. They were obtained for calculating ∂E_z/∂x (Fig.9a) and ∂E_z/∂y (Fig.9 b) on (9) and (8), respectively. As previously discussed, it is possible to see in Fig.9 that most values of C_ois 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_owere obtained (red points means C_o≈ 100, which is the maximum value assumed by C_oin 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 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_ocan 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_oincreased 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.

VI. ACKNOWLEDGEMENTS

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.

received 15 July 2011; for review 26 July 2011; accepted 12 Dec. 2011

• [1] A. Taflove and S. C. Hagness, Computational Electrodynamics, The Finite-Difference Time-Domain Method, 3rd ed. Artech House, 2005.
• [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. 302307, 1966.
• [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. 16231648, 2002.
• [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.
• [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.
• [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.
• [7] G. R. Liu and Y. T. Gu, An Introduction to Meshfree Methods and Their Programming Springer, 2005.
• [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. 6165.
• [9] G. E. Fasshauer, Meshfree Approximation Methods with MatLab World Scientific Publishing Company, 2007.
• [10] 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.
• [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. 244247.
• [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. 16301639, 1996.
• [13] S. C. Chapra and R. P. Canale, Numerical Methods for Engineers, 6th ed. McGraw-Hill, 2010.
• [14] C. A. Balanis, Advanced Engineering Eletromagnetics John Wiley and Sons, New York, 1989.

Publication Dates

• Publication in this collection
  09 Mar 2012
• Date of issue
  Dec 2011

History

• Accepted
  12 Dec 2011
• Reviewed
  26 June 2011
• Received
  15 June 2007