COMPUTER-AIDED NUMERICAL INVERSION OF LAPLACE TRANSFORM

This paper explores the technique for the computer aided numerical inversion of Laplace transform. The inversion technique is based on the properties of a family of three parameter exponential probability density functions. The only limitation in the technique is the word length of the computer being used. The Laplace transform has been used extensively in the frequency domain solution of linear, lumped time invariant networks but its application to the time domain has been limited, mainly because of the difficulty in finding the necessary poles and residues. The numerical inversion technique mentioned above does away with the poles and residues but uses precomputed numbers to find the time response. This technique is applicable to the solution of partially differentiable equations and certain classes of linear systems with time varying components.

The expected value of some function f of the random variable T is given by E[Z(t)] I(t)pm,,,(a, t)dt Obtaining binomial expansion of a portion of the density function gives (l e_at)n   Absolute error vs time FIGURE 3 Broken lines indicate bounds of errors.
changing order of the summation and integration, the result is 10-4 The integral can be identified as Laplace transform off(t) which is denoted by F(s)  where the function F(s) is defined as If the limit is taken as m and n approach infinity, the density function becomes very peaked and behave like a Dirac delta function located at the mode (peak) of the density and hence the density tends to sample the function f(t) at t=ln(n/m)/ m or f[ln((n + m)/m) ] m,n-oolim Elf(t)] Absolute error Vs ime FIGURE 7 Broken lines indicate bounds of errors for other N.
for convenience sake 'a' is assigned the value ln((m + n)/m) a= t A combination of (7) and (9)  If the S-domain function F(s) is known, then (11) may be employed to evaluate the corresponding time domain function for finite m and n, the estimate of the inverse transform f,,n(t) may be written as: Absolute error Vs lime Expression ( 13) is called the Gaver   and N must be an even integer.By recombination of similar terms in ( 14) we find that the averaged approximate inverse transform is given by fl(t) ln 2 ' VF ( ln 2 ) Vi is given by and k is computed using integer arithmatic.This expression is called the stehfest algorithm.The analytical form of F(s) is not limited to polynomial ratios.Here F(s) can include transcendental function or transport lag terms exp (-ts) which arise in the analysis of distributed parameter systems.
The method of "extrapolation to the limit" which Gaver used, leads to less accurate results for the same N, because not so many powers of 'n' cancel out.Moreover, with this method N must be a power of 2, so that in general one cannot make the best use of the available computer precision.Theoretically fn becomes more accurate for large values of N.
Practically, rounding errors worsen the results if N becomes too large because Vi with greater and greater absolute values occurs.For given F(s) and T, N at which the accuracy is maximal increases with the number of significant figures used.For fixed computer precision the optimal value of N is smaller i.e., maximum accuracy is greater and fn converges to f(t) faster.It was also found that with increasing N, the number of correct figures first increases linearly and then owing to rounding errors, decreases linearly.proportional to the nmber of digits the machine is working with.Evaluating an unknown function from its Laplace transform, one should compare the results for different N, to see whether the function is smooth enough (less oscillatory), what accuracy can be reached and what the optimum N is.One should also make sure that the unknown function f(t) has not any discontinuities, salient points, sharp peaks or rapid oscillations.AbsoLute error vs time FIGURE 16 Broken lines indicate bounds of errors for other values of N.
If the Gaver algorithm is employed, the coefficients with alternating signs create a problem in that large numbers of very nearly equal size must be subtracted from one another, thus implying the need for retention of a large number of significant digits.As the value of N is increased this becomes a word length problem for any digital processor.Similarly if the averaged approximate inverse algorithm is employed, there will be a value N beyond which a degradation of the inverse will occur due to finite digital processor word length.The algorithm works very well when oscillatory solutions are not anticipated, otherwise the asymptotic character of the inverse function soon dominates.An alternative method on similar lines is: Consider the Laplace transform inversion formula 1 (A) v(tl where s is the complex frequency variable and C is the arbitrary positive constant such that Resi < C where Si are the poles of V(s).
Using the substitution (A) changes to z = St (B) Approximate the function e z in (C) by a rational function.The Pade approximation was used. e,,,( where PN(Z) and QM(Z) are polynomials of order N,M respectively.As is well known, the Pade approximation formally equates a rational function to some terms of a series the coefficients Ci for < M + N are the coefficients of the Taylor expansion of ez.Therefore the Pade approximation N,M(z) has the first M+ N+ terms of its Taylor expansion equal to the Taylor expansion of ez.The remaining terms of both expansions differ.
It is not necessary to solve the system of equations arising from (E).A closed form exists V,M(Z)=(M+N)'+(M+N-1)'(N) z+(M+N-2) N z2 2 + + M ( N ) zl / (M+N)'-(M+N-1)'(M) z+(+N-) z 2 The function N,M can be written in a doubly infinite Pade table 0,0 1,0 2,0 0 All poles of (F) are simple and for M not differing considerably from N, all are in the right half plane.Although this has been used for the sake of simplicity, it can be shown that left half plane poles of (F) do not alter the results.
Inverting (D) into (C), the approximation l?(t) of Y(t) becomes, f C'+yoo (t) .c,_joo V(z/t) cV'M(Z)dz" Integral (G) can be evaluated by residue calculus by closing the path of integration along an infinite arc either to the right or to the left.In order that arc does not contribute to the integral M,N is taken such that the function has at least two more finite poles than zeros.Then F(z)dz -+-27rj (residue at poles inside the closed path). (I)   where the positive sign applies when the path is closed in the left half plane (counter-clockwise), whereas the negative one applies for the other case.Important properties of the method will be derived by alternatively closing the paths in both half planes.For N < M, M N,M(Z) Kilzzi (J) i=1 where zi are the poles of N,M(g) and ki are the corresponding residues.Closing the path of integration around the poles of N,M(Z) in the right half plane, M f"(t) =-1/t_aKi V(zi/t) (K) U. KUMAR decreases linearly.The optimum N is approximately proportional to the number of digits the machine is working with.
(b) The nature of the Gaver-Stehfest algorithms with large coefficients and alternating signs makes it necessary to subtract very large numbers of almost the same size, which creates word-length computational problems.Particular choices of N were made because of the limitation on N due to the finite word length of the digital processor used.With 32 bits, we get a 7 decimal digit resolution.In cases where N 18, the inverse function will demonstrate considerable scatter due to truncation error in the processor.
(c) It is clear that functions with oscillatory inverses present difficulty to the method.
(d) There is considerable improvement due to the multiplication algorithm used.
(e) Transforms for lightly damped sinusoids and BesSel functions require large mainframe computers with long word-length.
The accuracy of the method was considered from the theoretical and empirical points of view.It was found that the accuracy is high for small time.Examples were used to illustrate the idea.
FIGUREThree parameter exponential probability density functions.

FIGURE 2
FIGURE 2 Broken lines indicate bounds of errors for other values of N.

FIGURE 5
FIGURE 5 Broken lines indicate bounds of error for other values of N.

FIGURE 6
FIGURE 6 Broken line indicates bounds of errors for other values of N.

FIGURE 11
FIGURE 11 FIGURE 12 Broken line indicates bounds of errors for other values of N. FIGURE13

FIGURE 14
FIGURE14 Broken lines indicate bounds of errors for other values of N.

FIGURE 17
FIGURE 17 Any sequence of these functions converges to e z for any z as lim (M+ N) --o.