Journal of Computational and Applied Mathematics 10 (1984) 113-132 113 North-Holland

**Algorithm 27 **

**A method for the numerical inversion ** **transforms **

**of Laplace **

G. H O N I G *

*Departamento de Inform4tiea, Pontifica Universidade Catrlica, Rio de Janeiro, Brasil *
**U. **H I R D E S

*lnstitut fi~r FestkOrperforsehung, Kernforschungsanlage Ji~lich, Germany, Fed. Rep. *

Received 19 February 1980 Revised 4 July 1983

*Abstract: *In this paper a numerical inversion method for Laplace transforms, based on a Fourier series expansion developed by
D u r b i n [5], is presented. The disadvantage of the inversion methods of that type, the encountered dependence of discretization and
truncation error on the free parameters, is removed by the simultaneous application of a procedure for the reduction of the
discretization error, a method for accelerating the convergence of the Fourier series and a procedure that computes approximately the
' b e s t ' choice of the free parameters. Suitable for a given problem, the inversion method allows the adequate application of these
procedures. Therefore, in a big range of applications a high accuracy can be achieved with only a few function evaluations of the
Laplace transform. The inversion method is implemented as a F O R T R A N subroutine.

*Keywords: *Numerical integration, Laplace transforms, numerical Laplace inversion.

**1. Introduction **

The significance of numerical Laplace inversion is obvious from the big range of applications. Well known in engineering, Laplace transformation methods are also used in order to solve differential and integral equations and to assist when other numerical methods are applied (see [7,10]).

A number of numerical inversion methods has been developed during the last few years. In what follows we confine ourselves to the methods using Fourier series approximations. Many tests (see [1,2,4,5,8] have demonstrated the simplicity and the accuracy of these methods.

Fourier series were first used by Dubner and Abate [4] in 1968 for the numerical inversion of Laplace transforms. Durbin [5] improved the method in 1973.

Other authors, Simon et al. [8] in 1972, Veillon [9] in 1974 and Crump [2] in 1976, used different acceleration methods in order to speed up the convergence of the Fourier series derived in [4,5]. Some of these methods at times achieve a considerable reduction of the truncation error. However, they fail in other cases and their efficiency heavily depends on the choice of the parameters of the methods of Durbin or Dubner and Abate, and the choice of these parameters was somewhat arbitrary.

The biggest disadvantage of the above mentioned methods is the dependence of the discretization and

* Present address: Deuzer Maschinenfabrik G M B H , Postfach 3165, D5902 Netphen 3, Germany, Fed. Rep.

0377-0427/84/$3.00 © 1984, Elsevier Science Publishers B.V. (North-Holland)

114 *G. Honi~ U. Hirdes / Laplace transform *

truncation errors on the free parameters: by a suitable choice of these parameters the discretization error becomes arbitrarily small, but at the same time the truncation error grows to infinity and vice versa. A first step in surmounting this problem was taken in [1], (see Section 3), where the so-called ' K o r r e k t u r ' - m e t h o d was presented. It allows a remarkable reduction of the discretization error without increasing the truncation error.

Nevertheless the accuracy of the ' K o r r e k t u r ' - m e t h o d also depends on a ' g o o d ' choice of the free parameters. The procedure introduced in Section 5 gets closer to the solution of the problem. It allows the a p p r o x i m a t e computation of optimal parameters for all above mentioned methods (for the definition of ' o p t i m a l ' see Section 5).

Section 4 describes a new method for the acceleration of convergence of the Fourier series. It is applicable if the infinite series for the approximation of the inversion integral does not alternate (see Section 4). In this case all other tested acceleration methods fail.

Finally, in Sections 2 and 4, those of the above mentioned inversion and acceleration methods that were found to be best by numerical tests as well as by the derived error estimates are summarized.

A new algorithm for the numerical Laplace inversion, taking into consideration the new, above mentioned, procedures and implemented as a F O R T R A N subroutine 1 is described in Section 6.

This subroutine can be used as a ' b l a c k box': the only input parameters are the t-values for which

*f(t) = L-I[F(s)] *

shall be computed and of course the Laplace transform *F(s). *

On the other h a n d - - f o r a
concrete p r o b l e m - - t h e user can make an optimal choice of all free parameters in order to have the most
efficient combination of the algorithms contained in the subroutine.
The efficiency of the inversion method is shown by the examples in Section 7.

**2. The method of Durbin **

The Laplace transform of a real function f : R ~ R with f ( t ) = 0 for t < 0 and its inversion formula are defined as

*F(s) = L[f(t)] *

= f 0 °~ e - ' t f ( t ) d t , *(1) *

1

### f~+i~ *eS,F(s)ds, *

*f( t)= L-l[ F(s)] = T~i o-i~ * (2)

with s = v + iw; v, w ~ R.

v c R is arbitrary, but greater than the real parts of all the singularities of

*F(s). *

The integrals in (1) and
(2) exist for R e ( s ) > a ~ R if
(a) f is locally integrable,

(b) there exist a t o >/0 and k, a ~ R, such that If(t)[ ~< k e ~t for all t >/t 0, (3) (c) for all t ~ (0, ~ ) there is a neighbourhood in which f is of bounded variation.

In the following we always assume that f fulfils the conditions (3) and in addition that there are no singularities of

*F(s) *

to the right of the origin. Therefore (1) and (2) are defined for all y > 0. The possibility
to choose v > 0 arbitrarily, is the basis of the methods of D u b n e r - A b a t e and Durbin. The latter is now
described.
Using the inversion formula (see [5])

*f(t) =--~-e°'f~[Re{F(S)~o *

} cos(wt) - I m ( F ( s ) } *sin(wt)]dw, *

(4)
1 See Appendix A, subroutine L A P I N (LAPlace INversion).

*G. Honig, U. Hirdes / Laplace transform * 1 1 5

which is equivalent to (2), and a Fourier series expansion of h ( t ) = e "'

*f(t) *

in the interval [0, 2T], Durbin
derived the approximation formula
f ( t ) = - ~ - - ½ R e { F ( v ) } + Y'~ Re F v + , - - ~ -

k = 0

### x c o s - ~ - t - ~ Im f v + i sin ~--t *- r l ( v , t , r ) , * (5)

k = 0

for 0 < t < 2T.

*Fl(v, t, T) *

is the discretization error, given by
*FI(v, t, T) = £ e-2"krf(2kT+ ,). *

^{(6) }

*k = l *

Since there are no singularities of

*F(s) *

in the right halfplane we have (see [3, Bd. 1]) a c > 0, rn >~ 0 and a
t o >/0, such that
[ / ( t ) ] ~<

*ct m *

for allot > t 0. (7)
From (6) and (7) the following estimates for the discretization error can be deduced (see [5]):

(a) m = 0,

