KRYLOV’S SUBSPACES ITERATIVE METHODS TO EVALUATE ELECTROSTATIC PARAMETERS

Most of the electromagnetic problems can be stated in terms of an inhomogeneous equation Af = g in whi ch A is a differential, integral or integro-differenti al operator, g in the exitation source and f is the unknown function to be determined. Methods of Moments (MoM) is a procedure to solve the equation above and, by  means of an appropriate choice of the Basis/Testing (B/T), the problem can be translated into an equiv alent linear system even of bigger dimensions. In this wo rk we investigate on how the performances of the ma jor Krylov’s subspace iterative solvers are affected by different choice of these sets of functions. More specifically, as a test case, we consider the algeb ric linear system of equations obtained by an elect rostatic problem of evaluation of the capacitance and electr ostatic charge distribution in a cylindrical conduc tor of finite length. Results are compared in terms of ana lytic l/computational complexity and speed of convergence by exploiting three leading iterative m ethods (GMRES, CGS, BibGStab) and B/T functions of Pulse/Pulse (P/P) and Pulse/Delta (P/D) type.


INTRODUCTION
As well known, capacitance C and electrostatic charge distribution q represent important parameters in electrostatic computation above all in applications in which problems of electrostatic charge can determine situations of danger (for example in aeronautics field). Because the computation of C and q leads to the formulation of the problem in terms of integral equations of type Equation (1): where, A is an integro/differential operator, g a known forcing and f to be determined (Sadiku, 2000), it is imperative the exploitation of a procedure of resolution which, through a reduced computational charge (particularly useful in real-time applications), reformulates the problem in the algebrical field in more accessible resolution terms. In such a context, the Method of Moment (MoM) can be considered an excellent candidate for the problem resolution because, by means of a suitable couple of Basis and Testing functions (B/T), it translates the integral equation into a correspondent matrix equation whose solvability requires an iterative numerical procedure above all in those cases in which the matrix dimension is too high to exploit direct procedures. The exploitation of MoM has become popular since in 1967 published a pioneering work about the MoM application in electrostatic problems and the scientific community has produced meaningful works of electrostatic modelling on different geometries exploiting MoM (Chakrabarty et al., 2002;Alad and Chakrabarty, 2012;Prarthan and Chakrabarty, 2001) in which the choice of the B/T couple was asked to the accuracy of the solution, the dimension of the matrix and the evaluation of the elements of the matrix and to problems linked to the conditioning of the same matrix and more, MoM has been applied to many problems in

Science Publications
AJAS several scientific domains. In Electromagnetics MoM has been applied to radiation problems, scattering problems, analysis and design of microstrips; analysis abd design of antennas; lossy structures, propagation problems (Sadiku, 2000;Angiulli and Versaci, 2002;2003;Misran et al., 2008;Elangovan and Perinbam, 2012;Sultan et al., 2012;Wu et al., 2012;Zhang et al., 2013;Mouheddin and Jamel, 2009;Apaydin and Sevgi, 2012;Kumar and Parvathi, 2012). Surely, as Harrington suggests, the choice of B/T functions of Pulse/Delta (P/D) and Pulse/Pulse (P/P) type leads, from one hand, to easy formulations of the elements of the matrix L but, from the other, to slow convergences and reduced levels of accuracy, so, even willing to restore the desired accuracy, it results imperative the increase of the number of subsections, leading sometimes tobad conditionings of the same matrix. In this study, the authors report a comparative study on P/D and P/P functions used in MoM in terms of complexity, convergence and conditioning of the matrix for the evaluation of C and q in a cylindrical conductor of finite length. The numerical data are presented for each B/T couple in terms of conditioning number and convergence. Furthermore, they have been compared three leading iterative methods for the solution of the respective linear systems in Krylov's subspace. Specifically, GMRES with and without restarting procedure, a biorthogonalization algorithm adapted from the bi-Conjugate Gradient Iteration (CGS) and bi-conjugate stabilized method (BicGStab) testing methods on matrixes with number of elements different from zero increasing step by step, evaluating the execution times and the number of necessary iterations. The study is organized as follows. Section 2 is dedicated to a quick overview of the MoM methodology and the analytical formulation and numerical translation of the electrostatic problem under examination. Section 3 describes the philosophy of the resolution of linear systems through Krylov's subspaces and the numerical results, obtained through the translation into MatLab® code, are presented in section 4. Dutiful conclusions and plausible future perspectives close the present work.

