A Comparative Analysis of Numerical Inverse Laplace Transform Methods for Electromagnetic Transient Analysis

Nowadays, transient analysis of electromagnetic systems is widely adopted as a virtual prototyping tool in the design of electrical and electronic devices. Time domain (TD) methods for the solutions of Maxwell’s equations represent the most obvious approach to tackling this type of problem. Frequency domain (FD) analysis techniques are also commonly used to return the TD response through the inverse Fourier transform. This approach suffers from the need to compute the frequency response over a wide frequency range that has to be finely sampled to avoid aberrations in the inverse transform. A valuable alternative is represented by the numerical inversion of the Laplace transform (NILT) that conjugates the capability to return the TD responses with the possibility to adopt complex frequency models that are free of most of the issues of TD techniques. The aim of this work is to perform a comparative analysis of two of the most popular methods for the computation of the numerical inverse Laplace transform, known as the NILT, and fast inversion of the Laplace transform (FILT), clarifying the issues and limitations behind these methods and pointing out possible solutions.


A Comparative Analysis of Numerical Inverse Laplace Transform Methods for Electromagnetic Transient Analysis
Fabrizio Loreto , Daniele Romano , Giuseppe Pettanice , and Giulio Antonini , Senior Member, IEEE Abstract-Nowadays, transient analysis of electromagnetic systems is widely adopted as a virtual prototyping tool in the design of electrical and electronic devices.Time domain (TD) methods for the solutions of Maxwell's equations represent the most obvious approach to tackling this type of problem.Frequency domain (FD) analysis techniques are also commonly used to return the TD response through the inverse Fourier transform.This approach suffers from the need to compute the frequency response over a wide frequency range that has to be finely sampled to avoid aberrations in the inverse transform.A valuable alternative is represented by the numerical inversion of the Laplace transform (NILT) that conjugates the capability to return the TD responses with the possibility to adopt complex frequency models that are free of most of the issues of TD techniques.The aim of this work is to perform a comparative analysis of two of the most popular methods for the computation of the numerical inverse Laplace transform, known as the NILT, and fast inversion of the Laplace transform (FILT), clarifying the issues and limitations behind these methods and pointing out possible solutions.

I. INTRODUCTION
T RANSIENT analysis of electromagnetic problems is an important task in the design of electrical and electronic systems.Indeed, it allows clarifying physical phenomena such as losses and dispersion that significantly affect system performances on a short time scale [1].Time domain (TD) solutions of Maxwell's equations have been widely used in the modeling of transmission lines and interconnects [2], [3], [4], [5], [6], antennas [7], grounding systems [8], and printed circuit boards (PCBs) performance analysis [9].
Several methods have been proposed to solve Maxwell's equations in the TD, in differential and integral forms.Among them, the most popular ones are the finite difference TD (FDTD) method [10], [11], the TD finite element method (TDFEM) [12], the TD integral equations method (TDIEM) [13], the partial element equivalent circuit (PEEC) method [14].TD methods have some indisputable advantages over the corresponding frequency domain (FD) counterparts.The most significant one relies on the fact that they allow characterizing an electrical system in a broadband frequency range with just one run, without the need to repeat the computation for different frequency samples, as required by FD techniques.Then, they also allow easy modeling of nonlinear materials or the incorporation of nonlinear devices, while this is not possible, or at least it is not so straightforward when using FD methods that assume the superposition principle is applicable [15].On the other hand, FD techniques are easier to be implemented, they do not have any stability issues, and they allow easy modeling of dispersive materials.
An alternative to both pure TD and FD methods is represented by the Laplace transform (LT) method.While the forward LT can be carried out easily and Maxwell's equations can be transformed in a straightforward way, the inverse LT is more challenging and a closed form is possible only in relatively simple cases.In general, it has to be computed numerically by properly choosing the Bromwich contour and exploiting Cauchy's theorem.This work aims to compare the two methods for the numerical computation of the inverse LT which have found application in electromagnetic modeling, namely the numerical inversion of the LT method (NILT) [16], and Fast Inversion of the LT (FILT) method [17], addressing their pros and cons.
The NILT method has been proposed for the first time in the 70s' [16] and has been applied mostly to high-speed interconnect modeling [18].It has received a renewed interest more recently thanks to some improvements [19], [20] that have been proposed to overcome its limitations.Furthermore, it has been adopted in the framework of the PEEC method [14] thus making it possible to use NILT to solve general 3-D EM problems [21], [22], [23].
The FILT method has been introduced in [17] relying on an expansion of the kernel of the inverse LT different from that used in NILT.It has been applied to the transient analysis of frequency-dependent media [24] and plasmonic systems in the framework of the FDTD method [25], [26], [27] and the boundary element method [28].Recently, it has been employed for the TD evaluation of the specific energy loss for pulse incidence on lossy media [29].
The common feature of NILT and FILT techniques is that they both provide the time-domain response at a time sample independently of another one; therefore, they are well suited for performing parallel computations, resulting in a significant reduction of processing time.Furthermore, the time samples can be arbitrarily chosen, thus making it possible to select them adaptively.Finally, the TD response obtained through these techniques is not affected by the accumulation error that occurs using time-stepping methods.
This work is organized as follows: Section II summarizes the NILT and FILT approaches, pointing out their distinctive features, pros, and cons.Section III illustrates how the NILT can be used along with the PEEC method and how underflow issues may arise for electrically large problems.Section IV presents the use of NILT and FILT approaches for the transient analysis of multiconductor transmission lines (MTLs).The observations in Sections III and IV are corroborated by the numerical tests presented in Section V for two pertinent case studies.The conclusions and future works are outlined in Section VI.