[ r a ( v , t, T)[ ~ < - - , C (8)

e 2 r e - - 1

(b) m > 0,

*iFl(v,t,T)l<K(2T)me_2,,r( a, *

+ " - + - - - - *a m + , ) *

, (9)
(2 v) m+'
where K, a I . . . a,,+l ~ R.

These estimates show that the discretization error can be made arbitrarily small by choosing v sufficiently large.

As the infinite series in (5) can only be summed up to a finite number N of terms, there also occurs a truncation error given by

**ImI ( + **

k=N+l ~ - - t - - i sin ~ - - t .(10)

Hence the approximate value for

*f(t) *

^{is }

= e " ' [ - 1 R e { F ( v ) } + ~ { R e { F ( o + i ~ - ~ ) } ~-~t

*fN(t) ~- *

^{k=0 }

^{COS }

*-Im{F(v+i~)}sin ~-~t}]. *

^{(11) }

**3. The reduction of the discretization error by the 'Korrektur'-method **

It is obvious from (8) and (9) that the discretization error can be made arbitrarily small if the product vT is sufficiently large.

Unfortunately, the truncation error (10) may diverge for large values of

*vT. *

*116 * *G. Honig, U. Hirdes / Laplace transform *

T h e ' K o r r e k t u r ' - m e t h o d allowes a reduction orf the discretization error without enlarging the truncation error. With (11), equation (5) can be written in the f o r m

*f ( t ) *

= f ~ ( t ) - F l ( v , *t, T ) . *

(12)
T h e ' K o r r e k t u r ' - m e t h o d uses the a p p r o x i m a t i o n formula

*f ( t ) = f ~ o ( t ) - e 2"vf~(ZT+ t ) - r Z ( v , * *t, T ) . *

(13)
As stated in the theorem below, the discretization error

*F2(v, t, T ) *

is m u c h smaller than *FI(v, t, T). *

Taking into consideration the truncation error (10), we find that the a p p r o x i m a t e value for f ( t ) , using the ' K o r r e k t u r ' - m e t h o d (13), is given by

*fuK ( t ) = fN( t ) -- e - Z~TfNo(ZT + t ). *

(14)
It should be m e n t i o n e d that the truncation error of the ' K o r r e k t u r ' - t e r m e - 2 " T f ~ ( 2 T

*+ t), *

is m u c h smaller
than F A ( N , v, t, T ) if N = N 0. W e can therefore choose N o < N (see Section 6), which means that only a
few additional function evaluations of *F ( s ) *

are necessary to achieve a considerable reduction of the
discretization error. This error reduction is pointed out in the following theorem.
T h e o r e m 3.1.

*Suppose f is a real function with f ( t ) = 0 for t < 0 that possesses the properties *

(3) *and * *F ( s ) = L [ f ( t ) ] its Laplace transform that has no singularities to the right of the origin, and suppose the *

*"Korrektur '-formula *

(13) *is used for the numerical inversion of F( s ), then *

(a) I F 2 ( v , t, T ) l <

*2e * *i f m = O , *

(15)
eZ~'r(e 2 v T - 1)

(b)

*, F 2 ( v , t , T ) l ~ 3 " e - 2 ~ T { K ( 2 T ) " e - 2 V r ( *

~ - ~ + . . - + *a' * *a"+, * *) ) *

(16)
( 2 v T ) " + '
*if m > 0 and *

( m ! ) / 2 " - 1 -%< *2vT. *

F o r the definition of c, m, K, a 1 .. . . . a,,+~ see equations (7), (8) and (9). F o r the p r o o f see [1] or [6].

Remarks. A c o m p a r i s o n with the m e t h o d of D u r b i n (see (8) and (9)) shows a reduction of the error b o u n d b y the factor 2e-2v T for m = 0 and 3" e - 2 v r for m > 0. T h e condition

*m ! / 2 " - 1 ~ 2 v T *

might be difficult
to fulfill for large m, and hence the application of the ' K o r r e k t u r ' - m e t h o d is not r e c o m m e n d e d in such
cases.
As m e n t i o n e d before, an adequate reduction of the total error can be obtained by the ' K o r r e k t u r ' - m e t h o d only if (for fixed N and T ) the p a r a m e t e r v is suitable. This is illustrated in Fig. 1. It shows the error curves

0 0015 o

00010 o

o

0 0005

00000 0

Fig. l.

### IF21+IFAI~t~

1

**~ **

I F I I * I F A I
I F 2 1 ~ Fll
~ , _

vl T Vo T 3 4 ~T

*G. Honig, U. Hirdes / Laplace transform * 117

of the method of Durbin ( F 1 and FA) and of the ' K o r r e k t u r ' - m e t h o d ( F 2 and FA). A successful application of the ' K o r r e k t u r ' - m e t h o d requires a v < %. For v > v 0 the truncation error dominates and the reduced discretization error F2 does not lead to a noticeable reduction of the total error.

The opposite holds if the methods for the acceleration of convergence, that we describe in the following part, are applied. Acceleration of convergence is only sensible if v > v 0 or, which is the same, if the truncation error dominates.

However, a simultaneous application of an acceleration method and the ' K o r r e k t u r ' - m e t h o d (as realized in the subroutine L A P I N ) is recommended, if the parameters are optimally chosen by the procedures introduced in Section 5.

**4. Acceleration of convergence **

In this section three acceleration methods used in the subroutine L A P I N (see Appendix A) are briefly described: the e-algorithm, the m i n i m u m - m a x i m u m method and a method based on curve fitting. The latter is new, whereas the e-algorithm was already used in [2,9], the m i n i m u m - m a x i m u m method in [1] in order to speed up the convergence of the series approximating the Laplace inversion integral. We have also tested other acceleration methods such as the Euler transformation, or Aitken's extrapolation p r o c e d u r e (see [8]). But these turned out to be less efficient than the above mentioned methods.

We can consider f h ( t ) (see (11)) for fixed t as a discrete function of N and define:

**Definition. **

*fN(t) *

as a function of N is monotonous if
*IfN(t) *

- - f ~ (t)] >~ *IfM(t) *

- - f ~ ( t ) [
for all N, M with N ~< M.
For a n o n - m o n o t o n o u s

*fN(t) *

the e-algorithm (EPAL) and the m i n i m u m - m a x i m u m method ( M I N I -
M A X ) in general significantly increase the rate of convergence. However, they f a i l - - a s do the Euler
transformation and Aitken's extrapolation p r o c e d u r e - - f o r a monotonous *fN(t). *

The method based on
curve fitting (CFM) leads in this case to a considerable improvement in the results. With
c k , = - ~ - Re F v + i cos T t - I m we can replace (11) by

N

*fN(t)=½Co + E Ck" *

k = l

The

*e-algorithm *

applied to (17) is defined as follows:
Let N = 2 q + 1,

*q~N, *

*S m := ~ * *Ck *
*k = l *

and

v+l--~-- s i n - f - t , k = 0 , 1 , 2 . . .

(17)

~°~+1) - e(pm') e~om' '= O, ( " ~ ' - (18)

(,,1) .__ ( r e + l ) + *1 / ( e p *

*Ep + 1 "-- Ep _ 1 * *~ * *E1 * *"-- Sm ' *

then the sequence e~ 1), E~ 1), or1) ~5 , ' ' ' , ~ 2 q + l oo) = e~ ), converges to

*foo(t)-Co/2. *

For a n o n - m o n o t o n e *fN(t) *

in
general it converges faster than the sequence of the partial sums *s,,, *

m = 1, 2 . . . . of the untransformed
series.
118 *G. Honig, U. Hirdes / Laplace transform *

The

*minimum-maximum method *

also increases the rate of convergence only in the case of a non-monot-
onous *fN (t). *

The method works as follows. Having found three neighbouring stationary values *of fN (t) *

as a
function of N, say a maximum at N = N1 and N = N3, and a minimum at N = N2, an interpolating
function for the maximum values at N = N1 and N = N3 is constructed. The mean value of the minimum
~it N = N2 and the function value of the interpolating function at N = N2 yields the new approximation for

*f~(t). *

For a more detailed description of M I N I M A X see [1].
The

*method based on curve fitting *

(CFM), applicable only for monotonous *fN (t), *

consists in fitting the
parameters of any function that has a horizontal asymptote *"/A(X) *

-- ~', by demanding that this function is
an interpolating function for the points *(N, fN(t)), 0 <~ N O <~ N <~ N 1. *

The function value of the asymptote ~,
is the desired approximation for *f ( t ) . *

With the simple rational function
r ( x ) : = 7 + - - + Y, (19)

for example, we achieved high accuracy already for small N t.

The C F M fills an important gap: now it is also possible to speed up the convergence of the Fourier series (5) in case of monotonous

*f , ( t ) 2. *

The subroutine L A P I N (see Appendix A) automatically chooses
between C F M and EPAL ( M I N I M A X ) depending on whether/N *(t) *

is a monotonous function of N. Tests
have shown that for non-monotonous *fN (t) *

EPAL mostly is superior to M I N I M A X in accuracy. However,
we found examples where the e-algorithm falsifies the results, but the m i n i m u m - m a x i m u m method did not.
In these cases, the subroutine L A P I N chooses by itself the M I N I M A X method to accelerate the convergence of the series (for more details, see Appendix A, significance of the parameter H). Furthermore, the e-algorithm is not applicable for arbitrarily large N in (17), because of overflows occuring in (18) for larger N.

**5. The choice of optimal parameters **

We already mentioned that a good choice of the free parameters N and C O N =

*vT *