CAPACITANCE AND CHARGE DISTRIBUTION IN CYLINDRICAL CONDUCTORS OF FINITE LENGTH
MoM is a procedure which allows to solve numerically non homogeneous equations of type (1). In Particular, expanding f with a series of functions defined in Lthrough Introducing a set of functions weight W = {w 1 , w 2 , ..., w n } in the range of L and supposing fixed an appropriate internal product <g 1 , g 2 > in the data space, the (2) can be reformulated as Equation (3) The possible implementations differentiate themselves for the choice of f n and w n and the order of the summation breaking off (Sadiku, 2000;2010). MoM consists of four basic steps: (a) formulation of an appropriate integral equation; (b) by operation of discretization, transformation of the integral formulation into the equivalent algebraic formulation by means of Basis (B) and Testing (T) functions; (c) computation of the matrix elements; (d) resolution of the corresponding linear system to derive the parameters (Sadiku, 2000). Let us consider a metal cylinder with a circle section of L length and d diameter. Expressing the charge distribution q through functions of expansion, we obtain Equation ( where, f n are known basis functions and α n unknown coefficients. So, the potential at position r due to any charge distribution on the cylinder at position r'can be expressed by Equation (6): with Ω as source region and B(y'-y) as basis functions. Under the hypothesis according to which the axis of the cylinder is place of source points and on the surface they Science Publications AJAS lay points of observation, the problem formulation is one-dimensional, so the basis functions B (due to their localized formulation) restrict the integral to be over the n-th segment. On each y m (computed as m∆y) it is centred a new function, the so called Testing Function T(y-y n ), for which, known the difference of potential V and the number of subsections N, the expression gets the following form Equation (7): T(y)dy with m =1,...,N. Indicating with Equation (8): n y m y mn 2 2 1/ 2 ( n 1) y (m 1) y The matrix of the coefficients (square of N order) and with Equation (9): Known the supposed forcer, the problem can be presented in the matrix form Equation (10): For which the coefficients of expansion α n , obtainable through the resolution of the system (10), give q and C values as follows Equation (11)

P/D Forms as B/T Functions
In this case, the basis functions, defined through rectangles of unitary height and basis localized on each subsection of width y n+1 -y n get the form of (Alad and Chakrabarty, 2012) Equation (12 and 13): So that the formulation of the elements of the matrix A and respective integration get the form Equation (14) with y n+1/2 =y n +∆y/2and y n-1/2 =y n -∆/2 while the forcer becomes Equation (15): and q and C calculated as in (11).

P/P Forms as B/T Functions
Similarly, as for the P/D case, the evaluation of the coefficients of the matrix A and relative integration lead to the expression (Alad and Chakrabarty, 2012) Equation (16) 1/ 2 m m 2 2 1/ 2 2 2 (n 1)/ 2 n 1 2 1/ 2 (n 1)/ 2 (m 1)/ 2 (n 1)/ 2 (m 1)/ 2* (m 1) / 2 (n 1)/ 2 m 1 n 1) 2 2 y y (d / 2) y y y y y (d / 2) y *log (y ) y (d / 2) y y where, y n+1/2 = y n+ ∆y/2 and y n-1/2 = y n -∆y/2, y m+1/2 = y m+ ∆y/2 and y m-1/2 = y m -∆y/2. Finally, the forcer gets the form of Equation (17): In both cases it appears imperative to use an appropriate procedure of resolution of linear systems. Cramer's method cannot, obviously, be suggested because of the prohibitive costs of execution (O(n+1)!) and the direct methods, even they factorize the matrix Ain the product of triangular matrixes, in presence of Science Publications AJAS random matrixes it's not obvious that they are random, too. So, we drive our attention towards iterative methods and, in particular, towards projective methods based on Krylov's subspaces, because they fit well the resolution of linear systems of big dimensions even in presence of scattering of the coefficient matrix.

