# A Semi-Analytical Solution of Inverse Laplace Transform

Authors: Shuang Luo; Fu-yao Zhao
Journal: Journal of Mathematics
 (2022)

---

## Abstract

We propose a general method for constructing the semi-analytical solution of the inverse Laplace transform, realized through the powerful exponential approximation invented by Wang et al. in 1993. Bearing their credits, this method inherits all the merits such as analytical expression, avoiding free parameters, simple calculation with high accuracy, and the availability of error estimation. Illustrating calculations indicate the potential applications to the vast problems in the fields of mathematical physics as well as engineering and medicine.


---

## Body

## 1. Introduction

Laplace transform, bearing rigorous theory, is a common tool for solving initial value problems in mathematical physics, demonstrating wide applications in engineering as well as medical and pharmic calculations. For a functionf (x) defined at x ≥ 0, Laplace transform and its inverse transform are defined as follows (see 5”.3.1 of [1]).(1)Ls=∫0+∞fxe−sxdx,(2)fx=L−1Ls=12πi∫σ−j∞σ+j∞Lsesxds,where σ is a sufficient large positive number, j2 = −1.In actual calculation,L (s) in equation (1) can be obtained by solving the transformed differential equation, but f (x) in equation (2) can hardly be expressed analytically in general. Sometimes special skills are needed (for example, reference [2]), but most of the time, only numerical inverse transform can be carried out. One of the most widely used numerical methods was proposed by Stehfest in [3, 4]. Nevertheless, for problems in mathematical physics, analytical solutions are generally more desirable and welcomed. As a compromise, a lot of efforts have been devoted to finding the semi-analytical solution of equation (2), posing the marriage between analytical expression and approximate expansion. While the specific calculation techniques vary, the basic ideas are similar. For example, many research studies attempt to decompose f (x) into an expansion of some known function ψi (x):(3)fx=∑i=0+∞ciψix,where ci are undetermined constants and can be determined by equation (1). This method greatly depends on the selection of ψi (x) and the properties of f (x), while the resulting inversion is not always stable (or the inversion is stable, but the errors from the real solution are not always small enough) [5]. Other researchers [6–8] have proposed several inverse transform methods based on Fourier series, yet these methods rely on two “free parameters” which are unknown before calculation, posing a great influence on the calculation results. Somewhat complicated, Fioravanti and Viano have proposed an inverse transform method based on Pollaczek polynomial [9], which can provide accuracy estimation for some functions.Based on the exponential approximation method, we attempt to provide an approximate analytic solution of the inverse Laplace transforms. The basic idea is, instead off (x), to decompose L (s) into an expansion of a certain “basis function,” and the inverse transform of this “basis function” has an analytic solution. Many expansion methods can analyze errors and convergence (e.g., reference [2]); therefore, the error of equation (2) may also be estimated in many cases.In this paper, exponential expansion is adopted. Now there are many different methods on exponential approximation [10]. The common ones include Prony method [11], Prony method based on generalized eigenvalues [12], nonlinear least square method [13], and so on. Based on the consideration of accuracy and simplicity, this paper will refer to the method proposed by Wang et al. [14]. The uniform convergence of this method is guaranteed, and the calculation of this method with high accuracy (<10−3) is simple. Details are presented in the next section.


## 2. Semi-Analytical Solutions of Inverse Laplace Transforms

### 2.1. Method