is not only important
for the accuracy of the results but also for the application of the 'Korrektur'-method and the methods for
the acceleration of convergence (see Fig. 1). These methods do not improve the results if the parameters are
chosen badly.
Two methods are now presented which approximately determine the 'optimal' v for fixed N and 7". The main difference between the two methods is the definition of 'optimal parameters'.

**Definition **A. For fixed N and T the parameter v of the method of Durbin (see (12)), with or without the
application of the 'Korrektur'-method or a method for the acceleration of convergence, is optimal if the
absolute values of discretization and truncation error are equal.

**Definition **B. For fixed N and T the parameter v of the above mentioned methods (see Definition A) is
optimal if the sum of the absolute values of discretization and truncation error is minimal (see Fig. 1, where
the optimal v is given by v 0 and vl).

**Method **A. Let

*fN(t), fNK(t) *

be the approximate values for *f ( t ) *

computed by the method of Durbin (12)
and the 'Korrektur'-method (13) respectively. The bar indicates that one of the methods for the
acceleration of convergence may have been applied.
Let

*R ( N , v, t, T) *

be the expression in the brackets on the right side of (10). Neglecting the dependence
of v in R we can write for fixed *t, T: R ( N , v, t, T ) - - R ( N ) . & *

where 3 E [ - 1 , + 1 ] indicates the
application of an acceleration method (6 = 1 if none of the acceleration methods is used). The truncation of
*f N ( t ) orfNK(t ) *

is therefore given by
eV' (20)

F A ( N ,

*v, t, T) =--~-R( N ) 8 . *

2 For instance, fN(to) often is a monotonous function of N, if f(t) is a step function with the discontinuity in t 0.

*G. Honig, U. Hirdes / Laplace transform *

With **(5), (6) and (20) **we find

*iN(t) * *m y ( t ) - ~-~ R( N)3 + O(e-2"T). *

Let vl, 02 be large and v I 4:v23, then

*f~v( t) --]2u( t ) ~ R( N )3. (e °2'- eV't)/T *

o r

• • • t •

*R ( N ) 3 - - T "~'*'t" *

*eVz t _ eO,t *

Similarly we get for the 'Korrektur'-method

*R ( N ) 3 --- T f/~K(t) *

*~ f ~ K ~ t ~*

*cO2 t _ colt *

The discretization error (6) of the method of Durbin can be written as (see Section 6) F l ( o , t, T) = e - ~ ° 7 ( 2 T + t) + O(e-4vT).

From (13) we find for the discretization error of the 'Korrektur'-method

*F 2 ( o , t , T ) = ~ e - 2 V r k ( f ( 2 k T + t ) - f ( 2 T ( 3 k - 2 ) + t ) } . *

k = 2

Hence

119

(21)

(22)

(23a)

(23b)

(24a)

*F2(v, t, T) = * *e 4vT( *

f ( 4 T + t) - f ( 8 T + t)} + O(e-6Vr). (24b)
Using Definition A, the equations (23) and (24) easily lead to the following equations for the optimal
parameter v = VoAvr:
1 In

*R ( N ) 3 *

I (method of Durbin)
vAvr-- 2 T + t

*TfN(2T+t ) * **I **

^{(25a) }

and

**1 ** **In ** **R ( N ) 6 ** (' Korrektur'-method). (25b)

V°ar'r = 4 T + t T{ f N ( 4 T + t ) - ) ~ ( 8 T + t)}

Because of Definition A an upper bound for the total absolute error I f ( t ) - - f u ( t ) l or I f ( t ) - - f u K ( t ) l for v = vAvr is approximately given by

/)ApT • l

T E R R -= 2 ~ I R ( N ) 3 I. (26)

Method B. We now use Definition B and assume that

*R(N, v, t, T) *

is no longer a constant function of v to
get a more accurate approximation formula for Vop T.
The dependence of v is assumed to be as follows

*(R(N, v, t, T ) = R(N, v)8 *

because t, T are fixed,
v~ > 0 arbitrary):
*R ( N , v ) 3 = - R ( N , va)3. *

Re F v + i c o s - ~ - - - t - I m F v + l - - T- sin--~-t
*/[Re(F(vt+i-~-)}cos--~-tN¢r *

_ i m { F ( v l + i _ _ ~ ) ) s i n L "n ] 1 t ] - i (27)
3 We found that for example v~ ,= 20, v 2 ,= v I - 2 is a good choice.

120 *G. Honig, U. Hirdes / Laplace transform *

The acceleration factor 8 can be computed very easily from (23a) without additional function evaluations:

8 = (28)

*f ~ ( t ) --UZN(t ) *

Let AR(v~, v, N ) be the expression in brackets on the right side of (27) and

*R(N, *

vl)8 ,= R ( N ) & where
*R ( N ) 8 *

is c o m p u t e d by Method A. then (27) is expressible in the form
*R( N, v)8 = R( N ) 8 . AR(vl, *

v, N ) . (29)
Definition B implies (for the method of Durbin);

**e. ** **] **

_{= 0. }

_{( 3 0 ) }*[R(N, v ) 8 1 ~ - + l F l ( v , t, *

T ) I , ... ,%T
Approximating the derivative of

*R(N, v) *

by finite differences, we find the following iteration procedure to
solve equation (30):
*R ( N , v) * *--- *

(31a)
u ( i - 1) __ *U~ i-1) *

v (°) = oAvr, V~ °, = v I , (318)

*v~ i) = v (i-t), *

i = 1 . . . s - 1. (31c)
*v(i) *

