Application of best rational function approximation for Laplace transform inversion

Application of best rational function approximation for Laplace transform inversion I. M. Longman (*)

ABSTRACT The application of a method of least squares Laplace transform inversion due to the author is discussed, and examples are given, which also serve to give us best trigonometric and exponential function approximations to some known functions. In particular the unit step function H(t-1) is considered, and the "best" approximations obtained for it would seem to have application in electrical network theory, including the design of delay lines. 1. INTRODUCTION In a recent paper the author [1] has developed a technique for approximate Laplace transform of a function f-(P)

7 e-Ptf(t)dt (1) o of the Laplace transform operator p. Essentially what is done is to approximate f(p) b y an nt.h order rational function approximant n

g-n(P) = Z A r / ( p + ~r)

(2)

where the n poles p = - % are supposed to be all different and not to lie in the right-hand half plane Re (p-) > o, so that we assume Re ( a r ) > 1 o ,

r = 1 , 2. . . . .

(3)

(4)

of gn (P) is made to approximate f(t) in a "best" manner by demanding t h a t oo

(5)

0 be a minimum, More precisely for a fixed "weight" w, whose choice in advance is based on a number o f considerations, the A r s and arS are chosen to minimize I n . An essential feature of the m e t h o d is that apart from a constant (independent o f the Ars and ¢~rs) In is expressible directly in terms o f f'(p) which is presumed to be a known function o f p. In fact we have n

r=l k=l

ArA k

0tr+rVk+W

+ K(w)

n

2 z ari~-(ar+W) r-_l

(6)

where K (w)

=7

e - W t I f ( t ) ]2 dt

f(t) = E l ( t ) = f (e-U/u) du.

(9)

t

In the present paper the method o f minimization o f I n (w) is further discussed, particular attention being paid to (the frequent) cases where Ars and ~ ' s occur in complex conjugate pairs although o f course I n (w) and gn (t) are real. As an example we consider the case f-(P) = e - P / p

(10)

H(t-1) =0

t 1 )"

(11)

Also considered is the case ~-(p) = ( p 2 + 1)-1/2

(12)

for which

In(w) = f e-Wt[gn(t) - f(t)]2dt

z

(8)

for which the inverse is the exponential integral

1

n -art gn (t) = ~ Are r=l

n

f-(v) = (l/p) log (1 + p)

for which the exact inverse is the step function n.

The inverse

In(w) = X;

The method has been applied by the author [1] to the example

(7)

O

does not depend on the A~s and arS or even on n.

f(t) = Jo (t),

(13)

the Bessel function o f the ftrst kind of order zero. 2. THE MINIMIZATION OF In(w ) From equation (6) it is evident that I n (w) is a function o f N = 2n parameters Ar, a r (r = 1 . . . . . n). The method o f minimization adopted is that due to Fletcher and PoweU [2] for finding unrestricted local minima o f a function ¢ ( x 1, x 2 . . . . . x N)

(14)

o f N real variables x I . . . . . x N. In order to apply this to In(W ) in which Ar's and ar'S may occur in complex conjugate pairs, these parameters were first expressed in terms o f real x's according to the scheme A 1 = x 1 + ix 2

a 1 =Xn+ l+ixn+

A 2 = x 1 - ix 2

'~2 = Xn+l - iXn+2

A 3 = x 3 + ix 4

a 3 = Xn+ 3 + iXn+ 4 (15a)

A 4 = x 3 - ix 4

a4 = Xn + 3

-

2

iXn+ 4

Department o f Environmental Sciences Tel-Aviv University, Israel Journal o f Computational and Applied Mathematics, volume I, no 1, 1975.

17

where if n is even we end with An_ 1 = Xn_ 1 + ix n

given by (4) we can construct a possible g n + l ( t )

an_ 1 = X2n_ 1 + iX2n

n+l

gn+l(t) = z A;.e An

= Xn_ 1 - ix n

an

= X2n_ 1 - iX2n

-art

r___l

(15b)

by taking

while if n is odd we end with a r = ar

An = xn

a n = X2n

r =

1, . . . , n