KRYLOV'S SUBSPACES FOR SOLVING LINEARSYSTEMS
From the linear system (10) and writing the residual r in the form V-Aα, the k-th step is linked to the initial residual r (0) in polynomial basis through (Liu, 2011) Equation (18): With γ as relaxation or acceleration parameter. Setting r (k) = pk (A)r (0) pk(A), results a polynomial of k grade. Defining the vectorial subspace through the linear envelope Equation (19): From (18), it is allowed to write r k ∈K k+1 (A, r (0) ). K m (A; z) is a vectorial space called Krylov's subspace of m order generated by all vectors u∈R n of the form u = p m-1 (A)z with ∂p m-1 ≤m-1. In absence of pre-conditioning, the iterated α (k) gets the form of Equation (20): Belonging to the W k = {z= α (0) +y: y∈K k (A; r (0) )} space. In general, it can be thought to methods building solutions approximated in the form Equation (21): Choosing q k-1 so that α (k) is the best possible approximation in W k.. dim (K m (A; z)) = min (m, ∂z) and, as a consequence, dim (K m (A; z) = f(m)) not decreasing. Increasing m the basis could become numerically unsteady because the vectors tend to be linearly independent so, m fixed, it can be computed an orthonormal basis for K m (A; z) through Arnoldi's procedure which, choosing z 1 = z/||z|| 2 generates an orthonormal basis {z i } through Gram-Schmidt's method and for k = 1,...,m computes With W k = 0 the process ends with a breakdown, or z k+1 = w k /||w|| 2 and k=k+1 and, in step m, z 1 ,…, z m form a base for K m (A; z) so it is solved the linear system computing α k ∈W k as a vector minimizing the Euclidean norm of the residual ||r|| 2 = min z∈ w k ||V-Az|| 2 generating the so-called procedure Generalized Minimum Residual (GMRES) (Liu, 2011).

GMRES
The basis vectors for K m (A; r (0) ) are memorized in columns of Z k with z 1 = r (0) /||r (0) || so the new iterate gets the form α (k) = α (0) +v k z (k) with z (k) minimizing with H matrix of Hessenberg being superior. In terms of computational cost, Arnoldi's procedure is of m 2 /2 order, while GMRES is of m 2 order. Obviously, the procedure can be improved computing the residual at regular intervals and interrupting the algorithm when K gets the sufficient dimension to the tolerance required. Because GMRES requires a considerable computational effort, they have been developed different variables among which also GMRES (m) based on the use of the restarting after m steps. based on the use of the restarting after n steps, it is not certain that the method is not subjected to stagnations and, in particular, we would like not to be necessary to memorize all the vectors of the basis to compute the solution, eliminating also the necessity of two harboured cycles. In such a context, another class of projective methods is based on the algorithm of Lanczos' bi-orthogonalization, in which the main benefits stands in memorizing only a reduced number of vectors (Liu, 2011).