1 In *~v]R(N'v)(°18 *

+ [ R ( N ' v *(i *

1))16.t
= , i = 1, 2 . . . s. (31d)

2 T + t 2 T 2 l f u ( 2 T + t) I

In (31b) v~zr and v 1 are defined by Method A (see equation (25) and (22) respectively). Equation (31d) follows from (30) substituting (31a) for the derivative of

*R(N, v) *

and (24a) for F l ( v , t, T).
Hence the optimal p a r a m e t e r v = VoBvr for fixed N, T, t is given by

**VoB r = v ( ' . ** **( 3 2 ) **

An analogous result holds for the ' K o r r e k t u r ' - m e t h o d : Replacing F1 by F2 in (30), we only have to substitute (31d) by

v(~)= 1 In , i = 1, 2 . . . s. (31e)

4 T +

*t * *4TZ[?N(4T+ t) * *--fN(8T+ *

t)l
It is sufficient to choose s = 1 or s = 2, as numerical tests have shown.
Again it is possible to compute approximately an upper bound for the total error, if the parameter v = v~v r is used for the computation of

*fN (t) a. *

We find
*v~°"t * *2°"°~rlfN(2T+ *

t)l. (33)
TERR----e ~ I R ( N , voBvr)l 8 + e -

Remarks. The simplifications m a d e in order to derive equations (25) and (31) do in general not cause any difficulties. Besides a few exceptions, which are discussed later, the optimal parameters VoApT and VoBvr are a good approximation for the true optimal parameters.

If the e-algorithm is used in order to accelerate the convergence of the series (5), it is not possible to c o m p u t e the acceleration factor 8 for large N with the desired accuracy (the e-algorithm breaks down after a certain N O < N depending on v, t and T; N O is not known in advance; see also the final remarks of Section

4 A similar equation holds

*forfNK(t). *

*G. Honig, U. Hirdes / Laplace transform * 121

4). Therefore the m i n i m u m - m a x i m u m method allows a more accurate determination of the optimal parameters by Method A or B, and hence of the total error (26) or (33).

**6. Remarks about the use of the algorithm **

T1, TN I M A N

First some input parameters of the subroutine L A P I N are described (see also Appendix A):

Lower and upper bound of the interval in which

*f ( t ) *

shall be approximated.
The choice of this parameter decides whether the free parameters of the subroutine are placed automatically or not:

I M A N = 0 automatical,

I M A N = 1 manual choice of parameters.

I L A P I N If I L A P I N = 1, the approximate value

*fN(t) *

(see (11)) is computed with T = t, if I L A P I N = 2
with T = TN. For T = t the computation of sine and cosine terms and of the imaginary parts of
*F(s) *

in (11) is cancelled, on the other hand R e { F ( s ) ) becomes dependent on t. Hence
I L A P I N = 1 is recommended if the inverse is computed only for a few t-values.
I K O N V 1KONV = 1 implies the use of the acceleration method M I N I M A X , I K O N V = 2 the accelera- tion of convergence with EPAL. In both cases the subroutine automatically chooses the curve fitting method (CFM) for

*monotonous fu(t ). *

I C O N If this parameter is zero ( I C O N = 0) the parameter v is not optimally chosen by Method A or B (see Section 5). For I C O N = 1 the subroutine L A P I N chooses an optimal v for t =

*tls/2 J *

^{(see }

Appendix A) and uses this v for the computation of f ( t ) in the whole interval.

If I C O N = 2 an optimal v is computed for all t k ~ (T1, TN) (see Appendix A for the definition of t,).

I K O R I K O R = 1 leads to the application of the 'Korrektur'-method, I K O R = 0 to the approximation o f f ( t ) without the 'Korrektur'-method.

Table 1 shows the possible combinations of the above described algorithms offered by the subroutine L A P I N . Method A for an optimal choice of the free parameters is used if I L A P I N = 2, Method B if I L A P I N = 1.

The application of the 'Korrektur'-method is not recommended first, in the case that the absolute value of the 'Korrektur'-term

*f ( 2 T + t) *

is small or equal to zero: ] f ( 2 T + t)l<< 1 (no noticeable improvement of
the accuracy occurs), and second, if *f ( t ) *

is a rapidly increasing function as t ~ oo (see the remarks to
Theorem 3.1).
The same holds for the application of the methods for an optimal choice of the free parameters if the auxiliary quantities (these are f ( 2 T + t) in (24a) and f ( 4 T + t ) f ( 8 T + t) in (24b)) are equal to zero.

Although the denominators in (25) do not vanish (because of the discretization and truncation errors), the computed optimal parameters might be misleading. In these cases, it is very helpful to have a global conception about the Laplace inverse

*f(t). *

If not already known, this can be obtained easily by the
subroutine L A P I N choosing either I M A N = 0 or I M A N = 1, I L A P I N = 2, I C O N = 0, I K O R = 0 and
I K O N V = 1. The latter combination of the parameters coincides with the method of Durbin, at which the
M I N I M A X - m e t h o d is used to speed up the convergence of the Fourier series.
Let NS1 be the number of function evaluations of

*F(s) *

used to approximate *f ( t ) *

and NS2 the number
of function evaluations used to approximate the 'Korrektur'-term *f ( 2 T + t). *

For given NS1 we found that
Table 1

I L A P I N I C O N I K O R I K O N V

I M A N = 1 1 0 / 1 / 2 0 / 1 1 / 2

I M A N = 1 2 0 / l 0 / 1 1 / 2

I M A N = 0 l 1 0 2

122

T a b l e 2

*G. Honig, U. Hirdes / Laplace transform *

T t o t o 2 t I + t o 4 t l + t o 8t 1 + t o

C O N 2 0 18 5 5 5

I L A P I N = 1 I L A P I N = 2 I L A P I N = 1

I C O N = 1 ! C O N = 1 I C O N = 2

*t o * *T I + ( T N - T 1 ) [ N ] / ( N + I ) * *T I + ( T N - T 1 ) [ N ] / ( N + I ) * *T I + k ( T N - T 1 ) / ( N + I ) , k = I ( 1 ) N * *a *
*t I * *T I + ( T N - T 1 ) [ 2 ] / ( N + I * *) * T N *T I + k ( T N - T 1 ) / ( N + I ) , k = I ( 1 ) N " *

a N = n u m b e r o f t - v a l u e s f o r w h i c h *f(t) *shall b e a p p r o x i m a t e d .

N S 1 / 2 ~< NS2 ~< NS1 if NS1 ~< 100 and N S 1 / 1 0 ~< NS2 ~< N S 1 / 2 if NS1 >~ 100 is a good choice for NS2.

As stated above, occuring overflows in (18) often make it impossible to transform a given n u m b e r NS1
of terms of the series (17) by the e-algorithm. If it breaks down after **N < **NS1 iterations the accuracy is
diminished with that the optimal parameters and the total error (33) are calculated. The M I N I M A X - a l g o -
rithm does not cause such problems.

Finally, we would like to mention that it can be necessary to know the s-values in advance for that *F ( s ) *
has to be evaluated (for instance, if *F ( s ) * isa not analytically known solution of a differential or integral
equation). These s-values are defined by T and C O N *= v . T ( s = v + i k v / T , * k = 0(1)NSUM). In Table 2,
the values of C O N and T, for which the 'auxiliary quantities' ( f N s l ( T ) ) must be computed, are listed.

If I C O N = 0 no auxiliary quantities are needed. If I C O N > 0, I K O R = 0 requires the computation of fNsl(t0) and fNm(2ta + to), I K O R = 1 the computation of fNsl(t0), fNSl(4tl + t 0) and f y m ( 8 t l + to) with

the above given values of CON, t o and q-

**7. Numerical examples **

**From the following examples one can get some idea of the accuracy of the subroutine LAPIN and the **
possibilities it offers to find an ' o p t i m a l ' solution for a given problem by a suitable choice of the different

T a b l e 3

F ( t ) = *t.sin(t)/2, *^{F ( s = }*s(s 2 + *^{1) - 2 }

M e t h o d D u r b i n L A P I N

C o d e ( v . T = 5) 2 - 1 - 0 - 1 I M A N = 0 1 - 2 - 1 - 1

N S 1 + N S 2 2 0 0 0 30 60 160 + 140

t R e a l a b s o l u t e e r r o r b y (33)

1 0.4 D - 0 2 0.2 D - 0 2 0.1 D - 1 1 0.351 D - 1 0 0 . 2 7 3 D - 1 0

3 0.6 D - 0 1 0.1 D - 0 1 0.1 D - 1 1 0.211 D - 0 9 0 . 2 3 6 D - 0 9

5 0.1 D - 0 1 0.2 D - 0 2 0.1 D - 1 2 0.573 D - 0 9 0 . 6 4 8 D - 0 9

7 0.8 D - 0 1 0.9 D - 0 2 0.2 D - 0 9 0 . 7 8 6 D - 0 9 0 . 9 7 9 D - 0 9

9 0.2 D - 0 1 0.9 D - 0 3 0.1 D - 1 0 0 . 7 7 4 D - 0 9 0.191 D - 0 8

11 0.1 D - 0 1 0.3 D - 0 3 0.2 D - 1 0 0 . 5 9 0 D - 0 9 0.273 D - 0 8

13 0.6 D - ) 1 0.2 D - 01 0.9 D - 1 0 0 . 1 0 7 D - 0 8 0 . 3 7 2 D - 08

15 0.5 D - 0 2 0.1 D - 0 1 0.2 D - 0 9 0.175 D - 0 8 0 . 3 5 6 D - 0 8

17 0.7 D - 0 2 0.2 D - 0 1 0.1 D - 0 9 0.225 D - 0 8 0 . 5 6 9 D - 0 8

19 0.1 D - 0 0 0.5 D - 0 2 0.4 D - 1 0 0.321 D - 0 9 0.633 D - 0 8

C P U - t i m e sec ( f o r

t = 1(1)20) 1.21 0.04 0 . 0 4 2 1.14

A l l c a l c u l a t i o n s w e r e p e r f o r m e d in d o u b l e p r e c i s i o n o n t h e I B M 3 7 0 / 1 6 8 c o m p u t e r o f K F A Jiilich.

*G. Honig, U. Hirdes / Laplace transform * 123

algorithms described above. The parameters used in the following tables are defined in Section 6. The code numbers refer to the parameters I L A P I N - I C O N - I K O R - I K O N V ( 1 - 2 - 0 - 1 , for example, means: I L A P I N

= 1, I C O N = 2, I K O R = 0 and I K O N V = 1). I M A N = 0 coincides with 1 - 1 - 0 - 2 .

Table 3 shows the errors occuring if

*F(s)= *

s(s 2 + 1)-2 is inverted. Three different parameter combina-
tions of L A P I N are compared with the method of Durbin. With only a few function evaluations of *F(s) *

and little CPU-time, L A P I N computes the Laplace inverse with a considerable accuracy. For instance, the first two columns show that the method of Durbin needs about 60 times more function evaluations and 30 times more CPU-time than L A P I N to get a comparable accuracy. Column 3 s h o w s - - a n d this was confirmed by many other e x a m p l e s - - t h a t using the subroutine L A P I N as a 'black box' ( I M A N = 0) yields excellent results. From a comparison between columns 5 and 6 it is obvious that the absolute error computed by (33) is a good approximation for the real error (such an error estimation is only possible if I C O N = 2; it is most accurate for I K O N V = 1).

As a second example, the Laplace transform of the step function 0 t < l O ,

*U(t- *

10),= 0.5 t = 10,
**1**t > 1 0 ,

was inverted (see Table 4). Comparing the results of Durbin with those obtained from L A P I N we come to the same conclusion as in the example above. The high accuracy at t = 10 is possible only if the curve-fitting method is used to speed up the convergence of the series. EPAL and M I N I M A X fail because fN(10) is m o n o t o n o u s in N. But also at points very close to the discontinuity at t = 10 a high accuracy is

obtained by L A P I N as demonstrated in the last column of Table 4.

From Table 5 one gets some idea of the effect of the 'Korrektur'-method. It can be seen that (with the same number of function evaluations) the 'Korrektur'-method (right part of Table 5) is clearly more accurate. This is not only caused by the reduction of the discretization error. The optimal parameter r o y r • T is smaller if the ' Korrektur'-method is applied. Hence the truncation error is reduced too.

As a final example, we show in Table 6 the dependence of the optimal parameter Vov r on the number NS1 of function evaluations (Oov r also depends on t). It is obvious that primarily the possibility to

T a b l e 4

*f ( t ) = U(t - *10), F ( s ) = exp( - *l O s ) / s *

M e t h o d D u r b i n L A P I N L A P I N

C o d e ( v . T = 5) 2 - 1 - 0 - 1 1 M A N = 0 1 - 2 - 1 - 2

NS1 + N S 2 2000 50 60 350 + 150

t R e a l a b s o l u t e e r r o r t R e a l a b s o l u t e e r r o r

5 0.6 D - 0 2 0.4 D - 0 5 0.5 D - 0 6 9.0 0.1 D - 1 3

6 0.6 D - 0 2 0.1 D - 0 4 0.5 D - 0 6 9.2 0.3 D - 1 2

7 0.6 D - 0 2 0.6 D - 0 4 0.5 D - 0 6 9.4 0.8 D - 1 4

8 0.7 D - 0 2 0.5 D - 0 3 0.5 D - 0 6 9.6 0.1 D - 0 9

9 0.7 D - 02 0.1 D - 01 0.6 D - 06 9.8 0.6 D - 05

10 0.6 D - 02 0.5 D - 04 a 0.6 D - 06 a 10.0 0.3 D - 09 a

11 0.5 D - 0 2 0.3 D - 0 4 0.8 D - 0 5 10.2 0.1 D - 0 5

12 0.6 D - 0 2 0.6 D - 0 2 0.5 D - 0 6 10.4 0.8 D - 1 0

13 0.6 D - 0 2 0.1 D - 0 2 0.5 D - 0 6 10.6 0.2 D - 0 9

14 0.6 D - 0 2 0.1 D - 0 2 0.5 D - 0 6 10.8 0.8 D - 1 0

15 0.6 D - 0 2 0.2 D - 0 2 0.5 D - 0 6 11.0 0.2 D - 1 2

C P U - t i m e

sec ( f o r C P U - t i m e

t = 1(1)20) 2.02 0.06 0.32 sec 13.70

a C F M w a s u s e d to a c c e l e r a t e c o n v e r g e n c e .

124 *G. Honig, U. Hirdes / Laplace transform *

Table 5

*f ( t *= ( - t 3 + 9t 2 - 18t + 6)/6, *F(s) = (s - *1)3/s 4

Method LAPIN LAPIN

Code 1 - 2 - 0 - 2 (without ' Korrektur') 1 - 2 - 1 - 2 (with "Korrektur')

NS1 + NS2 100 60 + 40

t Absolute error Vop v - T Absolute error roy r • T

Real By (33) Real By (33)

**1 ** **0.1 D - 1 0 ** **0.2 D - 1 0 ** 12.2 0.3 D - 1 2 0.4 D - 1 5 9.4

3 0.9 D - 1 0 0.1 D - 1 0 14.7 0.6 D - 1 2 0.5 D - 1 3 9.9

6 0.1 D - 0 9 0.2 D - 1 0 15.4 0.2 D - 1 2 0.2 D - 1 2 10.1

0 0.5 D - 1 0 0.7 D - 1 0 16.0 0.8 D - 1 3 0.1 D - 1 2 10.6

CPU-time

s e c

( t = 1(1)10) 1.49 0.81

Table 6

*f ( t ) *= 1 - e x p ( t ) erfc(~/t); *F(s) = 1/(sVrs + *1))

Method LAPIN

Code 1 - 2 - 0 - 2

NS1 Absolute error at t = 1 Vop x . T

Real By (33)

10 0.103 D - 0 4 0.416 D - 0 5 6.61

20 0.352 D - 10 0.900 D - 1 0 11.96

30 0.356 D - 1 1 0.801 D - 1 1 13.16

c a l c u l a t e Vop T f o r a g i v e n N S 1 ( a n d N S 2 ) m a k e s t h e i n v e r s i o n m e t h o d s b a s e d o n F o u r i e r s e r i e s e x p a n s i o n s a c c u r a t e a n d e f f i c i e n t .

**Appendix A. The F O R T R A N subroutine L A P I N **

*A.1. Purpose *

S u p p o s e *f ( t ) * is a r e a l f u n c t i o n a n d *F ( s ) *i t s L a p l a c e t r a n s f o r m , t h a t h a s n o s i n g u l a r i t i e s t o t h e r i g h t o f
t h e o r i g i n , t h e n t h e s u b r o u t i n e L A P I N c a l c u l a t e s f r o m *F ( s ) *f o r a p r o g r a m m e r s - c h o s e n set o f t - v a l u e s , t h e
c o r r e s p o n d i n g v a l u e s f o r *f ( t ) , * u s i n g t h e a b o v e d e s c r i b e d a l g o r i t h m s .

*A.2. Usage *

C A L L L A P I N ( F , T 1 , T N , N , I M A N , I L A P I N , I K O N V , N S 1 , N S 2 , I C O N , I K O R , C O N , **H, E, **I E R ,
1 O U T )

*A.3. Description o f the parameters *

T y p e : I N T E G E R * 4 N , I M A N , I L A P I N , I K O N V , N S 1 , N S 2 , I C O N , I E R , 1 O U T D O U B L E P R E C I S I O N T 1 , T N , C O N , H , E

I n p u t p a r a m e t e r s : T 1 , T N , N , I M A N , I L A P I N , I K O N V , N S 1 , N S 2 , I C O N , 1 K O R , C O N O u t p u t p a r a m e t e r s : H , I E R

*G. Honig, U. Hirdes / Laplace transform * 1 2 5

*A. 3.1. Significance of the parameters *

T1, T N N

I M A N

I K O N V NS1

N S 2 I C O N

I K O R C O N H

Laplace t r a n s f o r m - - e x t e r n a l subroutine, to declare as E X T E R N A L in the calling p r o g r a m a n d written b y the user,

S U B R O U T I N E F (SR, SI, F R , FI), D O U B L E P R E C I S I O N SR, SI, F R , FI, SR real part of s,

SI i m a g i n a r y part of s,

F R real part of Laplace transform, F I i m a g i n a r y part of Laplace-transform,

lower a n d u p p e r b o u n d s of the interval in which *f ( t ) *shall be a p p r o x i m a t e d ,
n u m b e r of t-values for which f ( t ) shall be c o m p u t e d , the t-values are given b y

t k = T I + ( T N - T 1 ) k / ( N + I ) , k = l ( 1 ) N ,

= 0 all further parameters are placed automatically, except NS1, which has to be set equal to 60;

= 1 a m a n u a l choice of the following parameters is possible,

I L A P I N = 1 implies the application of the a p p r o x i m a t i o n f o r m u l a with T = t, I L A P I N = 2 with T = T N ( T = 2 T N for the ' K o r r e k t u r ' - t e r m s ) ,

if I K O N V - 1 the m i n i m u m - m a x i m u m method, if I K O N V - - 2 the e-algorithm is used to accelerate the convergence of the series,

n u m b e r of function evaluations of *F(s) *used to a p p r o x i m a t e
*f ( t ) ( f N ( t ) *= f N s x ( t ) ) (NSa = 60 if I M A N = 0),

n u m b e r of function evaluations of F(s) used to a p p r o x i m a t e the ' K o r r e k t u r ' - t e r m *f ( 2 T + t) *
*( f u ( 2 T + t ) = f N s z ( 2 T + *t)),

= 0 n o o p t i m a l choice of the free parameters,

= 1 optimal choice of the free parameters for t = *tiN~2], *

= 2 optimal choice of the free parameters for t = t k, k = 1(1)N,

= 0 no application of the ' K o r r e k t u r ' - m e t h o d ,

= 1 application of the ' K o r r e k t u r ' - m e t h o d ,

b y the choice o f C O N = *vT, *the free p a r a m e t e r v is determined in case of I C O N = 0,
matrix of d i m e n s i o n 6 by N, contains on return in the

*1st row *the a p p r o x i m a t i o n values for f ( t k ) , k = l(1)n,
*2nd row *work area,

*3rd row *the c o m p u t e d optimal p a r a m e t e r C O N o v r = r o y r • T for t = t k, k = I ( 1 ) N , ( C O N o v r = 18
if the truncation error (23) is zero, C O N o v r = 1 if the discretization error (24) is zero or of very
small absolute value),

*4th row *the code for the used acceleration m e t h o d :

= 0 n o acceleration of convergence,

= 1 M I N I M A X ,

= 2 E P A L ,

= 4 C F M ,

= P >~ 4 E P A L , overflow occurred after P - 2 iterations (NS1 = P - 2 values of *F ( s ) *are used for
the a p p r o x i m a t i o n of f ( t ) ) ; note: the o u t p u t p a r a m e t e r on the 4th row of H m a y be different
f r o m the input p a r a m e t e r I K O N V . T h e subroutine L A P I N always uses

(a) C F M *i f f u ( t ) *is m o n o t o n o u s i n N ,

(b) M I N I M A X if the application of E P A L falsifies the results,

(c) n o acceleration m e t h o d if only 1 or 2 stationary values of *f N ( t ) *as a function of N are
found,

*5th row *absolute error calculated by f o r m u l a (33) if I C O N = 2 a n d I L A P I N = 1; if I C O N = 1 the
error estimation is valid only at t ^{= } *tiN/21, *

*6th row *contains a control n u m b e r of up to 6 digits *( n l n 2 n 3 n a n s n 6 ) , * the digits refer to the used

126 *G. Honig, U. Hirdes / Laplace transform *

E I E R

I O U T

auxiliary quantities,

the definition of t o and t~):

n I ~ f r ~ s l ( t 0 ) , ( C O N = 20), n 2 --)fN2s,(to), ( C O N = 18), n 3 --)fNs,(2t, + to),

