# Numerical Real Inversion of the Laplace Transform by Using Multiple-Precision Arithmetic.

1 Introduction

The aim of the present paper is numerical realization of the real inversion of the Laplace transform

F(s) = L f (s)= [[integral].sup.[infinity].sub.0][e.sup.-st]f (t) dt,

which appears in various fields such as mathematical science, mechanics, engineering, and mathematical finance. Particularly we extend the framework proposed by Saitoh  to weighted Hilbert spaces for the sake of calculating practical problems. From the viewpoint of numerical computations, the growth of computational errors is crucial for reliable computations since the real inversion of the Laplace transform is unstable. To overcome the difficulty multiple-precision arithmetic also plays essential role in our algorithm.

The real inversion of the Laplace transform is the problem to find the original function f = [L.sup.-1] F from the image function F(s) which is given on the non-negative real axis {s [greater than or equal to] 0}. In general, [L.sup.-1] F is solved by using Laplace transform tables or the Bromwich integral [1, 13]. However the former is applicable only to cases where the image function F is analytically calculated, and the latter requires a priori information of the distribution of poles of F in the complex plane.

On the other hand, Saitoh  proposed an inversion formula which required only the information of F on the non-negative real axis. However it is applicable only to a few problems since the derivative f' of the original function is assumed to decay exponentially and the image function F belongs to [L.sup.2](0, [infinity]) in . To overcome the limitations, we introduce weighted Hilbert spaces.

The proposed method also adopts the Tikhonov regularization as similar as Saitoh's strategy . To detect the singularities of the exact solution by the regularization method, we need a small regularization parameter to realize accurate approximation. It also approximates instability of the problem at the same time. That means computational errors grow rapidly in its numerical processes and computation fails. To overcome the difficulty we apply multiple-precision arithmetic, and thus achieve reliable numerical inversions.

2 Instability of Real Inversion of the Laplace Transform

For a positive number [omega], the Laplace transform of the function [f.sub.[omega]](t) = [square root of [omega]]sin [omega]t is

[F.sub.[omega]](s) = L[[F.sub.[omega]]](s) = [[omega][square root of [omega]]]/[[s.sup.2][[omega].sup.2]], and it satisfies

[mathematical expression not reproducible]

Therefore [||[F.sub.[omega]]||.sub.[infinity]] [right arrow] 0 as [omega] [right arrow] [infinity] while [||[F.sub.[omega]]||.sub.[infinity]] [right arrow] [infinity]. This means that inversion of the Laplace transform [L.sup.-1] is not stable with the supremum norm. This leads to numerical instability of its direct discretization, and numerical computation fails with the standard computer environments due to rapid increase of various kinds of computational errors.

3 Tikhonov Regularization on Reproducing Kernel Hilbert Spaces

Let X and Y be Hilbert spaces, and L : X [right arrow] Y be an injective and bounded linear operator. For a linear equation Lf = g and a fixed positive number [alpha], there exists a unique function [f.sub.[alpha]] [member of] X which is the minimizer of the Tikhonov functional [J.sub.[alpha]](f) = [alpha] [||f||.sup.2.sub.X] + [||Lf - g||.sup.2.sub.Y], and it is also a unique solution to ([alpha]I + L*L)[f.sub.[alpha]] = L*g .

In addition, if X admits a reproducing kernel, we have the next theorem.

Theorem 3.1. Let X be a reproducing kernel Hilbert space of functions on a set A, and let [K.sub.t] = K(*, t) be its reproducing kernel. Then, for any positive number [alpha] and any point t [member of] A, there exists a unique solution [H.sub.[alpha],t] = [H.sub.[alpha]](*, t) [member of] Y to the equation

([alpha]I + LL*)[H.sub.[alpha],t] = [LK.sub.t], (3.1)

and the Tikhonov regularized solution [f.sub.[alpha]] [member of] X to Lf = g is given by

[f.sub.[alpha]](t) = [(g, [H.sub.[alpha],t]).sub.Y], t [member of] A.

Proof. By virtue of the reproducing property, we have

[f.sub.[alpha]](t) = [([f.sub.[alpha]], [K.sub.t]).sub.X] = [([([alpha]I + L*L).sup.-1]L*g, [K.sub.t]).sub.X] = [(g,L[(aI + L*L).sup.-1][K.sub.t]).sub.Y].