In this paper, we substitutes = p2 in equation (2) and define the function g (p) as follows:(4)gp=gs=sLs.Forg (p) being continuous in [0, +∞), we assume that an exponential approximation to equation (4) holds as follows.(5)gp≈∑r=0nCre−rp=∑r=0nCre−rs,where Cr are constants. According to the Laplace transform table (see section 81 in reference [15] or many textbooks on integral transforms), equations (2), (4), and (5) can be used to get the following equation:(6)fx≈∑r=0nCr12πi∫σ−j∞σ+j∞1se−rsesxds=∑r=0nCrerfcr2x.Equation (6) is the semi-analytical solution proposed in this paper. Obviously, equation (6) can be seen as some kind of expansion, where “erfc” is the complementary error function:(7)erfcx=2π∫x+∞e−η2dη.Obviously, the validity of this method depends on equation (5). Based on reference [14], assume that there is a continuous function g (x) ∈ [0, +∞), and g (x) is in existence when x ⟶ ∞; then, g (x) can be approximated by the function gn (x).(8)gnx=∑r=02n−1Cre−rx,where Cr can be determined by the following equation:(9)Cr=a02+∑k=1n−1kak+∑k=1n−1−1k+n1−knan+k,r=0,∑k=rn−1k−r2kk+rCk+rk−r22r−1ak+∑k=1n−1−1n+k−r1−kn2n+kn+k+rCn+k+rn+k−r22r−1an+k,1≤r≤n,∑k=r−nn−1−1n+k−r1−kn2n+kn+k+rCn+k+rn+k−r22r−1an+k,n+1≤r≤2n−1,where Cmn is the combinatorial number, and according to Gauss–Chebyshev quadrature formula, ak (k = 0, 1, 2, …, 2n − 1) can be calculated as follows:(10)ak=2π∫−11g−ln1+y2coskarccosy1−y2dt≈2N∑l=1Ng−ln1+yl2cosk2l−12Nπ,where yl is the zero point of Chebyshev polynomial of N-th order:(11)yl=cos2l−12Nπ,l=1,2,…,N.Usually,ak will be accurate enough when N = 20. Here, for r > 0, more concise combinatorial notations are adopted for the coefficients of ak and an+k (we have corrected minor clerical errors in the original formula which can be easily checked by readers).


### 2.2. Convergence and Accuracy

AssumeL (s) is to be known and equation (5) is true, and its maximum relative error αe can be found to satisfy(12)Ls−1s∑r=0nCre−rs≤αeLs.Note that we have not chosen a specific exponential approximation method at this point. Since Laplace transform is linear, we have(13)fx−∑r=0nCrerfcr2x≤αeL−1Ls.So, we have the following theorem about equation (6).Theorem 1.
The convergence and accuracy of equation (6) satisfy equation (13).Proof.
To perform the inverse Laplace transform to the term in the absolute value sign at the left side of equation (12) directly, we have(14)fx−∑r=0nCrerfcr2x=L−1Ls−1s∑r=0nCre−rs.
Take the absolute value of equation (14) and pay attention to the property of absolute value, and equation (13) could be founded.
According to Theorem1, the property of the function L (s), the convergence and accuracy of exponential approximation equation (5) directly affect the convergence and accuracy of equation (6). Usually, when researchers propose an exponential approximation method, they will give theorems related to convergence and accuracy. Therefore, the convergence and accuracy of equation (6) can be estimated in advance. For example, if equation (8) will be used, there is a theorem on convergence and error estimation about equation (8) in reference [14].Theorem 2.
Ifφ (t) = g (−ln [(1 + cost)/2]) ∈ LipMα and t ∈ [0, π], then there is |gn (x) − g (x)| ≤ 4παM/nα uniformly on x ∈ [0, +∞).The proof of Theorem2 can be found in Section 3 in reference [14].Theorems1 and 2 can be used for convergence analysis or error estimation. For example, when |L (s)| ≥ Mf > 0, ReLs⋅ImLs ≠ 0 and the condition for Theorem 2 is true; according to Theorem 1,(15)fx−∑r=0nCrerfcr2x≤4παMMfnαfx.Therefore, as long as the property of functionL (s) is good, equation (6) can have a good accuracy.


## 2.1. Method