n 4 ~ f N S a ( 4 t l + to),

the ' K o r r e k t u r ' - t e r m fNS2 (2tl + t0) a n d to f N s l ( t ) as follows (see Table 2 for

n5 ----~fNS1 ( 8 t l q" 10),

n 6 = n61 -4- n62 , /-/61 --*fNsl(t0), n62 ~ f N s z ( 2 t l + t o ) .

If one of the digits is zero the application of the e-algorithm falsifies the results of the c o r r e s p o n d i n g auxiliary quantity ( M I N I M A X is used). Otherwise the digits are equal to 1 (besides n62 , which is equal to 2).

N o t e : n 3 = 0 always if I K O R = 1, H a : n 5 = 0 always if I K O R = 0.

W h e n e v e r a zero in the control n u m b e r indicates the falsification of the results, it is r e c o m m e n d e d to increase NS1 or NS2.

matrix of dimension 3 by NS1, work area, error parameter, control of input data:

= 0 no error,

**= 1 ** **T N < **T1,

= 1 0 N < I ,

= 100 I M A N < 0 or I M A N > 1,

= 1000 I L A P I N < 1 or I L A P I N > 2,

= 10 000 I K O N V < 1 or I K O N V > 2,

= 100000 NS1 < 1 or (NS2 < 1 and I K O R = 1) or (NS1 4= 60 and I M A N = 0),

= 1000 000 I C O N < 0 or I C O N > 2 or ( I C O N = 2 a n d I L A P I N = 2),

= 1 0 0 0 0 0 0 0 I K O R < I o r I K O R > l ,

= 1 0 0 0 0 0 0 0 0 C O N ~< 0.

E.g. I E R = 11 means that 1 a n d 10 occured, o u t p u t unit number.

*A. 4. Subroutines *

L A P I N calls the subroutine L A P I N 2 to calculate the Laplace-inversion with optimal parameters and in case of w r o n g input d a t a the subroutine E R R O R .

T h e subroutine F m u s t be a d d e d by the user.

**A c k n o w l e d g m e n t **

T h e authors wish to thank K. Mika for valuable suggestions a n d c o m m e n t s a n d for careful revision of the manuscript. W e thank U. F u n k and G. H o f e m a n n for helpful assistance with the numerical calcula- tions.

F o r m a k i n g possible the stay at the D e p a r t a m e n t o de I n f o r m h t i c a of Pontificia Universidade Cat61ica (PUC), Rio de Janeiro, the first a u t h o r wishes to thank the director of the department, Prof. Carlos Jos6 Pereira de Lucena, the Conselho N a c i o n a l de D e s e n v o l v i m e n t o Cientificio e Tecnolbgico (CNPq), Gesell- schaft for M a t h e m a t i k u n d D a t e n v e r a r b e i t u n g Bonn ( G M D ) , P U C Rio de Janeiro and K F A JiJlich.

*G. Honig, U. Hirdes / Laplace transform * 127

**References **

*[1] P. Albrecht and G. Honig, Die numerische Inversion der Laplace-Transformierten, Angew. lnformatik 8 (1977) 336-345. *

*[2] K.S. Crump, Numerical inversion of Laplace transforms using a Fourier series approximation, J. A C M 23 (1) (1976) 89-96. *

*[3] G. Doetsch, Handbuch der Laplace-Transformation, Bd. L II, I I I (Birkh~iuser, Basel, 1950). *

[4] H. Dubner and J. Abate, Numerical inversion of Laplace transforms by relating them to the finite Fourier cosine transform, J.

*A C M 15 (1) (1968) 115-123. *

*[5] F. Durbin, Numerical inversion of Laplace transforms: an effective improvement of Dubner and Abate's method, Comput. J. 17 *
(4) (1973) 371-376.

[6] G. Honig, Zur numerischen LOsung partieller Differentialgleichungen mit Laplace-Transformation, Diss. Univ. Dortmund D290, Berichte der Kernforschungsanlage Jtilich-Jiil-1550, 1978.

*[7] K. Kehr, G. Honig and D. Richter, Stochastic theory of spin depolarization of muons diffusing in the presence of traps, Z. Phys. *

B 32, (1978) 49-58.

[8] R.M. Simon, M.T. Stroot and G.H. Weiss, Numerical inversion of Laplace transforms with application to percentage labeled
*experiments, Comput. Biomed. Re{;. 6 (1972) 596-607. *

[9] F. Veillon, Quelques m6thodes nouvelles pour le calcul num+rique de la transform~ inverse de Laplace, Th. Univ. de Grenoble, 1972.

[10] G. Warzee, Application de la transform~e de Laplace et de la m6thode des 61+ments finis ~ la r6solution de l'6quation de
*conduction thermique instationnaire, C.R. Acad. Sei. Paris Sbr. A278 (1974) 1265-1266. *

**Subroutine LAPIN **

S U B R O U T I N E L A P I N C F , T I , T N , N , IMAN, ILAPIN, I K O N V , N S 1 , N S 2 , 1 C O N ,

x I K O R , C O N , H , E , IER, IOUT)

D O U B L E P R E C I S I O N A B R N , A B S F , C O N , C O N O P T , C O N I , C O N 2 , D E L , E , x F N , F N S I , F R E A L , F I M A G j H , P I , R A C C , R N S U M , RNSUMK, x T , T A , T B , T K ' , T N , T 0 , T I , V I , V 2 , W

I N T E G E R H M O N O E X T E R N A L F

DIMENSION HC6,N),EC3,NSI)