II. COMPLEX FD TECHNIQUES OF INVERSE LT
The LT is a powerful analytical tool for the solution of linear systems of differential equations, allowing the elimination of one or more variables and moving in a transformed domain of a complex variable s.The transformed system is usually easier to handle.In EM modeling and circuit theory often the suppressed variable is the time t.This permits avoiding simple and partial derivatives that are a constituent part of the TD model.
Applying the LT methods to linear time-invariant EM systems results in an algebraic system of equations being A(s) the characteristic matrix of the system, X(s) the unknowns vector and U(s) the vector containing the sources acting in the system.The TD solution x(t) of (1) requires computing the inverse LT of X(s) [30] that in general reads where the complex exponential e st is the kernel of the inverse transformation.
The analytical solution is possible only for a limited set of cases including RLC time-invariant circuits [15], [31] of reduced size.In general, the integral in (2) has to be performed in the complex plane and this requires considerable analytical efforts that could lead to closed-form solutions only for a restricted class of functions.Fortunately, many approximation algorithms widely used in engineering modeling applications have been proposed in the last century.Among them, we remind the Weeks [32] method, the Dubner and Abate [33] method, the Stehfest [34] method (widely used in groundwater flow and petroleum reservoir applications) and the Durbin [35] method.A detailed comparison and description of all the approximate Laplace inverse transform methods were provided by Davies and Martin [36].As pointed out in [37], it is difficult to recommend just one inversion method, since the performance is strictly dependent on the function type and, hence, on the application.
In the field of TD EM modeling and circuit simulation, there are two methods that have shown the best compromise between computational efficiency and accuracy.
1) The NILT method was proposed by Nakhla et al. [16] and is based on the Padé approximation of the kernel in (2); 2) The method proposed by Hosono [17], more recently referred to as the Fast Inverse LT (FILT) method, is based on the approximation of the kernel in (2) by a trigonometric hyperbolic function.Although these two methods are suitable for both EM and circuit modeling areas, the former has become more popular in TD circuit analysis and the latter in the transient EM modeling of antennas and microwave devices.
Moreover, the NILT/FILT techniques represent a valid alternative to the well-established inverse fast Fourier transform (IFFT) method.Indeed, in order to restore aberrations-free TD results starting from an FD (s = jω) discrete representation of the system (1), it is required to compute the system response solution over a pertinent frequency range using an appropriate number of frequency samples on the imaginary axis of the complex plane.These requirements often lead to significant computational efforts.On the contrary, NILT/FILT techniques return the TD response by taking advantage of Cauchy's theorem that requires computing the response only on the singularities of the kernel.This implies that a reduced number of computations need to be performed in the complex plane.