CGS and BicGStab
Lanczos' bi-orthogonalization produces a couple of orthogonal bases for the subspaces Equation ( solving both Aα = V and A T α* = V*. The algorithm is efficient if it is really necessary to solve both the systems and it becomes excessively expensive if it is used only for the resolution of Aα = V. Infact, the matrix A is sometimes available in the matrix-vector form and, in this case, it is not possible to get also A T out in the same form. This way, they're searched some modifications called Transpose-Free Variants which won't use the transposed matrix. A first solution is Conjugate Gradient Squared (CGS) which writes the residual through a Science Publications AJAS squared polynomial φ j (A). Such a choice, in case of irregular convergence, provokes an increase of rounding off errors and possible overflows. This leads to prefer the BicGStab method (Biconjugate Gradient Stabilized) based on the writing of the residual as product of two polynomials Ψ j (A)φ j (A) (Liu, 2011). At any case, it is not possible, in general, to say which of the two methods is more efficient. However, on the basis of the estimations of the computational costs, decreasing the number of the elements different from zero present in the matrix, BicGStab should be quicker. Infact, even in BicGStab there are two matrix-vectors products, these should have a limited cost as to the products carried out m 2 times and to GMRES' high cost of memorization. In this work, the idea is to test methods on matrixes with a number of elements different from zero increasing step by step, considering the time employed and the number of necessary iterations. At any case, dealing with iterative methods, it is necessary to fix the stop conditions, for which, from the recourse relation on the error e (k+1) = Be (k) with B matrix of iteration, we obtain ||e (k+1) ||≤ ||B|| ||e (k) || and applying Schwarz's inequality we find ||e (k+1) ||≤||B||(||e (k+1) -α (k) ). If then ||B||≤1, it stands out the well known rise ||α-α (k+1) particularly useful for the estimation of the number of iterations for the condition satisfying ||e (k+1) ||≤∈ with ε fixed tolerance. To increase the efficiency, the stop test will be carried out on the normalized residual ||r (k)| | ||r (0) ||≤∈ or ||r (K) || ||b||≤∈ while the control on the error occurs through the sequence of rises

NUMERICAL RESULTS
The matrix elements A mn have been computed by means of (14) and (16) by B/T functions of P/D and P/P type respectively. Above the construction of the single matrixes, the correspondent linear systems have been solved through iterative procedures of GMRES, CGS and BicGStab type and, as a comparison, with the Gauss-Jordan matrix inversion. In all cases the forcers used were equal to the excitation voltages obtaining, through the (11), the distributions of electrostatic charge on each single subsection and the values of the capacitance of the charged conducting cylinder. Interested in assigning subdivisions leading to conditions of stability of the solutions respect the disturbance of the entry results, the choice of the number of N subsections has been derived from the comparison between the two criteria. The first, for each value of the length/diameter ratio, assesses the minimums of the trends of the conditioning numbers K(A) = ||A|| ||A -1 ||≥1 (graded index of sensitiveness of the solution to the disturbance of the entry data with values as close to the unit as less sensitive it is the system to the disturbance), N increasing for each A matrix drawn per P/D and P/P functions. In certain cases under examination, the drawn trends, even presenting instability for little values of N,settle themselves step by step N increases of value, to produce, then, strong conditions of instability over a certain N given value. By way of example, Fig. 1 shows the trend of the matrix condition number according to the number of subsections for each couple of B/T functions as to the value 15, while Fig. 2 shows the trend of the capacitance per unit of length according to the number of subsections for each B/T couple as to the value 6 of the length/diameter ratio (no restarting for GMRES). In particular, to the i-th minimum, min i (K(A)), it corresponds an i-th value of N, N i , so that the first criterion gives Nthe minimum value of the N i , min(N i ). The second criterion assesses the trends of the capacitance (pF) when increasing the number of subsections for each B/T couple and for each resolutive procedure of the linear system until the capacitance values are not affected by evident oscillations due to numerical instability. Labelling N j the values for N over which, at any case, it occurs a numerical instability, the second criterion gives N the minimum value of the N j , min (N j ). So, the final assignation for N will be Equation (24) Keeping us, this way in a safety advantage because really far from conditions of instability, but always sufficiently close to the conditions of minimum for K(A). In the case under examination, according to the (24), the final N is tested around 45 units. The comparison for capacitance of specimen (pF), obtained with different B/T functions and for increasing values of length/diameter ratio, have been compared with the analytical solution (Alad and Chakrabarty, 2013) and are reported in Table 1 for each procedure of resolution of the linear system.