In this paper, we substitutes = p2 in equation (2) and define the function g (p) as follows:(4)gp=gs=sLs.Forg (p) being continuous in [0, +∞), we assume that an exponential approximation to equation (4) holds as follows.(5)gp≈∑r=0nCre−rp=∑r=0nCre−rs,where Cr are constants. According to the Laplace transform table (see section 81 in reference [15] or many textbooks on integral transforms), equations (2), (4), and (5) can be used to get the following equation:(6)fx≈∑r=0nCr12πi∫σ−j∞σ+j∞1se−rsesxds=∑r=0nCrerfcr2x.Equation (6) is the semi-analytical solution proposed in this paper. Obviously, equation (6) can be seen as some kind of expansion, where “erfc” is the complementary error function:(7)erfcx=2π∫x+∞e−η2dη.Obviously, the validity of this method depends on equation (5). Based on reference [14], assume that there is a continuous function g (x) ∈ [0, +∞), and g (x) is in existence when x ⟶ ∞; then, g (x) can be approximated by the function gn (x).(8)gnx=∑r=02n−1Cre−rx,where Cr can be determined by the following equation:(9)Cr=a02+∑k=1n−1kak+∑k=1n−1−1k+n1−knan+k,r=0,∑k=rn−1k−r2kk+rCk+rk−r22r−1ak+∑k=1n−1−1n+k−r1−kn2n+kn+k+rCn+k+rn+k−r22r−1an+k,1≤r≤n,∑k=r−nn−1−1n+k−r1−kn2n+kn+k+rCn+k+rn+k−r22r−1an+k,n+1≤r≤2n−1,where Cmn is the combinatorial number, and according to Gauss–Chebyshev quadrature formula, ak (k = 0, 1, 2, …, 2n − 1) can be calculated as follows:(10)ak=2π∫−11g−ln1+y2coskarccosy1−y2dt≈2N∑l=1Ng−ln1+yl2cosk2l−12Nπ,where yl is the zero point of Chebyshev polynomial of N-th order:(11)yl=cos2l−12Nπ,l=1,2,…,N.Usually,ak will be accurate enough when N = 20. Here, for r > 0, more concise combinatorial notations are adopted for the coefficients of ak and an+k (we have corrected minor clerical errors in the original formula which can be easily checked by readers).


## 2.2. Convergence and Accuracy

AssumeL (s) is to be known and equation (5) is true, and its maximum relative error αe can be found to satisfy(12)Ls−1s∑r=0nCre−rs≤αeLs.Note that we have not chosen a specific exponential approximation method at this point. Since Laplace transform is linear, we have(13)fx−∑r=0nCrerfcr2x≤αeL−1Ls.So, we have the following theorem about equation (6).Theorem 1.
The convergence and accuracy of equation (6) satisfy equation (13).Proof.
To perform the inverse Laplace transform to the term in the absolute value sign at the left side of equation (12) directly, we have(14)fx−∑r=0nCrerfcr2x=L−1Ls−1s∑r=0nCre−rs.
Take the absolute value of equation (14) and pay attention to the property of absolute value, and equation (13) could be founded.
According to Theorem1, the property of the function L (s), the convergence and accuracy of exponential approximation equation (5) directly affect the convergence and accuracy of equation (6). Usually, when researchers propose an exponential approximation method, they will give theorems related to convergence and accuracy. Therefore, the convergence and accuracy of equation (6) can be estimated in advance. For example, if equation (8) will be used, there is a theorem on convergence and error estimation about equation (8) in reference [14].Theorem 2.
Ifφ (t) = g (−ln [(1 + cost)/2]) ∈ LipMα and t ∈ [0, π], then there is |gn (x) − g (x)| ≤ 4παM/nα uniformly on x ∈ [0, +∞).The proof of Theorem2 can be found in Section 3 in reference [14].Theorems1 and 2 can be used for convergence analysis or error estimation. For example, when |L (s)| ≥ Mf > 0, ReLs⋅ImLs ≠ 0 and the condition for Theorem 2 is true; according to Theorem 1,(15)fx−∑r=0nCrerfcr2x≤4παMMfnαfx.Therefore, as long as the property of functionL (s) is good, equation (6) can have a good accuracy.


## 3. Example