A. TD Simulation via NILT
As mentioned before, the NILT technique is based on the approximation of the exponential kernel in (2) through its Padé approximant [38], leading to being where P N (z) and Q M (z) are polynomials of order N and M, respectively; z i and K i denote the poles and residues of the rational function P N (z)/Q M (z).The choice M = N + 2 has been shown to guarantee the solution stability [19].In (3) the change of variables z = st is assumed.The evaluation of the integral in (3) is performed through Cauchy's residues theorem applied to a path in the right-hand half-plane.It leads to where M is assumed to be an even integer.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.It is clear from (5) that, for each time sample t, it is required a number of M/2 evaluations of X(s) on the complex plane, over known points defined by the Padé poles.An illustrative pole trajectory on the complex plane is depicted in Fig. 1, considering a time interval [1], [2], [3], [4], [5] ns and assuming M = 12.
Generally speaking, all the techniques to compute the inverse LT are affected by a progressive loss of accuracy occurring with the time increasing.This is also the case with the NILT and FILT methods.The NILT approximation error [19] and given by the following equation: where The Padé approximation of the exponential function in terms of the rational function of order N /M matches the first (N + M + 1) terms of its Taylor-series expansion.As a result, higher values of M and N lead to better accuracy in TD.Unfortunately, it is not feasible to employ arbitrarily large values of M without that rounding errors in the residues K i affect significantly the solution.This issue is easily explained through a simple example regarding the NILTbased computation of the function f (t) = sin(ωt), with ω = 6.28 • 10 9 rad/s.In Fig. 2 it is shown that employing M = 24 for the inverse transform guarantees a good accuracy until to 6 ns but is not enough for larger times.
In order to preserve the accuracy over the entire time window, an option is to increase further the number of terms in the summation (5) choosing M = 28, but round errors in the high Padé residues blow up the solution (in this case the NILT results are scaled by a 10 −11 factor).Hence, for this example, the choice is limited to values M < 28.Obviously, this is valid for general inverse LTs: every inversion has its expansion limit when approached through the Padé method.

B. TD Simulation via FILT
Following the procedure introduced by Hosono [17], the inversion integral in ( 2) is approximated as: by the introduction of the approximate inversion kernel being α a parameter that has to be properly chosen.Introducing ( 9) in ( 8), it is possible to obtain as follows: where and K is the truncation number of the series.The FILT approximation error ϵ F = |x(t) − x ap (t)| is given by the following equation: being L a constant and E a finite vector [39].In conclusion, the error is of exponential order with α.Nevertheless, the parameter α cannot be chosen arbitrarily large, because the numerical evaluation of the terms X n (α, t) in (10) leads to overflow in finite precision.Consequentially, also in this case, the error keeps growing with time.For this reason, very often it is useful to speed up the convergence of the truncated series by applying the Euler transformation presented in [17] and [28].This results in adding a few additional terms to the original series, for the purpose of limiting K as much as possible for the convergence.Thus, (10) becomes Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where p is the number of terms of the Euler transformation and It is clear from (11) that the state vector variables must be sampled, for each t, on K + q different points on the complex plane, which share all the same real part α/t.Fig. 3 shows the first ten complex points per time sample over which the state vector has to be evaluated assuming α = 3.The trajectory is considered in the time interval [1], [2], [3], [4], [5] ns for illustrative purposes.
With the aim of clarification, the example regarding the FILT computation of the sin function employed in Section II-A is here reconsidered.In particular, in Fig. 4 it is evident that contrary to NILT, employing the parameters combination: (K = 20, α = 2) it is possible to reproduce satisfactorily the exact sin waveform in the overall time window.Unfortunately, it is also confirmed that it is not possible to choose α arbitrarily large (here α = 5 is chosen) to reduce the approximation error in (13) as much as possible.The sin inverse transforms in Fig. 5 complete the comparison of the two methods.The waveform is computed utilizing 12 terms for NILT (M = 24) and for FILT (K = 12, α = 2).It is evident that by fixing the number of terms for each series, the convergence is more easily reached by NILT when treating this kind of waveform (this is true also for typical circuit waveforms).
In conclusion, if long transients have to be evaluated, the FILT technique permits a successful evaluation in the overall time window (no matter how large it is), at the cost of an increase in the computational burden.If shorter transients are needed, employing NILT guarantees a faster convergence and, hence, a significant saving in computational costs.
It is just the case to observe that the two techniques can be combined to solve the same problem, exploiting the advantages of both of them.