(15c) an+l=an=an

The above scheme also allows for the possibility of some (or all) of the Ars, arS being real. For example if x 2 and Xn+ 2 are zero we have A I =A 2=x I

a I =a 2=xn+

1

Ar' = Ar

r= 1,...,n-1

and any choice o f A n, A n + 1 satisfying A n + A n + 1 = A n.

so that A1, a I are real. The presence of A2, a 2 does not disturb, since then, as we can see from equation (2) or equation (4) we can absorb two terms in g-n (P) or in gn (t) having the same a into one term. For the procedure o f Fletcher and Powell we need to start with an initial vector

Evidently then g n + l ( t ) is identical with gn(t), and of course g n + l ( P ) is identical with gn(P) as given by (2). This type of local minimum which is characterized by two or more equal a r s is o f course immediately detectable.

x = (Xl, X2 . . . . . ,XN) and proceed in steps towards the minimum of ~ using also the vector

3. CHOICE OF THE WEIGHT w In the integral (5) the value o f w is more or less at our disposal The expression "more or less" is used since for some functions f(t) the integral (7) with w = 0 will not converge, and then we cannot be sure that our minimization procedure (5) has meaning. This is the case for example when q t ) = Jo(t). For such cases we can often ensure convergence by using a positive value of w. The value of w also has a role in emphasising the fit o f gn(t) to f(t) for small values o f t, since a large positive value o f w deemphasises the contribution o f the large values o f t to In(w ) . On the other hand if (7) converges with w = 0, we may prefer this value o f w in order to get the "best" fit over the whole range o f t. In some cases, e. g. (8) even negative values of w may give results o f interest. At the expense o f only a relatively slight complication o f our formulas we can even employ a "window" ( e - W l t - e - w 2 t ) in place of e -wt in the integral (5). w 1 and w 2 can be adjusted to emphasise a particular range o f values of t. It would be desirable to choose w 2 > w I so that the integral (5) remains positive. In certain problems, for example when r ( p ) is given by (10), f(t) may tend to a constant c when t -~ 0% so that (7) diverges for w = 0. We may nevertheless still obtain the "best overall fit" advantage o f w = 0 by subtracting c from f(t) (which amounts to subtracting c/p from it(p)), doing the inversion, and finally adding back c to our result. This method was used in the approximate inversions o f (10) below. It is o f interest here to note that if we solve our problem for given f-(p) and w, we have also automatically solved the problem for f ( p + b ) (having inverse e - b t f ( t ) ) w i t h weight w - 2 b , if we add b to our arS. This. is evident on inspection of equations (6), (7). In particular from the solution At, a r (r = 1 . . . . . n) for ~-(p) using weight w = 0, we obtain immediately the solution for ~-(p-a)(having inverse eatf(t))with weight w = 2a, ff we subtract a from our a r s.

g = ' grad ~. At each stage, then, from the current vector x__we calculate the ArS and a r s from (15), and hence obtain $ = In(w ) from (6), and the vector g by appropriate (analytical) differentiations of equation (6). The procedure was programmed for the computer, and results obtained for a number of examples. In cases where K(w) is known we are able to give the value o f rain In(w ). In other cases, however, we are still able to follow the relative changes in In(w ) as the minimization proceeds. Also since K(w) is independent o f n we can still record the changes in the minimum value as we increase n. It is as well to point out here that our minimization procedure (depending upon our starting vector x ) is liable to give us a local minimum o f In(w ) but not necessarily the absolute minimum. This however will not in general disturb us since we make a comparison o f results for different values of n, and if on proceeding from n to n + l our minimum increases we should assume that we do not have the absolute minimum o f I n +1 (w). Then we can either use a different starting vector x t o search for the absolute minimum, or simply proceed to n+2. Such increase in minimum as we go from n to n + l is immediately detectable o f course even ff we do not know the value o f K(w), since as noted above K(w) is independent o f n. In our first example given below such a local minimum was found during the course of the calculations for the case n = 7, and is included in the results as a matter o f interest. In any case a decreasing sequence {min I n } is still useful even in the (unlikely) case that each member is not an absolute minimum. It is o f interest to mention that some local minima o f In(w ) may correspond to the absolute minimum o f Im(w ), where m < n. For example from gn(t) as