Let [G.sub.[alpha],t] = [G.sub.[alpha]](*, t) = [([alpha]I + L*L).sup.-1][K.sub.t] and [H.sub.[alpha],t] = [LG.sub.[alpha],t]. Then [f.sub.[alpha]](t) = [(g, [LG.sub.[alpha],t]).sub.Y] = [(g, [H.sub.[alpha],t]).sub.Y]. Applying L to the both sides of ([alpha]I + L*L)[G.sub.[alpha],t] = [K.sub.t] leads to ([alpha]I + LL*)[LG.sub.[alpha],t] = [LK.sub.t], which is (3.1).

Remark 3.2. The equation (3.1) shows that [H.sub.[alpha],t] is the Tikhonov regularized solution to the equation L*H = [K.sub.t].

4 Real Inversion of the Laplace Transform on Reproducing Kernel Hilbert Spaces

Let w(t) be a positive continuous function defined on t > 0. We will denote by [H.sub.w] the set of all absolutely continuous functions f defined on [0, [infinity]) satisfying f (0) = 0 and

The space [H.sub.w] is a reproducing kernel Hilbert space whose reproducing kernel is

[K.sub.t](s) - [[integral].sup.min(s,t).sub.0]w[([xi]).sup.-1]d[xi].

In fact, since

[mathematical expression not reproducible]

the reproducing property of [K.sub.t] is shown as

[mathematical expression not reproducible]

For a positive function u [member of] C[0, [infinity]), we consider the weighted [L.sup.2] space

[mathematical expression not reproducible]

We consider the operator L : [H.sub.w] [right arrow] [L.sup.2.sub.u] defined by