III. COMPLEX FD PEEC MODELS FOR THE ACCURATE TD SIMULATION OF FAR-APART STRUCTURES
The PEEC method is an integral equations method that is able to provide equivalent circuits of EM problems which can be analyzed in both the time and FDs [14].The equivalent circuit is built through the definition of the partial elements, describing the elementary (or partial) electromagnetic interactions related to the EM model.After that, an appropriate 3-D discretization (mesh) of the structure under analysis is carried out.The continuously increasing complexity of modern electronic devices requires more and more full-wave descriptions of the EM phenomena related to the operation of the interconnects.Thus, it is desirable to include the propagation factors inside the model in order to reach a high degree of accuracy.Obviously, in the modeling of radiating devices, the full-wave representation is mandatory, especially when the EM response of two or more separated far-apart structures is needed.
In the standard PEEC formulation, the partial inductance [40] describes the magnetic interaction between the currents flowing in elementary volumes of the mesh.The full-wave partial inductance between two volumes V m , V n is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
defined in the Laplace domain as follows [14]: where ûm , ûn are the normal unit vectors associated with the cross sections S m and S n , r m and r n are the observation points in the two volumes over which the integration has to be performed and τ = |r m − r n |/c 0 is the time delay in the background medium characterized by a propagation speed c 0 .
Similarly, the electric interactions between the charges on two elementary surfaces A ℓ , A m of the mesh are described by the coefficient of potential defined, in the Laplace domain, as follows: The exponential e −sτ in ( 16) and ( 17) accounts for the propagation delay due to the finite value of the speech of light in the background medium.
The enforcement of Kirchhoff laws to the PEEC equivalent circuit leads to the following equations in the modified nodal analysis (MNA) form [41]: The unknowns in (18) are represented by the node potentials (s), and the branch currents I(s).Besides the partial elements matrices P(s), L p (s), the circuit model entails the incidence matrix A, the lumped elements matrix Y le (s) and the self impedance matrix Z(s).I s (s), V s (s) are the impressed node currents and branch voltages, respectively.
It is evident that the PEEC MNA representation (18), expressed in the Laplace domain, is of the form (1). Hence, the transient response of PEEC models can be recovered from the complex FD using the techniques presented in Sections II-A and II-B.

A. Minimum Delay Extraction for the Modeling of Separated Structures
In general, the TD response evaluated on a victim device due to a source located far away from it has to satisfy the causality principle [42] and, thus, it cannot occur before the minimum time delay t d defined as follows: being d min the minimum distance between the interacting objects and c the phase velocity in the background medium where they are located.If t d is too large (from tens to hundreds of ns) the approximate inverse transform solution of the far-away part of the system computed using NILT and FILT methods will be affected by a significant error unless a countermeasure is taken.To this purpose, the Bromwich integral in (2) can be recast through the introduction of the delay exponential term e st d with t ′ = t − t d > 0. It is then possible to apply the inversion techniques to the Bromwich integral considering the delayed kernel e st ′ .The NILT expression (5) becomes while, the FILT expression ( 14) can be rewritten as follows: where being h a generic index representing n or q.In conclusion, the effective inversion time is t ′ < t, and the inaccuracies introduced by long simulation times are dramatically reduced.