Journal o f Computational and Applied Mathematics, volume I, no 1, 1975.

18

4. ACCURACY There is not known to the author a theorem stating that lim rain I n (w) = 0 1-1--).oo

or under what conditions this will be so. However in view o f the high accuracy usually achieved by the present author [1], [3]-[9] and by other workers (Luke [10], [11], [12, pp. 255-268], Akin and Counts [13], [14], Counts and Akin [15], Murphy [16])in the application o f rational approximations to Laplace transform inversion, we expect even better results using approximations designed to be the "best" for this purpose. Even the proof o f such a theorem is far from being all that we need. More important is the rate o f convergence, and our method is useful if for small values o f n, say n < 10 or even smaller, we can attainhigh accuracy. At any rate it is evident that abs min I n (w) must decrease as n increases, and since all In(w ) are bounded below by zero, it follows that the above limit must exist. We may reasonably conjecture that for most "smooth" functions f(t) the limit will be zero. Of course the most superficial consideration must convince us that discontinuous functions are not readily approximated by continuous ones, so that if f(t) is discontinuous, as in the first example below, we should not expect high accuracy by this method of inversion. If the discontinuities of f(t) are known in magnitude and location (as they usually are) from physical or other considerations, then they should be subtracted out o f f(t) before we start the inversion, and then replaced afterwards. More precisely, we subtract from f-(p) the Laplace transform o f a staircase function constructed from the discontinuities o f f(t), and then correct our approximate inverse afterwards. The second example in which qt) : Jo(t) shows the very high accuracy which is readily attainable when f(t) is continuous. 5. EXAMPLES In this section two examples are considered, namely those specified in equations (107 , (12). The exact inverses are given in equations (117 , (137 respectively. Of course the interest here is not in achieving the inversion, but rather as examples o f the method. We consider first (10). Here some interest attaches to the results since we obtain trigonometric - exponential function approximations to the unit step function which would seem to have important application in electrical network theory and the design of delay lines. In order to be able to use w = 0 for this case 1/p has been subtracted from f-(p) as explained above. In table I below are given values o f the At., ar for weight w = 0 and n = 1 (1) 7. Since for this case we know that K(0) = 1, we were able to compute the In(0 ) , and these are included in the table. In table 2 are given approximate exponents a r and amplitudes A r for a local minimum for the case Journal o f Computational and

n = 7 which was obtained during the course of the computations. We see that the value of I n obtained is larger than that for the case n = 6. Graphs o f the approximants gn(t) for n = 1 (17 7 are shown in fig. 1, 2, 3 where they are compared with the exact inverse H(t-1). We now turn to our second example as given in equation (127, the exact inverse being Jo(t). In contradistinction to the previous example we have here a continuous inverse and we may expect high accuracy for low values of n. However now the integral (77 for K(w) diverges for w = 0, and we have to choose a positive value for w. This o f course causes weighting o f our approximation in favor of the smaller values o f t, so that our accuracy decreases as we increase t. Table 3 shows exponents % and amplitudes A r for w = 1, 1.5, 2, 4 for n = 2 and n = 3. Taking advantage of the formula oo rt/2 f e-Wt [Jo (t)]2 dt = k__ f (1-k2sin20) -1/2 dO O

7/

O

where k = 2 (w 2 + 4) -1/2

(see [17, p. 250, no. 561.01] or [18, V. 183, no. 12"]) we are able to express K(w) in terms of the complete elliptic integral o f the first kind. Then for the above values o f w, K (w) was obtained with the aid of elliptic integral tables [Abrawowitz and Stegun, 19, Ch.17]. Thus we are able to calculate In(w ) also for this second example. The very small values of I n attained for low values o f n testify to the accuracy of the method for this example where f(t) is continuous. Table 4 shows a comparison between the exact inverse Jo(t) and g2(t), g3(t) for w = 1 and w = 4, for a few selected values o f t. 6. CONCLUSION The practical application o f a method of Laplace transform inversion by means of best rational function approximations has been demonstrated. It is suggested that the use of an appropriate "window" may facilitate accurate inversion for a given range o f t. It is hoped to develop further this concept in a later publication. ACKNOWLEDGMENTS The author is grateful to E. Natanel and D. Levin for drawing his attention to the minimization method o f Davidon as modified by Fletcher and Poweil, and to D. Levin for useful discussions and assistance in programming for the computer. The computations for this paper were carried out on the CDC 6600 computer at the Computation Centre o f Tel-Aviv University.

* Note error in this entry where LI2 (1/2 at)] 2 appears in place o f Llo(1/2at)] 2.

AppiiedMathematics, volume I, no 1, 1975.

19

TABLE 1 E x p o n e n t s ~r a n d amplitudes A r for f ( p ) = ( e - P - 1 ) / p having inverse f(t) = H ( t - 1 ) - l . •

Weight w _- 0.

K(0) _ o ~ ° [ f ( t ) ] 2 d t = 1.

n

~r

1

1.2564312086

2

1.5109396678 + 1.8261810223i 1.5109396678 - 1.8261810223i

3

4

5

6

7

t

A p p r o x i m a t e inverse is gn(t) = r_-I Are

Ar

In

- 1.4306637259

1.4486431360 + 4.1507410634i 1.4486431360 - 4.1507410634i 2.2466035635

" In = o~[gn(t)-f(t)]2d

.1854712448

.3533631161 - 1.1622826495 .3533631161 + 1.1622826490i .6880990747 + .6880990747 - 2.6064053780

.0860095029

.0646253880i .0646253880i .0529953681

1.3780882745 1.3780882745 2.3837990030 2.3837990030

+ + -

6.7422396356i 6.7422396356i 2.0875319714i 2.0875319714i

.3431240936 .3431240936 .0601336237 .0601336237

+ .2466989810i .2466989810i - 2.4126151283i + 2.4126151283i

1.3302795660 1.3302795660 2.2840537292 2.2840537292 2.8712498099

+ + -

9.4611323327i 9.4611323327i 4.5214951648i 4.5214951648i

.1247157114 .1247157114 1.3670835978 1.3670835978 -4.1527734398

2.9456997965 2.9456997965 2.1723107680 2.1723107680 1.2978874762 1.2978874762

+ + + -

2.2220563329i 2.2220563329i 7.1860430028i 7.1860430028i 12.2535299705i 12.2535299705i

.5295000919 - 3.7795276709i .5295000919 + 3.7795276709i .9435086350 + .2031908768i . 9 4 3 5 0 8 6 3 5 0 - .2031908768i .0101241432 + .2181536055i .0101241432 - .2181536055i

.0230233972

2.8206778021 2.8206778021 2.0898287291 2.0898287291 1.2749123805 1.2749123805 3.3205398634

+ 4.7169330622i 4.7169330622i + 9.9666952430i 9.9666952430i + 15.0950400613i - 15.0950400613i

1.9182576247 + 1.5939814919i 1.9182576247 - 1.5939814919i .5186153362.4188242928i .5186153362 + .4188242928i .0483875771 .1662336720i .0483875771 + .1662336720i - 5.9154864513

.0191865946

+ + -

.0374362771

.2633060391i .2633060391i .6840268300i .6840268300i .0286201202

TABLE 2 A p p r o x i m a t e e x p o n e n t s ctc and a m p l i t u d e s A r for a local m i n i m u m o f I n for the case n = 7 for f-Cp) = (e - p - 1)/p with w = 0 n

ar 2.3572 + 4 . 6 0 1 0 i 2.3572-. 4.6010i 1.3950 + 9 . 6 5 9 0 i 1.3950 9.6590i .7646 + 23.2368~ •7646 - 2 3 . 2 3 6 8 i 2.9452

Ar 1.4025 1.4025 .1882 .1882 .0469 .0469 - 4.4360

In

+ + +

.8523i .8523i .2530i .2530i .0133i .0133i

J o u r n a l o f C o m p u t a t i o n a l a n d Applied Mathematics, v o l u m e I, no 1, 1975.

.0259

TABLE 3 Exponents a r and amplitudes A r for f-(p) = (p2 + ! ) - 1 / 2 , with inverse f(t) = Jo(t). n Are_~r t " Approximate inverse is gn(t) = r _-~1 Case 1.

3

~r

Ar

In

.2476743906 + .8312538756i . 2 4 7 6 7 4 3 9 0 6 - .8312538756i

.4900495491 + . 2 1 8 3 3 8 2 9 5 2 i .4900495491 - . 2 1 8 3 3 8 2 9 5 2 i

2.4960x10 -4

.1126198382 + .9352523370i .1126198382 - .9352523370i .6180339888

.2972499381 + . 1 6 0 6 3 5 0 7 0 4 i .2972499381 - .1606350704i .4074549432

2.4605x10 -6

Case 2.

w = 1.5, K(w)

=

.50809968004853.

.2182458366 + .7905694150i . 2 1 8 2 4 5 8 3 6 6 - .7905694150i

.4959502923 + . 1 7 6 7 5 3 6 8 8 0 i .4959502923 - .1767536880i

3 . 3 0 6 5 x 1 0 -7

.1026012409 + .9140957959i .1026012409 - .9140957959i .5000000001

.3108985870 + .1351616683i .3108985870 - .1351616683i .3787150419

1.3338x10 -7

Case 3.

w = 2, K(w) = .41731342083704.

.1892071150 + .7653668647i .1892071150 - .7653668647i

.4981290830 + .1464425238i .4981290830 - .1464425238i

5 . 8 2 1 2 x 1 0 -6

.0906837538 + .9002271654i . 0 9 0 6 8 3 7 5 3 8 - .9002271654i .4142135626

.3183319443 + .1144043159i .3183319443 - . 1 1 4 4 0 4 3 1 5 9 i .3634974975

1.0871x10 -8

C ~ e 4.

3

oo e_Wt 2 K(w) = of [Jo(t)] dt.

w = 1, K(w) = . 6 4 2 6 3 7 6 8 1 7 6 7 1 8 .

n 2

oo - w t " 2 [ g n ( t ) - f ( t ) ] dt. of e

In

w = 4, K(w) = .23625158273244.

.1147425269 + .7265425280i .1147425269 - .7265425280i

.4998055590 + .0834626086i .4998055590 .0834626086i

3 . 5 7 1 4 x 1 0 -8

.0565823680 + .8777636960i .0565823680 - .8777636960i .2360681107

.3286218015 + .0671968188i . 3 2 8 6 2 1 8 0 1 5 - .0671968188i .3427618181

6.9400x10 -12

-

TABLE 4 Comparison between Jo(t) and the approximations gn(t) for n = 2 , 3 with w = 1 and w = 4. w=l t

w~4

Jo(t)

g2(t)

g3(t)

0 1.000000 0.5 .938470 1.0 " .765198 1.5 .511828 2.0 .223891 2.5 - .048384

.9801 .9480 .7675 .5007

1.00195 .93765 .76594 .51267 .22330 - .05004

.2103 -.0509

g2 (t)

g3 (t)

.99961 .93829 .76505 .51393 .22510 -.06087

1.000005 .938474 .765187 .511844 .223935 .048507

Journal o f Computational and Applied Mathematics, volume I, no 1, 1975.

21

REFERENCES 1. I.M. Longman, "Best rational function approximation for Laplace transform inversion", SIAM J. Math. Anal., _55 (1974), pp. 574-580.

11. ~ , "Approximate inversion of a class of Laplace transforms applicable to supersonic flow problems", Quart. J. Mech. Appl. Math., 1_7 (1964), pp. 91-103. 12. ~ , "The special functions and their approximations": Vol. 2, Academic Press, New York, 1969.

2. R. Fletcher and M. J. D. Powell, "A rapidly convergent descent method for minimization", Computer J., 6 (1963-4), pp. 163-8.

13. J. E. Akin and J. Counts, "The application of continued fractions to wave propagation in a semi-infinite elastic cylindrical membrane", J. Appl. Mech., 36 (1969), pp. 420-424.

3. I.M. Longrnan, "The application of rational approximations to the solution of problems in theoretical seismology", Bull. Seism. Soc. Am., 56 (1966), pp. 10451065.

14. - - - - , "On rational approximations to the inverse Laplace transform", SIAM J. Appl. Math., 1_~7(1969), pp. 1035-1040.

4. ~ "The numerical solution of theoretical seismic problems", Geophys. J. Roy, Astr. Soc., 13 (1967), pp. 103-116.

15. J. Counts and J. E. Akin, "The application of continued fractions to wave propagation problems in a viscoelastic rod", DEMVPI Research Rep. 1-1, Dept. of Eng. Mech., Virginia Polytechnic Institute, Blacksburg, 1968.

5. , "On the numerical inversion of the Laplace transform of a discontinuous original", J. Inst. Maths. Applics. 4 (1968), pp. 320-328. 6.

16. J. A. Murphy. "Certain rational function approximations to (1 q-x)-1/2,,, j . Inst. Maths. Applics.,__7 (1971), pp'. 138-150.

. , "Computation of theoretical seismograms", Geophys. J. Roy. Astr. Soc., 2_11(1970), pp. 295-305.

7. ~ , " N u m e r i c a l Laplace transform inversion of a function arising in viscoelasticity", J. Computational Phys., 100 (1972), pp. 224-231.

17. P. F. Byrd and M. D. Friedman, "Handbook of elliptic integrals for engineers and physicists", Springer-Verlag, Berlin, 1954.

8. - - , "Approximate Laplace transform inversion applied to a problem in electrical network theory", SIAM J. AppL Math., 23 (1972), pp. 439-445. 9.

18. A. Erd~lyi, W. Magnus, F. Oberhettinger and F. G. Tricomi, "Tables of integral transforms", Vol. 1, Bateman Manuscript Project, McGraw-Hill, New York, 1954.

, "On the generation of rational approximations for Laplace transform inversion with an application to viscoelasticy", SIAM J. Appl. Math., 2_44(1973), pp. 429-440.

t

19. M. Abramowitz and L A. Stegun (Eds.), "Handbook of mathematical functions with formulas , Graphs and Mathematical Tables", Dover, New York, 1968.

10. Y. L. Luke, "On the approximate inversion of some Laplace transforms", Proc. 4th U.S. Nat. Cong. on Appl. Mech., 1962, pp. 269-276.

n-Z 1.0

4~rl

I

,2

0.8 0.6 0.4 0.2

n-2

0 -0.2

-0.4~ 0

t

I 02

~

I 0.4

~

I 0.6

f

I 0.8

~

I 1.0

i

I L2

~ ~t

I L4

Fig. 1. G r a p h s o f the a p p r o x i m a n t s gn(t) for n = 1,

~

I L6

I

I L8

;

I I 2.0

I 2.2

i

I 2.4

I

I 2.6

;

I 2.8

2, 3 as c o m p a r e d with t h e e x a c t inverse H ( t - 1 ) .

J o u r n a l o f C o m p u t a t i o n a l and A p p l i e d Mathematics, volume I, n o 1, 1975.

22

n,,6

- 9n 1.0

i'

08 06

\

n'~S -

n-5

I

0.4 02

n-4

_o.o_

I

~'4

I o.2

"o

~

I 0.4

;'5

I

n.6n-4

I o.6

I

I o.e

~

I Lo

~

I

~

t2

I t4

s

I

I

t6

I

I

Le

I 2.0

I

I ~-2

I

I

J

2.4

I

I

2.6

I 28

~t

Fig. 2. Graphs of the approximants gn(t) for n = 4, 5, 6 as compared with the exact inverse H(t-1).

1.0

v

-"'--"~ "~-~" ~

(18 0.6 G4 0.2

0

----t

Fig. 3. Graph of the approximant g7(t) as compared with the exact inverse H(t-1). Journal of Computational and Applied Mathematics, volume I, no 1, 1975.

23