Lf (s) = s L[f](s) = L[f'](s) = [[integral].sup.[infinity].sub.0](t)[e.sup.-st]dt, f [member of] [H.sub.w], s > 0.

Proposition 4.1 (). If positive functions w, u [member of] C(0, [infinity]) satisfy

I = [[integral].sup.[infinity].sub.0][[integral].sup.[infinity].sub.0][e.sup.-2st]w[(t).sup.-1]u(s)dt ds < [infinity]

then L : [H.sub.w] [right arrow] [L.sup.2.sub.u] is bounded.

Proof. It follows from the definition of the operator L that

[mathematical expression not reproducible]

Therefore

[mathematical expression not reproducible]

It is also known that L : [H.sub.w] [right arrow] [L.sup.2.sub.u] is compact under (4.1), which also explains ill-posedness of the problem.

Theorem 4.2. If positive functions w, u [member of] C(0, [infinity]) satisfy (4.1), then for any positive numbers [alpha] and t, there exists a unique solution [H.sub.[alpha]](*,t) [member of] [L.sup.2.sub.u] to the integral equation of the second kind

[alpha][H.sub.[alpha]](s, t) + [[integral].sup.[infinity].sub.0][H.sub.[alpha]]([sigma], t)L[1/w](s + [sigma])u([sigma])d[sigma] = [LK.sub.t](s), s > 0, (4.2)

and the Tikhonov regularized solution to Lf = g, g [member of] [L.sup.2.sub.u] is given by

[mathematical expression not reproducible]

Proof. We need only to show that

LL*y{s) = [[integral].sup.[infinity].sub.0] y([sigma]) [??][1/w] (s + [sigma])u([sigma]) d[sigma], y [member of] [L.sup.2.sub.u], s > 0.

For y [member of] [L.sup.2.sub.u] and s > 0 fixed,

L(L*y)(s)= [??][(L*y)'](s)

= [[integral].sup.[infinity].sub.0](L*y)'(t)[[e.sup.-st]]/w(t)dt (4.3)

(=) [[integral].sup.[infinity].sub.0] (L*y)'(t)d/dt [L.sub.[sigma][[K.sub.t]([sigma])](s)w(t)dt (4.4)

=[(L*y(t), L[[K.sub.t]](s)).sub.[H.sub.w] (as functions of t)] (4.5)

= (y, [L.sub.t][L[[K.sub.t]](s)]).sub.[L.sup.2.sub.u], (4.6)

where [L.sub.[sigma]] means that L operates with respect to [sigma]. From (4.3) to (4.4) we are using

[L.sub.[sigma]][[K.sub.t]([sigma])]{s) = [[??].sub.[sigma]] [d/d[sigma][K.sub.t]([sigma],[>>])] (s) = [[??].sub.[sigma]] [1.sub.[0,T]([sigma])/w([sigma])] (s) = [[integral].sup.t.sub.0] [e.sup.-s[sigma]]/w([sigma])d[sigma].

We also notice that L[[K.sub.t]](s) belongs to [H.sub.w] as a function of t because L[[K.sub.0]](s) = 0 and

[[integral].sup.[infinity].sub.0][|d/dt L[[K.sub.t]](s)|.sup.2] w(t)dt= [[integral].sup.[infinity].sub.0] [e.sup.-2st][w(t).sup.-1] dt < [infinity]

by the condition (4.1). Therefore (4.4) can be written as (4.5). Next we compute [L.sub.t][L[[K.sub.t]](s)]. For [sigma] > 0 fixed,

[mathematical expression not reproducible]

Consequently,

L(L*y)(s) = (y,[L.sub.t][L[[K.sub.t]](s)])L.sup.2.sub.u] = [[integral].sup.[infinity].sub.0] y([sigma]) [??] [1/w] (s + [sigma])u([sigma]) d[sigma].

Finally a real inversion formula of the Laplace transform on reproducing kernel Hilbert spaces is stated as follows.

Theorem 4.3. Let [H.sub.[alpha]] be the solution to (4.2). If F [member of] [??]([H.sub.w]), then the following inversion formula of the Laplace transform holds.

[mathematical expression not reproducible]. (4.7)

We give a remark on the meaning of the regularized solution in the case [f.sup.[dagger]] [??] [H.sub.w] and [g.sup.[dagger]] (s) = s [??][[f.sup.[dagger]]](s) [member of] [L.sup.2.sub.u], for which there exists the Tikhonov regularized solution [f.sup.[dagger].sub.[alpha]] = [L.sup.-1.sub.[alpha]] [g.sup.[dagger]] [member of] [H.sub.w] with [L.sub.[alpha]] = [alpha]I + L*L. This [f.sup.[dagger].sub.[alpha]] [member of] [H.sub.w] gives an approximation of [f.sup.[dagger]] [??] [H.sub.w] in the sense that it attains the minimum of the Tikhonov functional.

5 Weight Functions

In the proposed real inversion scheme, setting of function spaces is quite important We give three useful pairs of weight functions w and u satisfying the condition (4.1) If the growth rate of f is known a priori, then we can choose suitable spaces.

Example 5.1 ([8, 11]). w(t) = [e.sup.t]/t,u(s) = 1. For the pair, I = 1/2 and L[K.sub.t](s) = 1/[(s+1).sup.2]{l - [e.sup.-t(s+1)(t(s + 1) + 1)}, [??] [1/w] (s) = 1/[(s+1).sup.2].

We note that for any non-zero polynomial p [??] 0, we have p [??] [H.sub.w].

Example 5.2. w(t) = [(t +1).sup.-2n], u(s) = [e.sup.-s-1/s]. For the pair, I [approximately equal to] 0.38 (n = 1), 6.39 (n = 2), 828.95 (n = 3), 458661.07 (n = 4) and

L[K.sub.t](s) = (2n)!//[s.sup.2n+1]{[e.sub.2n](s) - [e.sub.2n](s(1+t))[e.sup.-st]), [??][1/w] (s) = (2n)!/[s.sup.2n+1][e.sub.2n](s), where

[e.sub.k](s) = 1 + s + [s.sup.2]/2! + ... + [s.sup.k]/k!.

We note that for any polynomial p with deg p [less than or equal to] n and p(0) = 0, we have p[member of] [H.sub.w].

Example 5.3. w(t) = [e.sup.-[square root of l], u(s) - [e.sup.-s-1/s]. For the pair, I [approximately equal to] 0.2463 and

L[K.sub.t](s) = 1/s[1-[e.sup.-st+[square root of t] + 1/[square root of s][e.sup.1/4s]{Ert(1/2[square root of s]) + Erf([square root of st] - 1/2[square root of s])}], [??][1/w] (s) = 1/s{1 + 1/[square root of s][e.sup.1/4s] Erfc(-1/2[square root of s]))},

Numerical Real Inversion of the Laplace Transform where

[mathematical expression not reproducible]

We note that for any polynomial p with p(0) = 0, we have p [member of] [H.sub.w].

6 Relaxing the Condition of [H.sub.w] at the Origin

Since the space [H.sub.w] requires a rigid condition at the origin, Theorem 4.3 is quite limited in application. We relax it by introducing an additional approximation with the mollifier.

Let p [member of] C[0, [infinity]) with supp p [subset] [0,2]. For a positive number [??] > 0, we let

P[??](x) = 1/[??][rho](x/[??]) x [greater than or equal to] 0

and

[f.sub.[??]](t) = [[integral].sup.t.sub.0][rho][??](s)f (t - s)ds, f [member of] [L.sup.[infinity]](0, [infinity]).

Then we have [f.sub.[??]] [member of] C[0, [infinity]), [f.sub.[??]](0) = 0 and L [f.sub.[??]] = L[f] * L[[rho][??]]. Moreover if f is continuous at t = [t.sub.0] > 0,

[mathematical expression not reproducible] (6.1)

For a given Laplace image F = L f, we let [F.sub.[??]] = F * L [rho][??] = L [f.sub.[??]].*If [f.sub.[??]] [member of] [H.sub.w], then the inversion [L.sup.-1] [F.sub.[??]] = [f.sub.[??]] is obtained from Theorem 4.3. The equation (6.1) indicates that f gives an approximation to f with small [??].

Finally we state an effctive real inversion of the Laplace transform. Suppose that F [member of] [L.sup.2.sub.u] is given as an image function of the Laplace transform. Let [H.sub.[alpha]] be the unique solution to (4.2). Then an approximate real inversion to F is given by

[f.sub.[alpha],[??]](t)= [L.sup.-1.sub.[alpha],[??]]F(t)= [[[integral].sup.[infinity].sub.0]. sF(s) L[[rho][??]](s)[H.sub.[alpha]](s,t)u(s) ds. (6.2)

An example of [rho] is given by

[mathematical expression not reproducible]

For this [rho], we have

L[[rho][??]](S) = 1/[[??].sup.2][S.sup.2][([e.sup.-[??]s]-l).sup.2].

We use this mollifier in our numerical examples given in the next section.

7 Numerical Examples

The rest of this paper presents some numerical examples to show effctiveness of the proposed inversion formula (6.2) in both typical and practical problems.

The following computations are processed with 200 decimal digits or 600 decimal digits precision with multiple-precision arithmetic environment 'exflib' . As stated before, a regularized equation is well-posed in the corresponding function spaces, and thus its numerical computations are expected to be stable. On the contrary, in actual numerical processing of the regularized equation, computational errors may grow rapidly and computation fails due to numerical instability, especially taking extremely small parameters for accurate approximation. Multiple-precision arithmetic enables us to realize computations virtually free from rounding errors by using enough precision in representations and arithmetic of real numbers, hence direct treatments of numerically unstable problems are achieved on digital computers. To reduce discretization errors, the integral equation (4.2) is discretized by the double exponential numerical quadrature rule and the collocation method with 6000 collocation points.

7.1 Typical Examples

The sixteen examples in Table 1 given in the review paper  are examined as typical problems.

Firstly we consider the choice of weight functions in [F.sub.6] = L [f.sub.6]. Saitoh's original framework  corresponds to w(t) = [e.sup.t]/t and u(s) = 1, and [f.sub.6] [??] [H.sub.w] and [F.sub.6] [??] [L.sup.2.sub.u]\(0, [infinity]). Then the numerical inversion shown in Fig. 1(a) oscillates with large amplitude. On the other hand, if we adopt w(t) = [(t+1).sup.-2] and u(s) = exp(-s-1/s), then [F.sub.6] [member of] [L.sup.2.sub.u](0, [infinity]) and numerical results shown in Fig. 1(b) shows good agreement with the exact solution [f.sub.6](t) = t.

Secondly, we examine the validity of the mollifier in [F.sub.5] = L [f.sub.5]. Numerical results by (4.7) without the mollifier are shown in Fig. 2(a). On the other hand, results by (6.2) with [??] = 0.1 and 0.01 are shown in Fig. 2(b). Comparing them, it is clear that the use of the moillifier is efficient to obtain stable and accurate numerical solutions even though the original function [f.sub.5](t) = 1 does not belong to [H.sub.w].

Next, the relation between the regularization parameter [alpha] and discontinuities of the original function is examined in [F.sub.12] = L [f.sub.12]. Computations with [alpha] > [10.sup.-12] are shown in Fig. 3(a), and those with [alpha] = [10.sup.-100] and [10.sup.-400] are shown in Fig. 3(b). The latter shows good agreement with the original functions including singularities. It implies the necessity of multiple-precision arithmetic for the sake of accurate approximation with small regularization parameters.

Other numerical examples in Table 1 are shown in Fig. 4 (the original functions are entire) and Fig. 5 (the original functions each have a singularity).

7.2 Applications to Practical Examples

In Fig. 6 we show several results of the proposed scheme (6.2) applied to inverse Laplace transform arising from certain practical problems in engineering, mechanics, and finance. The weights and the mollifying function are the same as the preceeding examples.

Example 7.1. In a circuit problem, the image function

F(S) = 1/s(s+c)(1/2sh - [e.sup.-2sh]/1-[e.sup.-2sh])

appears , where c and h are fixed numbers.

Example 7.2. The image function

F(s) = 1/s exp(-r[square root of s(1 + s)/1 + cs] (7.2)

appears in , where c = 0.4 and r = 0.5.

Example 7.3. The problem to find the original function for

F (s) =(100*-l)sinh([square root of s]/2)/s(s sinh[square root of s] + [square root of s]cosh [square root of s])

appears in .

Example 7.4. In , the image function of the Laplace transform

F(s) = exp(-2[[psi].sub.s])/s, cosh[[psi].sub.s] = [square root of 1 + [s.sup.2] + [s.sup.4]/16] (74)

appears.

Example 7.5. The Laplace transform of the Gaussian distribution

F(s) = exp([lambda]/[mu]{1-[(1+[2[mu].sup.2]s/[lambda]).sup.1/2]}) (7.5)

appears in  with [micro] = 5 and [lambda] = [10.sup.6].

Example 7.6. The problem to find the original function for

F{s) = 1-r/s-r(1-[e.sup.-s]) (7.6)

appears in , where r = 0.7,0.8,0.9, and 0.95.

These numerical results show availability of the proposed algorithm with weighted function spaces, the mollifier, and multiple-precision arithmetic to numerical real inversions of the Laplace transform.

Acknowledgments. The authors wish to thank Professor Saburou Saitoh for his valuable comments. They also thank Professors Tsutomu Matsuura and Yoshihiro Sawano for fruitful discussions with them. This work is partially supported by JSPS KAKENHI Grant Number 26400198.

References

 A. M. Cohen, Numerical methods for Laplace transform inversion, Springer, 2007.

 B. Davies and B. Martin, Numerical inversion of the Laplace transform: a survey and comparison of methods, J. Comput. Phys., vol. 33, 1979, pp. 1-32.

 H. Fujiwara, exflib--multiple-precision arithmetic software library, http://www-an.acs.i.kyoto-u.ac.jp/~fuj iwara/exflib

 H. Fujiwara and N. Higashimori, Tikhonov regularization on reproducing kernel Hilbert spaces and its application to numerical computation of real inversion of the Laplace transform (in Japanese), RIMS Kokyuroku, no. 1638, 2009, pp. 137-145.

 P. D. Iseger, Numerical transform inversion using Gaussian quadrature, Probab. Eng. Inf. Sci., vol. 20, 2006, pp. 1-44.

 R. Kress, Linear integral equations (1st ed.), Springer-Verlag, 1989.

 N. W. McLachlan, Complex variable theory and transform calculus with technical applications (2nd Ed.), Cambridge University Press, 1953.

 T. Matsuura, A. Al-Shuaibi, H. Fujiwara, and S. Saitoh, Numerical real inversion of the Laplace transform by using a Fredholm integral equation of the second kind, J. Anal. Appl., vol. 5, 2007, pp. 123-136.

 M. J. P. Musgrave and J. Tasi, Shock waves in diatomic chains - I. linear analysis, J. Mech. Phys. Solids, vol. 24, 1976, pp. 19-42.

 T. Sakurai, Numerical inversion for Laplace transforms of functions with discontinuities, Adv. Appl. Prob., vol. 36, 2004, pp. 616-642.

 S. Saitoh, Approximate real inversion formulas of the Laplace transform, Far East J. Math. Sci., vol. 83, 2004, pp. 727-733.

 Y. Sawano, H. Fujiwara and S. Saitoh, Real inversion formula of the Laplace transform on weighted function spaces, Complex Anal. Oper. Theory, vol. 2, 2008, pp. 511-521.

 D. Sheen, I. H. Sloan and V. Thomee, A parallel method for time discretization of parabolic equations based on Laplace transformation and quadrature, IMA J. Num. Anal., vol. 32, 2002, pp. 269-299.

 G. Spada and L. Boschi, Using the Post-Widder formula to compute the Earth's viscoelastic Love numbers, Geophys. J. Int., vol. 166, 2006, pp. 309-321.

 R. I. Tanner, Note on the Rayleigh problem for a visco-elastic fluid, Z. Angew. Math. Phys., vol. 13, 1962, pp. 573-580.

 T. C. T. Ting and P. S. Symonds, Longitudinal impact on viscoplastic rods, J. Appl. Mech., vol. 31, 1964, pp. 199-207.

Hiroshi Fujiwara

Graduate School of Informatics, Kyoto University, Kyoto, 606-8501, Japan

E-mail: fujiwara@acs.i.kyoto-u.ac.jp

Center for the Promotion of Interdisciplinary Education and Research, Kyoto University, Kyoto, 606-8501, Japan

```Table 1: Test problems in 

[f.sub.1](t) = [J.sub.0](t)        [F.sub.1](s) =
[([s.sup.2] + 1).sup.-1/2]
[f.sub.2](t) = [([pi]t).sup.-1/2]  [mathematical expression
cos([2t.sup.1/2])                  not reproducible]
[f.sub.3](t) = [e.sup.-t/2]        [F.sub.3](s) = [(s + 1/2).sup.-1]
[f.sub.4](t) = [e.sup.-0.2t]       [F.sub.4](s) =
sin(t)                             [([(s + 0.2).sup.2] +1).sup.-1]
[f.sub.5](t) = 1                   [F.sub.5](s) = [s.sup.-1]
[f.sub.6](t) = t                   [F.sub.6](s) = [s.sup.-2]
[f.sub.7](t) = [te.sup.-t]         [F.sub.7](s) = [(s + 1).sup.-2]
[f.sub.8](t) = sin(t)              [F.sub.8](s) =
[([s.sup.2] + 1).sup.-1]
[f.sub.9](t) = [([pi]t).sup.-1/2]  [F.sub.9](s) = [s.sup.-1/2]
[f.sub.10](t) = H(t - 5)           [F.sub.10](s) = [s.sup.-1][e.sup.-5s]
[f.sub.11](t) = -[gamma] - ln(t)   [F.sub.11](s) = [s.sup.-1] ln(s)
[f.sub.12](t) = square wave        [F.sub.12](s) = [(s(1 +
[e.sup.-s])).sup.-1]
[f.sub.13](t) = t cos(t)           [F.sub.13](s) = ([s.sup.-2] - 1)
[([s.sup.-2] + 1).sup.-2]
[f.sub.14](t) = ([e.sup.-t/4]      [F.sub.14](s) = [(s + 1/2).sup.-1/2]
- [e.sup.-t/2])                    - [(s + 1/4).sup.1/2]
[([4[pi]t.sup.3]).sup.-1/2]
[f.sub.15](t) = 2[e.sup.-4/t]      [mathematical expression not
[([[pi]t.sup.3]).sup.-1/2]         reproducible]
[f.sub.16](t) =[t.sup.-1]sin(t)    [F.sub.16](s) = Arctan([s.sup.-1])
```
Author: Printer friendly Cite/link Email Feedback Fujiwara, Hiroshi; Higashimori, Nobuyuki Libertas Mathematica Report 1USA Nov 1, 2014 3881 To Professor Saburou Saitoh on the occasion of his 70th birthday. Exploring the spectra of some classes of paired singular integral operators: the scalar and matrix cases. Arithmetic functions Laplace transformation Laplace transforms