B. Underflow Issues in Complex FD PEEC Models (Far-Field Problems)
When the techniques already described are applied in modeling separate distant structures, it is necessary to be careful in the numerical evaluation of the delayed partial elements on the complex plane, especially if the distance between the objects is significant.Let us focus on a single generic partial element describing an interaction (magnetic/electric) between two mesh elements, each located on a different object.Their interaction assuming a centerto-center (CC) propagation delay [14] can be written in the complex FD as follows: where H 0 m,n describes the static interaction and is the propagation time delay in the background medium between the two elements, being R CC m,n their CC distance and c 0 the background medium (free-space for the standard PEEC method) phase velocity.The expression (24) has to be evaluated on different sets of points over the complex plane for each value of t ′ , using both the NILT and FILT techniques.The argument of the exponential can be written as follows: where z is a complex number that depends on the technique employed.For sufficiently large values of d min and small values of t ′ , it is very likely that the exponential e −sτ m,n reaches an underflow condition if the absolute value of its argument exceeds a fixed value K max .This value can be easily found through the knowledge of the smallest positive normalized floating-point number in IEEE double precision.Hence, in order to avoid the underflow of all the "far away" partial elements in the PEEC system, the following condition should apply: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Obviously, the latter condition is more severe at the first computation times, if the distance d min between the two elements is significant.The FILT method uses a set of poles whose real part is highly adjustable in this sense.Indeed, for each computation time t ′ , it is sufficient to enforce the condition still being able to maintain a reasonable accuracy at the first computation points.On the contrary, when employing the NILT method, the only way to modify the real part of the Padé poles at the first time samples t ′ is to reduce the order M. Unfortunately, also using small values for M does not avoid the underflow condition if the structures are very farapart.To better explain the underflow issues, Fig. 6 represents the behavior of the absolute value of the delay exponential function in (24) for two far-apart elementary volumes, when computed over the FILT points on the complex plane and over the Padé poles having a maximum real part, assuming M = 20.The exponential analysis is carried out considering the distances: 10, 20, 40 m, and varying the quantity t ′ = t − t d .It is clear that, for any of the three distances, the delay exponential reaches the underflow condition for a larger time interval when it is computed over the Padé poles rather than using the FILT complex points.When this case occurs, for sufficiently small values of t ′ , the partial elements describing the mutual interactions between two or more far-away volumes or surfaces become zero, compromising the accuracy of the computation.In Fig. 7, it is shown that using small values of order M, e.g., M = 4 for the NILT technique, the time interval where condition (28) is not matched and the underflow occurs is reduced, but is still present for the considered distances.
In this regard, despite the small values of the complex exponential in Figs. 6 and 7, it is important to point out that, in expressions ( 21) and ( 23), a compensation factor of the same order of magnitude as the complex exponential is introduced in the final solution due to the exponential term e st d .Hence, this observation further confirms the importance to avoid the underflow condition.