COMMON /CLAPIN/ TA,TB,T0,CONOPT,ABSF,LVAL,HMONO I N I T I A L I Z E ARRAY H

DO 10 I = I , N DO 10 J = l , 6 H C J , I ) = 0.O0

10 CONTINUE

IER **= 0 **

IF C T N . L T . T 1 ) IER = 1 IF C N . L T . I ) IER = IER+I0 IF CIMAN.EQ.0) GO TO 100 IF CIMAN.EQ.I) GO TO 110 IER = IER+100

GO TO 120

PARAMETERS FOR IMAN = 0

100 IF CNSI.NE.60) IER : IER+I00000 IF C I E R . N E . 0 ) GO TO 120

ILAPIN = 1
IKONV = 2
ICON = 1
IKOR = 0
NS2 = 0
**GO ** TO 200
IMAN = 1

110 IF C I L A P I N . L T . 1 .OR. I L A P I N . G T . 2 ) IER = IER+I000 IF CIKONV.LT.1 .OR. IKONV.GT.2) IER = IER+10000 IF C N S I . L T . I . O R . CNS2.LT.1 . A N D . IKOR.EQ.1)) IER = IF CICON.LT.0 .OR. ICON.GT.2 . O R . CICON.EQ.2 .AND.

:: IER = IER+1000000