According to the first two sections, the solving method of the approximate analytical solution in this paper can be as follows: firstly, getL (s) in equation (1) by solving the differential equation of the original problem; secondly, get g (p) by equation (4); thirdly, get ak by equations (10) and (11); then get Cr by equation (9); finally, the approximate analytic solution of equation (6) can be obtained (the solving process is similar when using other exponential approximation method). A specific example will be given below.In geothermal energy extraction engineering, there are two important models: single fracture model and parallel fracture model. Single fracture model assumes that there is an infinitely long linear fracture with a width of 2b in the infinite rock plane, and the water with the initial temperature Tin flows into the fracture at a constant velocity of vw and has heat exchange with rock. Parallel fracture model is similar, but it assumes that there are an infinite number of parallel fractures, and the distance between fractures is 2D. Assuming that the temperature of water is Tw and taking time t as the target variable, the solutions corresponding to differential equations of these two models, after the Laplace transforms to Tw, are [16](16)L1x,s=Tin−T0sexp−xvws+λrbρwcw1αr/s+λr/h+T0s,(17)L2x,s=Tin−T0sexp−xvws+λrbρwcwρrcrsλrtanhρrcrsλrD+T0s,where L1 corresponds to the single fracture model and L2 corresponds to the parallel fracture model. The meanings and values of these parameters in this example are as follows: Tin—initial temperature of water, 40°C; T0—initial temperature of rock, 120°C; b—half of the fracture width, 0.1 mm; vw—flow velocity of water, 0.1 m/s; λr—thermal conductivity of rock, 3 W/(m·K); ρr—density of rock, 2700 kg/m3; cr—specific heat capacity of rock, 800 J/(kg·K); ρw—density of water, 1000 kg/m3; cw—specific heat capacity of water, 4200 J/(kg·K); αr = λr/(ρrcr), thermal diffusion coefficient; x—distance of observation point, 100 m in this example; h—convective heat transfer coefficient in the single fracture model, 10∼100 W/(m2·K); and D—half of the fracture spacing in the parallel fracture model, 3∼30 m.Substitute equations (16) and (17) into (6) and consider the delay theorem of Laplace transform (see reference [15] or many textbooks on integral transforms), and the result is the temperature of water Tw (x, t):(18)Tw1x,t=Tin−T0Ht−xvw∑i=02n−1Cierfciα2t−x/vw+T0,Tw2x,t=Tin−T0Ht−xvw∑i=02n−1Ci′erfciβ2t−x/vw+T0,where H is the Heaviside step function and α and β are as follows (in order to obtain higher accuracy, equations (16) and (17), respectively, adopt αs and βs as independent variables in exponential approximation):(19)α=xλrvwbρwcwαr,β=xλrvwbρwcwρrcrλr.In addition, we compare this approximate analytical solution with the finite element method (FEM) and Stehfest method. The accuracy of finite element method is very high and can be used as a benchmark. The formula of the Stehfest method is as follows:(20)ft≈ln2t∑i=1NViLln2ti,where(21)Vi=−1N/2+i∑k=inti/2+1/2mini,N/2kN/22k!N/2−k!k!k−1!i−k!2k−i!.The calculation results of the three methods are compared in Figures1(a), 1(b), 2(a), and 2(b).Figure 1
(a) Single fracture model (h = 10 W/(m2·K)). (b) Single fracture model (h = 100 W/(m2·K)).
(a)(b)Figure 2
(a) Parallel fracture model (D = 3 m). (b) Parallel fracture model (D = 30 m).
(a)(b)As can be seen from the figures, the method in this paper can have a very good accuracy. This example can provide reference for geothermal energy extraction engineering.


## 4. Discussion

The method in this study has a limitation. Obviously, equation (6) always tends to a finite value when x ⟶ +∞, so it can only be used on the function which has a finite limit value at infinity (such as problems where the Laplace function has only algebraic decay) and not applicable to functions which oscillate at infinity (such as sine function). But the unsolved functions in many practical problems tend to be finite values at infinity. According to computational experience, when f (x) tends to a finite value at infinity, the influence of Hadamard ill-conditioned of the inverse Laplace transform [17] can be neglected on equation (6).Besides, for the inverse Laplace change of a specific problem, the exact answer is usually unknown. So, the error is unknown in fact (especially for methods that have free parameters). But ifL (s) is known, the accuracy of equation (5) can be directly compared. This is useful when the exact solution is unknown. For many engineering problems [18–23], medicine problems [24], and integro-differential equations of second-order problems [25–27], this method might be useful.



---

*Source: 9129727-2022-08-19.xml*