IV. COMPLEX FD MODELS FOR ELECTRICALLY LONG MTLS
Approximate Inverse LT techniques have been widely used for the TD characterization of transmission line (TL) structures [20], [43].Recently, an extension of the NILT method including nonlinear circuits has been proposed in [44].
Neglecting incident EM fields effects, the equations of MTLs composed by n conductors in the Laplace domain can be cast in a standard state-space form as follows [45]: where x is the line abscissa and is the state vector containing the n line voltages and currents.The state matrix Â(s) reads where Z(s) = R(s) + sL(s) and Y(s) = G(s) + sC(s).
It is well known that the end voltages and currents can be related by the so-called chain parameters matrix ˆ (L, s) as follows [45]: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where the 2n × 2n chain parameters matrix ˆ (L, s) is Pertinent boundary conditions need to be enforced at the line ends x = 0, L. Assuming linear time-invariant terminations, they read where Vs0 (s), VsL (s), Ẑ0 , and ẐL are the multiport Thevenin sources and impedances of the circuits connected to the MTL at x = 0 and x = L, respectively.The resulting problem described by ( 32) and ( 34) is well-posed for the end voltages and currents V(0, s), Î(0, s)), V(L, s), Î(L, s).
The transient port voltages and currents can be obtained through the application of inverse LT techniques to (32) and (34), upon the evaluation of the chain parameters matrix over appropriate points of the complex plane, depending on the chosen inverse transform technique.
The exponential of a matrix has been studied in depth for many years and can be computed in several different ways.Among the others, methods involving approximation theory, differential equations, matrix eigenvalues, and matrix characteristic polynomials have been proposed [46], behaving differently in terms of computational stability and efficiency.
In order to better illustrate the problems that may arise in computing (33) in (32), we consider, without loss of generality, the case of perfect conductors (R = 0) in a homogeneous lossy dielectric medium.The chain parameters submatrices can then be expressed in a closed form as follows [45]: where I n is the n × n identity matrix and Ẑc is the characteristic impedance matrix defined as follows: The coefficient γ (s) = (sµσ + s 2 µε) 1/2 is the propagation constant in the Laplace domain assuming known the electrical conductivity σ , the magnetic permeability µ, and the dielectric permittivity ε of the surrounding medium.
When the physical length of the lines is large (L > 10 m), the TD computation of the electric quantities in z = 0 can be very challenging.Indeed, the behavior of the complex hyperbolic functions in (35a) evaluated on the complex plane plays a crucial role.To explain this, we consider the generic complex point s = z/t, where t is the time evaluation value and z is a complex point that depends on the technique adopted to perform the inverse LT.In the limit of vanishing losses, the Fig. 8. Absolute value of the complex hyperbolic cosine computed for different line lengths (lossless medium) using NILT (M = 24, lines with markers) and FILT (adaptive α, lines without markers).
propagation constant can be written as γ (z/t) ≃ (z/t)(µε) 1/2 and, hence, the hyperbolic cosine function becomes being c = 1/(µε) 1/2 the propagation speed in the homogeneous dielectric medium.The function f c (s) reaches an underflow or an overflow condition if the absolute value of the exponential arguments exceeds the value K − or K + , respectively.This value can be found from the knowledge of the smallest and largest positive normalized floating-point number in IEEE double precision.Hence, to avoid a data loss, it is necessary that the following inequality is satisfied for each value of t: The condition (39) becomes more severe as t and c approach smaller values if the line length L is significant.The FILT method is more flexible in this sense since it is easy to enforce the condition providing an accurate TD solution at the initial computational points t.The same holds for the hyperbolic sinusoidal function.
Employing the NILT technique, the only way to try to satisfy the inequality (39) is to reduce as much as possible the Padé expansion number M, but unfortunately, with very small values of M it is not possible to match (39) when t is small.Hence, the TD computation of the MTL transients at z = 0, where the port quantities are not negligible even at the first instants, becomes very difficult.
In Fig. 8 it is shown the trend of the absolute value of function f c (s) when evaluated over the Padé poles assuming M = 24 and when computed over the FILT complex points complying with (40).The behavior is plotted for two values of the length of the line, L = 20 m, and L = 40 m.
It is evident that using the NILT method, significant portions of the time computational window are compromised,  as pointed out in Fig. 8.As expected, the time interval where the overflow occurs increases with the length of the line.On the contrary, the adaptive choice of α in the FILT technique allows outperforming the NILT method avoiding the overflow condition also for small values of t.This guarantees an accurate computation of the port transient voltages and currents at z = 0, even at the beginning of the computational time window.On the other hand, if the line length is intermediate (tens of meters) and the traveling signals are smooth enough, the reconstruction of far-end signals could be more convenient using NILT, for the convergence reasons explained in Section II .
To summarize the NILT and FILT features in TD MTLs modeling, a short resume is described in Table II.

