Numerical Inverse Laplace Transforms for Electrical Engineering Simulation

Numerical inverse Laplace transform (NILT) methods are widely used in various scientific areas, especially for a solution of respective differential equations. In field of an electrical engineering many various approaches have been considered so far, but mostly for a single variable (1D NILT), see at least (Brančík, 1999, 2007b; Cohen, 2007; Valsa & Brančík, 1998; Wu at al., 2001) from plenty of papers. Much less attention was paid to multidimensional variable (nD NILT) methods, see e.g. (Hwang at al., 1983; Singhal at al., 1975), useful rather for more complicated electromagnetic systems. The 2D NILT methods, see e.g. (Brančík, 2005, 2007a, 2007b; Hwang & Lu, 1999), can be applied for a transmission line analysis, or nD NILT methods, n ≥ 2, for a nonlinear circuits analysis, if relevant Laplace transforms are developed through a Volterra series expansion, see e.g. (Brančík, 2010a, 2010b, Karmakar, 1980; Schetzen, 2006), to highlight at least a few applications. This paper is focused on the class of NILT methods based on complex Fourier series approximation, their error analysis, their effective algorithms development in a Matlab language, and after all, on their selected applications in field of electrical engineering to show practical usefulness of the algorithms.


Introduction
Numerical inverse Laplace transform (NILT) methods are widely used in various scientific areas, especially for a solution of respective differential equations.In field of an electrical engineering many various approaches have been considered so far, but mostly for a single variable (1D NILT), see at least (Brančík, 1999(Brančík, , 2007b;;Cohen, 2007;Valsa & Brančík, 1998;Wu at al., 2001) from plenty of papers.Much less attention was paid to multidimensional variable (nD NILT) methods, see e.g.(Hwang at al., 1983;Singhal at al., 1975), useful rather for more complicated electromagnetic systems.The 2D NILT methods, see e.g.(Brančík, 2005(Brančík, , 2007a(Brančík, , 2007b;;Hwang & Lu, 1999), can be applied for a transmission line analysis, or nD NILT methods, n ≥ 2, for a nonlinear circuits analysis, if relevant Laplace transforms are developed through a Volterra series expansion, see e.g.(Brančík, 2010a, 2010b, Karmakar, 1980;Schetzen, 2006), to highlight at least a few applications.This paper is focused on the class of NILT methods based on complex Fourier series approximation, their error analysis, their effective algorithms development in a Matlab language, and after all, on their selected applications in field of electrical engineering to show practical usefulness of the algorithms.

Multidimensional numerical inverse Laplace transform
An n-dimensional Laplace transform of a real function f(t), with t = (t 1 ,...,t n ) as a row vector of n real variables, is defined as (Hwang at al., 1983) 1 00 () = () e x p ( ) where s = (s 1 ,...,s n ) and T means a transposition.Under an assumption |f(t)| < Mexp(αt T ), with M real positive and α = (α 1 ,...,α n ) being a minimal abscissa of convergence, and the nD Laplace transform F(s) defined on a region {s ∈ C n : Re[s] > α}, with c = (c 1 ,...,c n ) as an abscissa of convergence, and the inequality taken componentwise, the original function is given by an n-fold Bromwich integral (2) In the papers (Brančík, 2007a(Brančík, , 2007b(Brančík, , 2010b)), it was shown for the 1D, 2D, and 3D cases, the rectangular method of a numerical integration leads to an approximate formula whose a relative error is adjustable, and corresponds to the complex Fourier series approximation of a respective dimension.The method has been generalized for an arbitrary dimension n in the recent work (Brančík, 2011).

Complex Fourier series approximation and limiting relative error
Substituting s i = c i + jω i into (2), and using a rectangular rule of the integration, namely ω i = m i Ω i , and Ω i = 2π/τ i as generalized frequency steps, with τ i forming a region of the solution t ∈ [0,τ 1 ) ×...× [0,τ n ), an approximate formula is with As is shown in (Brančík, 2011), a limiting relative error δ M of (3) can be controlled by setting c = (c 1 ,...,c n ), defining paths of the integration in (2), namely for i = 1,…,n, and while keeping the equalities τ 1 (c 1 -α 1 ) =...= τ n (c n -α n ).The simplification in ( 4) is enabled due to small values δ M considered in practice.The last equation is used for setting up parameters of the nD NILT method relating them to a limiting relative error δ M required for practical computations.

Practical computational methods
It should be highlighted that the formula (4) is valid, and a relative error supposed is really achievable by the nD NILT (3), if infinite numbers of terms are used in the series.In practice, it cannot sure be fulfilled, but a suitable technique for accelerating a convergence of infinite series is usable, as is e.g. a quotient-difference (q-d) algorithm (Rutishauser, 1957).Besides, as has been already successfully used for cases of n ≤ 3, the formula (3) can be rearranged to enable using FFT & IFFT algorithms for an effective computation.

Partial ILTs evaluation technique
The technique of practical evaluation of the n-fold infinite sum (3) follows the properties of the n-fold Bromwich integral (2), namely we can rearrange it into the form ( ,,,) = (,,,) Although the order of the integration may be arbitrary on principle, here the above one will be used for an explanation.Similarly, (3) can be rewritten as with ,) .
As is obvious we need to use a procedure able to make the inversion of Laplace transforms dependent on another n-1 parameters, complex in general.Let us denote arguments in (8) by p i = (p 1 ,...,p n-1 ,p n ).Then the ILT of the type can be used n times, i = n,n-1...,1, to evaluate (8), with p n = (s 1 ,...,s n-1 ,s n ), p n-1 = (s 1 ,...,s n-1 ,t n ),..., p 1 = (s 1 ,...,t n-1 ,t n ), and p 0 = (t 1 ,...,t n-1 ,t n ), while p j = s j for j ≤ i, and p j = t j otherwise.A further technique is based on demand to find the solution on a whole region of discrete points.

FFT, IFFT, and quotient-difference algorithms utilization
As is shown in (Brančík, 2007a(Brančík, , 2010c)), the discretized formula (10) can be evaluated by the FFT and IFFT algorithms, in conjunction with the quotient-difference (q-d) algorithm for accelerating convergence of the residual infinite series, see following procedures.
To explain it in more detail, let us consider an r-th cycle in gaining the original function via (9), i.e.
The above stated formula can be decomposed and expressed also as www.intechopen.com11 0 11 00 0 0 where individual terms are defined as As is evident the first and the third finite sum of ( 12) can be evaluated via the FFT and IFFT algorithms, respectively, while 2P+1 terms from the infinite sums are used as the input data in the quotient-difference algorithm (Macdonald, 1964;McCabe, 1983;Rutishauser, 1957).We can replace the above infinite power series by a continued fraction as which gives much more accurate result than the original sum truncated on 2P+1 terms only.The q-d algorithm process can be explained based on a lozenge diagram shown in Fig. 1.The first two columns are formed as while the successive columns are given by the rules Then, the coefficients d m , m = 0,...,2P, in ( 14) are given by For practical computations, however, the recursive formulae stated below are more effective to be used (DeHoog et al., 1982).They are of the forms for m = 1,...,2P, k ∀ , with the initial values A -1 = 0, B -1 = 1, A 0 = d 0 , and B 0 = 1.Then, instead of the continued fraction ( 14), we can write The q-d algorithm is a very efficient tool just for a power series convergence acceleration, here enabling (7) to achieve a relative error near its theoretical value defined by ( 4), see the following examples.

Matlab listings and experimental errors evaluation
In this part experimental verifications of the nD NILT theory above will first be presented, for one to three dimensional cases, i.e. n ≤ 3.For such dimensions the Matlab functions have been developed and errors stated on a basis of some sample images with known originals.The Matlab listings of basic versions of the NILT functions are provided, together with examples of their right calling.Another Matlab listings will be discussed in more detail later, in the chapter with practical applications.

One-dimensional NILT
In case of the 1D inverse LT, a well-known Bromwich integral results from (2), namely where indexes 1 were omitted.By using the theory above a path of the numerical integration is stated according to (4), leading to In contrast to most other approaches, the 1D NILT method described here enables to treat complex images resulting in complex originals as no real or imaginary parts are extracted during an evaluation process.It can be useful in some special applications, not only in the electrical engineering.We can śhow it on a simple transform pair () Of course, when preprocessing the transform to arrange it to a Cartesian form, as is shown on the right sides in ( 22), the result could by get by inverting the real and imaginary parts separately, by using an arbitrary NILT method.Here, however, no symbolic manipulations are needed in advance, and F(s) enters the NILT function in its basic form as a whole.
A Matlab language listing is shown in Tab. 1, where the relative error needed is marked by Er and is subject to a change if necessary, similarly as the minimal abscissa of convergence (exponential order) α , alfa, numbers of points for the resultant solution, M, and for the q-d algorithm, P.  Another test functions are considered in Tab. 2, with numerical results shown in Fig. 3.As is again obvious from Fig. 3 and by using the theory above the paths of numerical integrations are stated based on (4) as A Matlab language listing is shown in Tab. 3, with all the parameters denoted by similar way as in the previous case.Nevertheless, the Laplace transform variables and the original variables have changed in this listing as s 1 → p, s 2 → q, and t 1 → x, t 2 → y, respectively, which simplified a writing.The parameters are then indexed in compliance with these new notations.Besides, numbers of points used to plot three-dimensional graphical results are set by xpl and ypl.
For the same reasons as at the 1D NILT, the 2D NILT method discussed here enables to treat complex images of two variables resulting in complex originals.We will śhow it on a simple transform pair The 2D NILT function is called from a command line as follows: nilt2c('F',xm,ym); where F is a name of another function in which the F(p,q) is defined, and xm and ym mark upper limits of the original variables x and y.In our case, and for ω 1 = ω 2 = 2π, this function can have a form function f=exp2c(p,q) f=1./(p-2*pi*j)./(q-2*pi*j); and its calling can look like nilt2c('exp2c',3,3); with graphical results in Fig. 4. As the originals are bounded by values ±1, and α 1 = α 2 = 0, we can see the errors satisfy (21)   very well (δ M = 10 -8 was considered), excluding beginnings of the 2D region.

Three-dimensional NILT
In case of the 3D inverse LT, a three-fold Bromwich integral results from ( 2 Here only experimental results will be shown to verify an accuracy of the method.A Matlab language listing looks similarly like for the 2D NILT case, but the partial NILT subfunction is called once more, and respective arrays dimensions are enlarged.Original functions corresponding to 3D Laplace transforms cannot be displayed graphically as a whole, of course.However, for one variable chosen as constant, it is posssible to display three respective two-dimensional cuts.It will be demonstrated on the example of 3D shifted unit step, with a Laplace transform pair with different values of shifts along respective coordinates so that correctness of results can easily be identified, see Fig. 6.Errors again correspond to theoretically expected ones.Error Fig. 6.Numerical inversion leading to shifted 3D unit step f(t 1 ,t 2 ,t 3 ) = 1(t 1 −1,t 2 −2,t 3 −3)

Application of NILT algorithms to electrical engineering simulation
In this chapter some examples of the application of the NILT algorithms developed relating to problems of electrical engineering simulation are presented.First, the 1D NILT method is applied for the solution of transient phenomena in linear electrical circuits with both lumped and distributed parameters.This well-known approach is usable wherever linear ordinary differential equations (ODE) are transformed into algebraic ones so that an inverse Laplace transform can be considered.Then the 2D NILT method is utilized to solve transient phenomena on transmission lines (TL) after relevant telegraphic equations (a type of partial differential equations (PDE)) are transformed into algebraic ones by a 2D Laplace transform.In this way voltage and/or current distributions along the TL wires can be determined in a single calculation step.Finally, the utilization of the 1D to 3D NILTs to weakly nonlinear electrical circuits solution is discussed.In this case the relevant nonlinear ODEs describing the circuit are expanded into Volterra series which respective NILTs are applied on.

One-dimensional NILT algorithm application 3.1.1 Preliminary example based on lumped parameter circuit
A simple example demonstrating the application of the basic 1D NILT algorithm in Tab. 1 is shown in Fig. 7.This really initiatory linear electrical circuit was chosen with an intention to be also considered later, in chapter 3.3.1,as a nonlinear circuit, with G 2 being a nonlinear element.In this way one will be able to compare results and make some conclusions.Here one more parameter depict is used to define a method of plotting individual items from a set of originals.The 1D NILT function is called as niltcv('F',tm,'depict'); where 'depict' is a text string 'p1', 'p2', or 'p3', see Tab. 6 for more details.

Application for transmission line simulation
Here, the 1D NILT algorithms will be used to simulate voltage and/or current distributions along transmission lines (TL), as shown on a Laplace-domain TL model in Fig. 9.As is well known, this model results from the application of a Laplace transform, with respect to time, on a pair of partial differential equations (telegraphic) of the form with R 0 , L 0 , G 0 , and C 0 as per-unit-length (p.-u.-l.) parameters, being constant for uniform TLs, and with a length l.
When considering zero initial voltage and current distributions, v(0,x) = 0 and i(0,x) = 0, and incorporating boundary conditions, we get the Laplace-domain solution in the forms In this case, the 1D NILT algorithm in Tab. 5 is called as niltcv('Vsx',2e-8,'p3'); that is the plott3 function is used for the plotting, see Tab. 6, and a time limit is half of that in Fig. 10 to get well-observable results.Again, the current distributions can be gained via the above function slightly modified according to (35).Both 3D graphs are depicted in Fig. 11.

Conclusion
The paper dealt with a specific class of techniques for the numerical inversion of Laplace transforms, namely based on a complex Fourier series approximation, and connected with a quotient-difference algorithm to accelerate the convergence of infinite series arising in the approximate formulae.The 1D to 3D NILT techniques have been programmed in the Matlab language (R2007b), and most important ones provided as the Matlab function listings.To guide readers all the algorithms were explained on selected examples from field of electrical engineering, including right callings of the functions.In contrast to most others the NILT methods here developed are utilizable to numerically invert complex Laplace transforms, leading to complex originals, which can be useful for some special purposes.As has resulted from error analyses the accuracies range relative errors from 10 -8 to 10 -10 without difficulties which is acceptable for most of practical needs.Based on Matlab functions presented, one could further generalize a vector version of the 1D NILT function towards a matrix one, enabling e.g. to simulate multiconductor transmission line systems, as is shown in (Brančík, 1999), where, however, an alternative technique, socalled ε algorithm, has been applied to accelerate the convergence of infinite series.
According to the author's knowledge, the paper presented ranks among few summary works describing multidimensional NILT techniques, covering Matlab listings beyond, based just on a complex Fourier series approximation, and in conjunction with the quotient-difference algorithm, which seems to be more numerically stable compared to the ε algorithm mentioned above. www.intechopen.com been considered, and τ r = M r T r .

Table 1 .
If only real transforms F(s) are considered the bottom line in the listing can be inactivated.The NILT function is called from a command line as follows: niltc('F',tm); where F is a name of another function in which the F(s) is defined, and tm marks an upper limit of the original variable t.In our case, and for ω = 2π, this function can have a form Matlab listing of 1D NILT method accepting complex arguments As is obvious, the Laplace transform must be defined to enable Matlab array processing, i.e. element-by-element array operators have to be used.Thus, the calling our function can look like niltc('expc',4); if the function is saved under the same name, expc, or it is placed inside the M-file with own NILT function (Tab.1), following always its body.Alternatively, the calling can look like [ft,t]=niltc('F',tm); if respective variables in the brackets are to be saved in the memory after the function finishes.

Table 2 .
the relative errors satisfy theoretical expectations, with an exception of vicinities of discontinuities.Test Laplace transforms for errors evaluation