IF CIKOR.LT.0 .OR. IKOR.GT.1) IER = IER+I00O0000 IF CICON.EQ.0.AND.CON.LE.0.D0) IER = IER+100OO0O00

IER+I00000 I L A P I N . E Q . 2 ) )

128 *G. Honig, U. Hirdes / Laplace transform *

IF ( I E R . E Q . 0 ) GO TO 200

120 CALL E R R O R ( I E R , T 1 , T N , N, IMAN, ILAPIN, I K O N V , N S I , N S 2 , 1 C O N , IKOR,

* CON, I O U T )

RETURN

200 PI = 4.D0*DATAN(I.D0) CON1 = 2 0 . D 0

CON2 = C O N 1 - 2 . D 0 ABSF = 0 . D 0

J3 = ( 3 - 1 C O N ) / 2 + N * ( I C O N / 2 ) TA = T I

TB = TN

DO 830 L3=1,J3 LVAL = L3 HMONO = 0 KOR1 = IKOR

COMPUTATION OF THE OPTIMAL PARAMETERS 210 IF ( I C O N - l ) 2 1 5 , 2 3 0 , 2 2 0

215 JUMP = 0

CALL LAPIN2(F,T1,TN,N, ILAPIN, IKONV,NSI,NS2,1CON, IKOR,CON,H,E,JUMP) GO TO 830

220 TA = T I + F L O A T ( L 3 ) * ( T N - T I ) / F L O A T ( N + I ) TB = TA

230 NH = N / 2

T = T A + ( T B - T A ) * F L O A T ( N H ) / F L O A T ( N + I ) TK = F L O A T ( 2 - 1 L A P I N ) * T + F L O A T ( I L A P I N - I ) * T B COMPUTATION OF THE TRUNCATION ERROR (RNSUM)

TO = T CON = CONI JUMP = 1

CALL L A P I N 2 ( F , T I , T N , N , I L A P I N , I K O N V , N S 1 , N S 2 p l C O N p l K O R , C O N , H , E , J U M P ) 240 FN : H ( I , L 3 )

FNS1 = E ( 1 , N S I ) CON = CON2 JUMP : 2

CALL L A P I N 2 ( F , T I , T N , N, ILAPIN, IKONV,NSI,NS2,1CON, IKOR,CON,H,E,JUMP) 250 IF ( F N . N E . H ( 1 , L 3 ) . A N D . F N S I . N E . E ( 1 , N S 1 ) ) GO TO 255

CONOPT = CON ABSF = 0.D0 GO TO 320

255 RNSUM = T K * ( F N - H ( I , L 3 ) ) / ( D E X P ( C O N I ) - D E X P ( C O N 2 ) ) I F ( I L A P I N . E Q . 2 ) GO TO 260

COMPUTATION OF THE ACCELERATION FACTOR ( D E L ) RACC = T * ( F N S 1 - E ( I , N S I ) ) / ( D E X P ( C O N I ) - D E X P ( C O N 2 ) ) DEL = RNSUM/RACC

260 IF ( I K O R . E Q . 1 ) GO TO 280 OPTIMAL PARAMETERS (METHOD A)

TO = 2.D0*TK+T CON = CONI/4.D0 JUMP = 3

CALL L A P I N 2 ( F , T I , T N , N , ILAPIN, IKONV,NSI,NS2,1CON, IKOR, CON,H,E,JUMP) 270 FN = H ( I , L 3 )

CONOPT = -TK/C2.D0*TK+T)*DLOG(DABS(RNSUM/(TK*FN))) GO TO 310

O P T I M A L PARAMETERS FOR THE KORREKTUR METHOD (METHOD A ) 280 TO = 4 . D 0 * T K + T

CON = C O N 1 / 4 . D 0 JUMP = 4

CALL L A P I N 2 ( F p T l p T N , N p l L A P I N , I K O N V s N S 1 , N S 2 p I C O N p l K O R p C O N p H , E, J U M P ) 290 FN = H C 1 , L 3 )

TO = 8 . D 0 * T K + T JUMP = 5

*G. Honig, U. Hirdes / Laplace transform * 129

C A L L L A P I N 2 ( F , T I , T N , N , ILAPIN, I K O N V , N S I , N S 2 , I C O N , I K O R , C O N , H , E , J U M P ) 300 FN = F N - H ( I , L 3 )

C O N O P T = - T K / ( 4 . D 0 * T K + T ) * D L O G ( D A B S ( R N S U M / ( T K * F N ) ) ) C

C O P T I M A L P A R A M E T E R S ( M E T H O D B) 310 IF ( I L A P I N . E Q . 1 ) GO TO 315

A B S F = D A B S ( D E X P ( C O N O P T ) * R N S U M * 2 . D O / T K ) GO TO 320

315 VI = C O N 1 / T V2 = C O N O P T / T W = P I * F L O A T ( N S I ) / T C A L L F ( V 2 , W , F R E A L , F I M A G ) R N S U M K = R N S U M * F R E A L C A L L F ( V I , W , F R E A L , F I M A G ) R N S U M K = R N S U M K / F R E A L

A B R N = ( R N S U M K - R N S U M ) / ( V 2 - V 1 )

C O N O P T = - D L O G ( D A B S ( ( A B R N / T + R N S U M K ) / ( T * F L O A T ( I K O R * 2 + 2 ) * F N ) ) ) /

* F L O A T ( 3 + 2 * I K O R ) V l = V2

V 2 = C O N O P T / T

C A L L F ( V 2 , W , F R E A L , F I M A G ) RNSUMK = R N S U M K * F R E A L C A L L F ( V I , W , F R E A L , F I M A G ) RNSUMK = R N S U M K / F R E A L

A B S F : D E X P ( C O N O P T ) / T : : D A B S ( R N S U M K ) +

* D A B S ( D E X P ( - 2 . D 0 * C O N O P T ) * F N * F L O A T ( I K O R - I )

* + D E X P C - 4 . D 0 * C O N O P T ) * F N * P L O A T ( I K O R ) ) C

3 2 0 I F ( C O N O P T . L E . 0 . D 0 ) CONOPT = I . D 0 J U M P = 6

C A L L L A P I N 2 ( F , T I , T N , N , ILAPIN, IKONV, N S I , N S 2 , 1 C O N , IKOR, C O N , H , E , J U M P ) 830 C O N T I N U E

R E T U R N E N D

S U B R O U T I N E L A P I N 2 ( F , T I , T N , N , ILAPIN, I K O N V , N S I , N S 2 , ICON, I K O R , C O N ,

* H, E, J U M P )

C

C L A P L A C E - I N V E R S I O N W I T H O P T I M A L P A R A M E T E R S C