A. Transmitting and Receiving Dipole Pair
To illustrate the concepts discussed in Section III, the two dipoles sketched in Fig. 9 are considered, whose geometric features are described in Table III.Precisely, the TD receiver voltage response is computed first considering the near-field interaction between the dipoles, and, subsequently, the receiver dipole is moved away from the transmitting dipole, in its farfield.For illustrative purposes, the transmitting dipole is driven by a voltage source in series to a resistor and we consider the receiver terminated on a load resistor.The signal source considered has a trapezoidal waveform with rising and falling times τ r = 3.2 ns and width τ w = 9.5 ns so that the spectrum is significant until f max ≃ 1 GHz.Taking this frequency as a reference, the dipoles are half-wave long.It is known that the transition from the near-field region to the far-field Fraunhofer region of the transmitter occurs at distances r such that [47] where r is the radial distance measured from the center of the transmitting dipole, D = 2ℓ is its maximum physical length, and λ = c 0 / f max is the wavelength in the free space at the maximum frequency of interest / f max .1) Near Field Analysis: In the near-field analysis, the distance between the two dipoles d is set as d = R f /2.The TD receiver voltage is computed through the inverse transform  techniques and compared to the results of a time-stepping reference solver [48] adopting the backward differentiation scheme of the second order (BD2).Fig. 10 shows the receiver voltage response assuming 50 terminations for both ports.In this case, the transient is quite long-lasting (more than 20 ns).The FILT solution is obtained by choosing K = 20 and p = 8 in ( 22) so that it is observed a good match between the FILT solution and the reference solution throughout the time window.For the Padé expansion-based inverse transform, the modified NILT technique (NILT2) [19] has been employed, to achieve a higher degree of accuracy.It is clear from Fig. 10 that the maximum exploitable expansion order, in this case, is M = 8, while larger orders cause the result to explode because of rounding errors in the residues K i that impact significantly the solution.In general, employing NILT2, this behavior is observed for relatively small values of M, because terms of the type: K 3 i , K 2 i are involved in the summation.Hence, it is not possible to achieve better accuracy using more terms in the series.For this reason, in this case, the recommended method is FILT, because it allows adjusting the accuracy by choosing the optimal number of series terms.The latter becomes an important feature when dealing with long-lasting transients.
Fig. 11 shows the receiver voltage response assuming 1 k termination for the two ports, where it is evident that the transient is more time-limited compared to the previous case.In this case, employing NILT2 with order M = 4 (two Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.series terms) ensures a very good accuracy over the entire time window, while for FILT the same accuracy is reached considering K = 10 and p = 5 in (22) for a total of 15 series terms.Hence, in this case, employing the NILT technique is more convenient in terms of computational efforts since a smaller number of system evaluations on the complex plane is sufficient to represent the response for each time sample compared to the FILT method.In particular, the cpu simulation time running NILT was 43.4 s, while it was 163.3 s for FILT, considering the same number of computation points for the two techniques.For completeness, the NILT2 and FILT relative errors at different time are in Table IV, assuming the same number of series terms (M = 4 for NILT2, K = 2 for FILT).It is seen that NILT2 outperforms FILT in terms of convergence.This confirms that employing NILT for moderately short transients offers the opportunity to minimize the computational burden.
2) Far-Field Analysis: In the far-field analysis, the receiver dipole is moved away from the transmitter to a distance d = 20 m.In this configuration, the minimum propagation time delay between the two structures is t d = 66.67 ns.The computation of the receiver voltage starts at a time t d = r/c 0 .The port terminations are assumed 50 .
As already explained, because of the underflow issues, the delayed PEEC model becomes inaccurate at the first computation time samples, when computed over the Padé poles using NILT.This is clear by observing a zoom of the receiver voltage response in Fig. 12, in which the NILT response is zero also after the propagation delay t d between the two dipoles.Coefficients of potential matrix pattern at the initial instants computed over the Padé poles.Fig. 13 shows the pattern of the coefficients of the potential matrix P(s), computed at one of the first computation instants (where the response is still wrongly zero) over the Padé poles.It is known that the coefficients of the potential matrix are full [14], but, when computed for small values of t ′ , the underflow issue causes the coefficients of potential describing the mutual interactions to be zero, resulting in null offdiagonal blocks.The same behavior is observed in the partial inductances matrix L p (s).
In the FILT series, with reference to (22), the K and p parameters have been set as K = 20 and p = 8.For the choice of parameter α, the adaptive criterion (28) has been adopted.
Since the source is piecewise linear (PWL), all the global responses are conveniently obtained, for each port, through the superposition of delayed versions of a unique ramp response.When the receiver voltage is computed, the wrong null initial portion affects considerably the final result.This is evident in Fig. 14, where a satisfactory agreement is observed between the FILT results and the reference time-stepping solution BD2, while the NILT results, obtained with M = 20, are inadequate.