AJAS
They have been traced, as marked in Fig. 3, the trends of the capacitance (pF/m) according to the length/diameter ratio for each B/T couple and resolutive procedure of the linear system (no restarting for GMRES), highlighting the instability area, too. The distribution of electrostatic charge (C/m) with value equal to 400 of length/diameter ratio for each B/T couple and relatively to each resolutive procedure of the linear system (no restarting for GMRES) is visualized in Fig. 4. In order to evaluate the convergence of the numerical data on the capacitance, Fig. 5 shows the trend of the capacitance (pF/m) with value equal to 400 of length/diameter ratio, for each B/T couple and as to each resolutive procedure of the linear system (no restarting for GMRES). To assess the convergence speed of the iterative procedures and considering that when they increase the number of the elements different from zero present in the matrix, BicGStab should be quicker than GMRES, the methods have been tested on A matrixes with number of not null elements, increasing step by step and evaluating the number of iterations required and the CPU-Time. Specifically, with reference to the GMRES procedure, when it increases the length/diameter ratio and for different B/T functions, they have been drawn the dimensions of theA matrix which guarantees the convergence and the correspondent number of iterations in absence/presence of restarting (for GMRES). In this last case, the number of iterations has been differentiated in internal/external iterations, correspondent to the cycles nesting. Finally, as to the CGS and BicGStab procedures, always for different B/T functions and for increasing length/diameter ratio, they have been drawn the dimensions of the A matrix guaranteeing the convergence and the correspondent number of iterations and the comparison too of the CPU-time respect the GMRES procedure. The results are presented in Table 2.

CONCLUSION
In the present work, upon the double research about the number of useful subsections and the instability areas got from the comparison between the two criteria of which in the previous paragraph ( Fig. 1 and 2), it has been presented a comparative study about P/D and P/P functions used in MoM in terms of complexity, convergence and conditioning of the matrix for the evaluation of C and q in a cylindrical conductor of finite length. The numerical data are presented for each B/T couple in terms of number of conditioning and convergence.   (Alad and Chakrabarty, 2013), the capacitance values (pF/m) obtained with GMRES (no restarting), CGS, BicGStab and with Gausse-Jordan's inversion when varying the length/diameter ratio and for P/P couples. Variation of capacitances of hollow cylinder with height by diameter ratio using P/D and P/P methods and analytical one have been carried out by means Gauss-Jordan, GMRES. CGS and BicGStab algorithms. In particular, as for all the procedures of resolution considered, P/P with BicGStab performs better the analytical results (Fig. 3). The charge distribution q (Fig.  4) presents maximums at the extremity of the specimen and, through marked inclinations, it is ranged the minimum in correspondence at the middle of the specimen. The convergence data show that BicGStab procedure with P/P couple that converges faster than other procedures (Fig. 5) presenting a better adherence to the analytical solution. The comparison between GMRES (with and without restarting), CGS and BicGStab in terms of matrix dimensions guaranteeing the convergence, the number of iterations (differentiated into internal and external) and CPU-time according to the length/diameter ratio and for different B/T couples, presented in Table 2, has highlighted that P/P BicGStab, even presenting an increase of the number of iterations, provides a decreasing trend of CPU-time. Obviously, what has been dealing with in the present work is desirable to be tested on specimen with more complex geometries and on more matrixes comparing BicGStab with GMRES with different restarting options of adaptive type. However, the choice of operating on a cylindrical specimen was dictated by the fact that more complex specimens can be treated modularly as a set of cylindrical elements. Obviously, the present study has the limitation of treating only two pairs of functions B/T. In the future, it would required a more detailed study of comparison between functions belonging to a set of pairs B/T characterizing by a bigger cardinality. Because of the bad conditioning of the matrix of the system, it is imperative the exploitation of suitable pre-conditioners.