D O U B L E P R E C I S I O N A, A B S F , B , C O N , C O N O P T , D E L T A , D I V I , E , E I , E 2 , E 3 , E I N S ,

* F A K T O R , F I M A G , F R E A L , H, P I , P I T , P I T E , R A L , S U I M , SURE,

* TA, TB, TE, TL, TM, TM1, TN, TT, T0, TI, XI, X2, X3,

* V , W , Y I , Y 2 , Y 3

I N T E G E R R I C H T , R ICHTA, H M O N O E X T E R N A L F

D I M E N S I O N H(6, N), E(3, N S 1 ) , E3 (3)

C O M M O N / C L A P I N / T A , T B , T 0 , C O N O P T , A B S F , L 3 , H M O N O C

PI = 4 . D O * D A T A N ( I . D 0 ) C

IF ( J U M P . E Q . 0 ) GO TO 360 IF ( J U M P . L T . 6 ) GO TO 370 C

TO = TA TT = TB 11 : L3

J1 = ( 2 - 1 C O N ) * N + ( I C O N - I ) * L 3 C O N = C O N O P T

K O R 1 = IKOR GO TO 380 C

360 TO = TI TT -- TN I i = i J l = N K O R I = I K O R J U M P : 6 GO TO 380

130 *G. Honig~ U. Hirdes / Laplace transform *

370 KOR1 = 0 TT = TO

I i = L3 J1 = L3 C

380 DELTA = ( T T - T 0 ) / F L O A T ( N + I ) NSUM = N S l

C

C COMPUTATION OF THE T - V A L U E S FROM T 0 , T T D O 8 0 5 K 2 = 1 , 2

IF ( I L A P I N . E Q . I ) G O T O 4 2 0 T E : P L O A T ( K 2 ) * T T

V = C O N / T E

C A L L F ( V , 0 . D 0 , P R E A L , F I M A G ) RAL = -0.SD0*PREAL PITE = PI/TE C

405 DO 410 L=I,NSUM W = FLOAT(L-1)*PITE CALL F(V,W,FREAL,FIMAG) E ( 2 , L ) = FREAL

E ( 3 , L ) = FIMAG 410 CONTINUE C

420 DO 800 K I = I 1 , J 1 C

TL = T0+FLOAT(K1)*DELTA IF ( I L A P I N . E Q . 2 ) GO TO 440 C

IF ( K 2 . E Q . 2 ) TL = 3.D0*TL V = CON/TL

FAKTOR = DEXP(V*TL)/TL

CALL F(V,0.D0,FREAL,FIMAG) R A L = - 0 . 5 D 0 * F R E A L

P I T = P I / T L E I N S = I . D 0 S U R E = 0 . D 0 C

C M E T H O D O F D U R B I N 4 2 5 D O 4 3 0 L = I , N S U M

W = F L O A T ( L - 1 ) * P I T C A L L F ( V , W , F R E A L , F I M A G ) S U R E = S U R E + F R E A L * E I N S E I N S = - E I N S

E ( I , L ) = F A K T O R * ( R A L + S U R E ) 4 3 0 C O N T I N U E

G O T O 4 6 0 C

4 4 0 IF ( K 2 . E Q . 2 ) T L = T L + T E F A K T O R = D E X P ( V * T L ) / T E S U R E = 0 . D O

S U I M = 0 . D 0 D O 4 5 0 L = I , N S U M W = F L O A T ( L - I ) * P I T E

SURE = S U R E + E ( 2 , L ) * D C O S ( W * T L ) SUIM = S U I M + E ( 3 , L ) * D S I N ( W * T L ) E ( I , L ) = F A K T O R * ( R A L + S U R E - S U I M ) 450 CONTINUE

C

C SEARCH FOR STATIONARY VALUES 460 NMAX = N S U M * 2 / 3

MONOTO = 0 K = 0

R I C H T A = D S I G N ( 1 . S D 0 , ( E ( I , N S U M ) - E ( I , N S U M - I ) ) ) DO 500 L = I , N M A X

J = NSUM-L

R I C H T = D S I G N ( 1 . 5 D 0 , ( E ( I , J ) - E ( 1 , J - 1 ) ) ) IF ( R I C H T . E Q . R I C H T A ) G O T O 5 0 0 K = K + I

E 3 ( K ) = E ( I , J ) R I C H T A = R I C H T

I F ( K . E Q . 3 ) GO TO 510

*G. Honig, U. Hirdes / Laplace transform * 131

500 CONTINUE

IF ( K . E Q . 0 ) G O T O 7 0 0 HCK2,K1) = ECI,NSUM) IF ( K 2 . E Q . I ) H C 4 , K 1 ) = 0 G O T O 790

5 1 0 K E = 2

IF ( ( E 3 ( K E ) - E ( 1 , J ) ) : : F L O A T ( R I C H T A ) . G T . 0 . D 0 ) G O T O 560 J M I N = N S U M / 3

J M A X = J - 1

D O 540 J J = J M I N , J M A X J = J - 1

R I C H T = D S I G N ( I . 5 D 0 , ( E ( I , J ) - E ( 1 , J - 1 ) ) ) IF (RICHT.EQ.RICHTA) GO TO 540 RICHTA = RICHT

KE - 3-KE

IF ((E3(KE)-E(I,J));CFLOATCRICHTA).GT.0.D0) GO TO 560 54-0 CONTINUE

550 MONOTO = 1

IF (IKONV.EQ.2) GO TO 630 C

C MINIMUM-MAXIMUM METHOD (MINIMAX)

560 H ( K 2 , K I ) = ( E 3 ( 1 ) + E 3 ( 3 ) ) / 4 . D 0 + E 3 ( 2 ) / 2 . D 0 IF ( K 2 . E Q . I ) H ( 4 , K I ) = 1

G O T O 7 9 0 C

C E P S I L O N A L G O R I T H M ( E P A L ) 6 3 0 K : O

N S U M M I = N S U M - 1 E 2 = E ( I , I ) D O 6 6 0 L = I , N S U M M I El = E ( I , I ) T M = 0 . D 0 L P 1 -- L + I D O 6 5 0 M = I , L M M = L P 1 - M T M I -- E ( I , M M )

D I V I -- E ( I ~ M M + I ) - E ( I , M M )

IF ( D A B S ( D I V I ) . G T . I . D - 2 0 ) G O T O 6 4 0 K = L

G O T O 6 7 0

6 4 0 E ( I , M M ) = T M + I . D O / D I V I T M = T M i

6 5 0 C O N T I N U E E 2 = E 1 6 6 0 C O N T I N U E C

6 7 0 IF CDABSCEI).GT.DABSCE2)) El = E2

IF ( D A B S ( E C I , 1 ) ) . G T . D A B S ( E 1 ) ) E ( 1 , 1 ) = E1 H ( K 2 , K 1 ) : E ( I , I )

IF ( K 2 . E Q . I ) H ( 4 , K 1 ) = K + 2 G O T O 7 9 0

C

C C U R V E F I T T I N G ( C F M ) 700 Xl : FLOAT(NSUM)-2.DO

X2 = FLOATCNSUM)-I.D0 X3 = FLOATCNSUM)-0.D0 Y1 = ECI,NSUM-2) Y2 = EC1,NSUM-1) Y3 = EC1,NSUM)

B = ( ( Y 3 - Y I ) * X 3 * X 3 * C X l + X 2 ) / C X l - X 3 ) x _ ( y 2 _ y i ) * X 2 , X 2 x C x I + X 3 ) / C x I _ X 2 ) ) / ( X 3 _ X 2 )

A = ( ( Y 2 - Y I ) - B * ( X I - X 2 ) / ( X l * X 2 ) ) * ( X 2 * X 2 * X I * X 1 ) / ( X I * X I - X 2 * X 2 ) H ( K 2 , K I ) = Y I - ( A / X I + B ) / X I

M O N O T O = 1

IF ( K 2 . E Q . 1 ) HC4, K1) = 3

790 H C 3 , K I ) = CON HC5,K1) = ABSF

HMONO : HMONO+K2*MONOTO*I0**C6-JUMP) IF (JUMP.LT.6) GO TO 800

HC6,K1) = FLOATCC2-K2)*HMONO) +