B. Three-Phase Cable
In Section IV, it has been pointed out that the TD modeling of MTLs through inverse LT techniques is feasible and easy to implement since the formulation is closely connected to the classic FD representation in terms of per-unit-length (p.u.l.) parameters.In the second example, a three-phase cable is considered.Its cross section is represented in Fig. 15.It is assumed that the shield is a perfectly conducting cylinder, filled with a uniform dielectric material characterized by a dielectric relative permittivity ε r = 3.0.The line length is assumed as L = 10 m, which is a standard value employed in motor drive systems applications.The cable geometrical quantities are summarized in Table V.
The cable can be modeled as a six ports MTL system.The p.u.l.inductances matrix of the MTL model is analytically known [45] Since the conductors are placed in a homogeneous dielectric medium, the per-unit-length capacitance matrix is easily obtained as follows: The same holds for the per-unit-length conductance matrix being σ the equivalent dc conductivity of the dielectric material.
For example purposes, only one conductor is fed at one end through a 19 voltage source, while the other five ports are terminated on 19 passive loads.The voltage waveform is chosen as a double exponential pulse with α = 20 ps and β = 33 ps.The near-end crosstalk voltage induced in one of the victim lines is shown in Fig. 16.The focus should be put on the NILT solution behavior in the rising portion of the signal, shown in detail in Fig. 17, where it is evident that the results are unavailable for a significant part of the rising edge.The reason resides in the fact that the chain matrix reaches an overflow condition due to the very large real part of the Padé poles at the beginning of the analysis, e.g., at the first time samples.The converse is true for the FILT results, where, as described in the previous section, an adaptive choice of the real part of the FILT complex points permits an accurate representation at the early times.For completeness, In Fig. 18 the far-end voltage response is sketched.In this example, M = 24 (12 terms in the NILT series) was found to be the minimum NILT expansion number that guarantees an acceptable accuracy.For the same purpose, FILT analysis is performed by adopting K = 12 and p = 8 in ( 14), for a total number of 20 evaluations per time sample of the MTL system.As a consequence of the chosen number of  series terms, the simulation CPU time for the NILT technique was 2.67 s, while for FILT was 3.78 s, suggesting that if only the far-end responses are needed the NILT technique is the most convenient since it reaches an acceptable accuracy with a lower computational effort.On the contrary, for the computation of the near-end responses the FILT technique becomes necessary, since, due to its flexibility, it provides accurate results also at early times without losing important data.To combine the benefits offered by each technique, should not be excluded the possibility to combine the two approaches using FILT to support NILT where it is not able to provide data.

VI. CONCLUSION
Well-known numerical techniques can solve Maxwell's equations in the TD.Despite their robustness, they have some limitations primarily related to their stability and the modeling of dispersive media.A possible alternative is represented by the NILT that returns the TD responses inheriting the advantages of the complex frequency modeling.This work has presented a comparative analysis of two of the most popular methods for the computation of the numerical inverse LT, known as the NILT and FILT.Issues and limitations behind these two methods have been addressed and the possible solutions pointed out.Two numerical tests have been presented to corroborate the pros and cons of the two methods.

Fig. 4 .
Fig. 4. FILT-based computation of the sin function for two settings of the parameters.

Fig. 6 .
Fig.6.Absolute value of the complex exponential computed for different distances between two elementary volumes using NILT (M = 20, lines with markers) and FILT (adaptive α, lines without markers).

Fig. 7 .
Fig. 7. Absolute value of the complex exponential computed for different distances between two elementary volumes using NILT (M = 4, lines with markers) and FILT (adaptive α, lines without markers).

Fig. 10 .
Fig. 10.Receiver near field voltage response in the case of 50 terminations (example V-A).

Fig. 11 .
Fig. 11.Receiver near field voltage response in the case of 1 k terminations (example V-A).

Fig. 12 .
Fig.12.Zoom of the receiver far-field voltage response at the initial instants (example V-A).

Fig. 13 .
Fig.13.Coefficients of potential matrix pattern at the initial instants computed over the Padé poles.

TABLE I NILT
VERSUS FILT: A SUMMARY Finally, to resume the advantages and the limitations of the two techniques discussed in Sections II and III, a summary is reported in Table where, for different scenarios for distances and transient duration, the most suited technique is recommended.

TABLE II NILT
VERSUS FILT: A SUMMARY FOR MTLS

TABLE III GEOMETRIC
FEATURES OF THE TWO DIPOLES SYSTEM , exhibiting diagonal elements L s