International audienceBased on least-squares approximation of the rectangular pulse [1] by exponential functions, this paper presents an alternative method for performing numerical inversion of the Laplace transform. It compares favourably with the celebrated Vlach's method

Bréhonnet, Pascale

Calvez, Léon-Claude

Tanguy, Noël

Vilbé, Pierre

HAL-Université de Bretagne Occidentale

An Alternative Method for Numerical Inversion ofLaplace TransformsPascale Bre´honnet, Noe¨l Tanguy, Pierre Vilbe´, Le´on-Claude CalvezTo cite this version:Pascale Bre´honnet, Noe¨l Tanguy, Pierre Vilbe´, Le´on-Claude Calvez. An Alternative Methodfor Numerical Inversion of Laplace Transforms. IEEE Transactions on Circuits and SystemsPart 2 Analog and Digital Signal Processing, Institute of Electrical and Electronics Engineers(IEEE), 2006, 53 (6), pp.434-437. <hal-00176168>HAL Id: hal-00176168https://hal.archives-ouvertes.fr/hal-00176168Submitted on 21 Feb 2014HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestine´e au de´poˆt et a` la diffusion de documentsscientifiques de niveau recherche, publie´s ou non,e´manant des e´tablissements d’enseignement et derecherche franc¸ais ou e´trangers, des laboratoirespublics ou prive´s.This article has been accepted for inclusion in a future issue.> TCAS-II 1585<  1  Abstract— Based on least-squares approximation of the rectangular pulse [1] by exponential functions, this paper presents an alternative method for performing numerical inversion of the Laplace transform. It compares favourably with the celebrated Vlach's method.  Index Terms— Approximation, Laplace transform, numerical inversion,  rectangular pulse   I. INTRODUCTION ARIOUS methods have been proposed for numerical inversion of the Laplace transform. Among them, methods based on approximation of delta or exponential functions [2]-[4], on Fourier series [5] and on orthogonal expansions [6] have been considered. In this paper, we propose a particularly easy to use and robust method as an alternative to the well-known Vlach's method [3].  II. OUTLINE OF VLACH’S METHOD It is based upon rational approximation of the exponential function in the complex inversion formula ( ) ( ) ∫∫ ∞+ ∞−∞+ ∞− == jc jc zjc jc st dzetzFtjdsesFjtf 12121 ππ , where ( )sF  denotes the Laplace transform of the time function ( )tf . Replacing ze  by some rational approximant ( )zE  yields the following approximate inverse Laplace transform ( ) = ∑= tFattf iNi i γ11~ , (1) where ia  and iγ  are respectively the residues and the poles of the approximant ( )zE , either real or occuring in complex-conjugate pairs. In each new problem, only the evaluation of  Manuscript received May 7, 2004. Authors are with the Telecommunication Systems and Electronics Laboratory (LEST UMR CNRS 6165), Université de Bretagne Occidentale, CS 93837, 29238 BREST Cedex 3, France. (Corresponding author : P. Bréhonnet; e-mail: Pascale.Brehonnet@univ-brest.fr).  Copyright (c) 2006 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending an email to pubs-permissions@ieee.org   the Laplace transform has to be recalculated. A number of methods for obtaining optimal sets of constants ia  and iγ  have been investigated. The accuracy of the results depends on both the approximant of the exponential function ze  and of the function under consideration. There are many known methods for finding a rational function which is an approximant to ze . First, Vlach [3] used a Padé approximation with the same degree of numerator and denominator. A numerator lower by one degree has also been used. A drawback of Padé approximation has been unstability [4], especially for low-order approximations. In [4], the Padé approximant to ze  is written ( ) ( )( )zQ zPzE NM= . The main difficulty lies in the choice of the polynomial degrees M and N. For an increasing N, we find an increasingly more accurate approximant to the exponential function. But in this case, residues are increasing too, generating numerical inaccuracies.  For each N-order model, it is possible to have different solutions by varying the degree of numerator. Best results are generally obtained when the numerator degree is small, provided the denominator is stable. For each particular Laplace transform to be inverted, it is possible to determine the optimal degree M to be used with a given denominator ( )zQN , but this is a relatively long and tedious work. To apply Vlach’s method, the following steps must be respected [4]: 1) Choose a sufficiently high number of poles N and adjust M in such a way that ( )( )zQ zPtzF NM  has at least two zeros at infinity. Compute the corresponding approximant ( )tf~  to ( )tf  via equation (1) and tabulated values of ia  and iγ  [3]. 2) Solve other approximant ( )tf~  for same N and different M. If the solution differs considerably from the previous solution, keep changing M until the results stabilise for different M.  In practice, Vlach’s method is often carried out with { }2,1, −−∈ NNNM  which does not necessarily correspond to the best choice. It is worth noting that with increasing order N, numerical problems appear because the residues are too An Alternative Method for Numerical Inversion of Laplace Transforms P. Bréhonnet, N. Tanguy, P. Vilbé, and L.C. Calvez V This article has been accepted for inclusion in a future issue.> TCAS-II 1585<  2great. So, Vlach's method is not applied beyond order 30, with calculators using double precision. The problem of the choice of degrees of numerator and denominator is a great drawback because it depends on the Laplace transform under consideration. Moreover, to prove the efficiency of a method, examples are tested for some Laplace transforms for which the exact inverses are known but this is not the case in common engineering practise so, it is very interesting to dispose of several methods to ratify the results. The method proposed in this correspondence has proved its efficiency. It can be considered as an alternative to Vlach's method, without the above-mentioned drawbacks.  III. BASIS OF PROPOSED METHOD It is based on the approximation of the rectangular pulse by a linear combination of damped exponentials. Let ( )t∏  denote the causal rectangular function defined as ( )01=∏ t   ifif  [ ][ ]1.0 1.0∉∈tt . (2) Assuming [ [∞∈ .01Cf  has a continuous derivative f ′  on [ [∞.0 , we can write ( ) ( ) ( ) ( )000+−=′=′Π ∫∫∞ ftfdfdxxftx t ττ  (3) which may be rewritten as  ( ) ( ) ( )∫∞ ′Π++= 00 dxxftxftf . Replacing the rectangular function by some approximant  ( ) ∑ = −=Π Nk tk keat 1~ γ , with the real part ( ) 0>ℜ kγ , and integrating by parts yield the following approximant to ( )tf  ( ) ( ) ( )∑=  +−++= Nk kkk ftFtaftf 1 00~ γγ  (4) provided ( ){ }kkf t γσ ℜ< min , where fσ  denotes the abscissa of convergence of the Laplace integral. Notice that ( ) ( )00~ +=+ ff  is guaranteed by (4), on account of the initial value theorem for the Laplace transform. To preserve the final value when it exists, the constraint 11 =∑ =Nk ka  has to be imposed. In this case, equation (4) reduces to ( ) ∑= = Nk kkk tFtatf 1~ γγ  (5) in which only the knowledge of  the constants ka , kγ relative to an approximant to the rectangular pulse are required to compute ( )tf~  from ( )sF . Like in Vlach’s method, determination of the poles of the Laplace transform to be inverted has been avoided. Residues ka  [Table I] are significantly smaller than in Vlach’s method [3] for a same order (the ratio which is about 510  for 10=N  increases with N) which makes our method less sensitive to numerical inaccuracies (see Fig. 3 and Fig. 4).  The preciseness of formula (5) depends on the accuracy of the approximation of the function ( )t∏ . In [1], an optimization technique is provided for the least-squares approximation of the rectangular pulse, owing to a near optimal starting point, the integrated squared error is reduced through few iterations of a Gauss Newton process, beyond all usual orders (up to order 200). This great accuracy is particularly useful for determining the required constants ( )kka γ,  which must be computed only once for a given order N and can be stored for permanent use.  IV. ILLUSTRATION We have checked many examples to illustrate the applicability of the suggested numerical inversion technique. For some of them, Vlach's method gives also good results but usually after several tests. First consider the irrational Laplace transform  ( ) ( )ssshsF += 11, which is the normalised transfer function of a RC distributed circuit. Fig. 1-2 represent the time signal obtained by Vlach's method and by our method; in both cases, the accuracy is acceptable. Consider the irrational Laplace transform  ( )sseesF −−+= 231 23  with a view to computing the time signal ( )tf , which is a very severe test. The numerical differences between the methods are shown in Table II and Fig. 3-6. In Vlach’s TABLE I  Negatives of Poles – Present method Residues – Present method 3.65623977941246 ±  2.33980675356003 j 3.04259335096752 ±  7.47702946613745 j 2.47879760012144 ±  12.94030319544558 j 1.94081484142228 ±  18.38586389094272 j 1.22236607997850 ±  23.52676351409633 j -2.02699639954484 ±  6.5556240442333 j 2.21838424179506 ±  0.27882389516450 j 0.48583161398352 ±  0.62262273410175 j -0.07636822980337 ±  0.30776792623453 j -0.10085122644020 ±  0.02574141374327 j  This article has been accepted for inclusion in a future issue.> TCAS-II 1585<  3method, for an N-order approximant to ze , M represents the most appropriate degree of the numerator obtained after several tests. On the other hand, for a given N order, our method yields directly the approximate time signal. In both cases, maxt  denotes the maximal time for which accuracy is acceptable. In practice, we used  22max ~~2ffftacceptable′Π−Π−≈ . 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.30.40.5 time transient response Fig. 1 : Vlach’s method, order10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.30.40.5 time transient response Fig. 2 : Present method, order 10  0 1 2 3 4 5 6 7 8 9 10−20−15−10−50510152025 timetransient response Fig. 3 : Vlach’s method, order 30, M = 29  0 1 2 3 4 5 6 7 8 9 10−20−15−10−505101520 timetransient response Fig. 4 : Vlach’s method, order 30, M = 25  0 1 2 3 4 5 6 7 8 9 10−20−15−10−505101520 timetransient response Fig. 5 : Present method, order 30  This article has been accepted for inclusion in a future issue.> TCAS-II 1585<  40 1 2 3 4 5 6 7 8 9 10−20−15−10−505101520 timetransient response Fig. 6 : Present method, order 60   The third example deals with the Laplace transform ( )sF  of an oscillating signal. For this example, the exact inverse is known to be the Bessel function 23J , so the quadratic error can be calculated. The normalised quadratic error criterion is defined as feQn =ˆ  where ffe −= ~ˆ  and  ( )∑ == nk ktff 1 22 ˆ . The time function is calculated up to st 150=  with a sampling period sT 05.0= , i.e. 3000=n . The best result obtained with Vlach's method corresponds to ( ) ( )30,26, =NM  and 4396.0=nQ ; in this case, the original signal can be built up to the th8  oscillation (Fig. 7). Our technique used with a 30-order model, yields directly an acceptable signal up to the th14  oscillation with 3569.0=nQ  (Fig. 8). Furthermore the time interval can be significantly extended using our method and a 60-order approximant to the rectangular pulse (Fig. 9). In this case, we obtain a quadratic error 01058.0=nQ . 0 50 100 150−0.4−0.3−0.2−0.100.10.20.30.40.50.6 timetransient response Fig. 7 : Vlach’s method, order 30, M = 26, Qn = 0.4396  0 50 100 150−0.4−0.3−0.2−0.100.10.20.30.40.50.6 timetransient response Fig. 8 : Present method, order 30, Qn = 0.3569  0 50 100 150−0.4−0.3−0.2−0.100.10.20.30.40.50.6 timetransient response Fig. 9 : Present method, order 60, Qn = 0.01058  V. CONCLUSION A numerical method particularly useful for time responses TABLE II  Method order maxt  Vlach 30, 29=M  1.5 Vlach 30, 25=M  2.5 Present method 30 3.5 Present method 60 8.5 This article has been accepted for inclusion in a future issue.> TCAS-II 1585<  5of high or infinite dimension systems is given in this paper. It is based on the approximation of the rectangular pulse. The poles and residues of the approximant were computed with high accuracy. In practice, a 30-order approximation is usually sufficiently accurate for most engineering applications but higher approximations are easily available if necessary. The technique does not require the determination of poles of the Laplace transform to be inverted which represents a great advantage, and only one iteration is necessary. So, the method can be implemented very easily in a program. For a given precision of computer, the residues [Table I] being significantly smaller than with Vlach’s method [4], the time signal can be built on a larger interval. This conclusion is supported by many numerical experiments demonstrating that our method is simple and efficient. As is well known, every Laplace transform numerical inversion method breaks down for some functions and therefore checking by different methods can greatly increase confidence in the results achieved.  REFERENCES [1] P. Vilbé, S. Azou, P. Bréhonnet, and L.C. Calvez, ''The L2-optimal time-delay rational  Laplace model revisited'', IEEE Trans. on Circuits and Systems, 1999, CAS1-46, (4), pp. 505-512. [2] V. Zakian, ''Numerical inversion of Laplace transform'', Electron. Lett., 1969, 5, (6), pp. 120-121. [3] J. Vlach, ''Numerical method for transient responses of linear networks with lumped, distributed  or mixed parameters'', J. Franklin Inst., 1969, 288, pp. 99-113. [4] K. Singhal, J. Vlach, ''Program for numerical inversion of the Laplace transform'', Electron.  Lett., 1971, 7, (14), pp. 413-414. [5] F.R.de Hoog, J.H. Knight and A.N.Stokes, ''An improved method for numerical inversion of  Laplace transforms'', SIAM J. Sci. Stat. Comput., 1982, 3, (3), pp. 357-366. [6] J.L. Wu, C.H. Chen, and C. F. Chen, ''Numerical inversion of Laplace transform using Haar wavelet operational matrices'', IEEE Trans. on Circuits and Systems, 2001, CAS1-48, (1), pp. 120-122    

An Alternative Method for Numerical Inversion of Laplace Transforms

https://hal.archives-ouvertes.fr/hal-00176168

An Alternative Method for Numerical Inversion of Laplace Transforms

Abstract

Similar works

Full text

Available Versions

HAL-Université de Bretagne Occidentale