A Semi-Analytical Solution of Inverse Laplace Transform

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.

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 function f (x) defined at x ≥ 0, Laplace transform and its inverse transform are defined as follows (see 5”.3.1 of [1]).where σ is a sufficient large positive number, j² = −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):where c_i 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 of f (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⁻³) is simple. Details are presented in the next section.

2. Semi-Analytical Solutions of Inverse Laplace Transforms

2.1. Method

In this paper, we substitute s = p² in equation (2) and define the function (p) as follows:

For (p) being continuous in [0, +∞), we assume that an exponential approximation to equation (4) holds as follows.where C_r 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:

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:

Obviously, the validity of this method depends on equation (5). Based on reference [14], assume that there is a continuous function (x)  [0, +∞), and (x) is in existence when x ⟶ ∞; then, (x) can be approximated by the function (x).where C_r can be determined by the following equation:where is the combinatorial number, and according to Gauss–Chebyshev quadrature formula, a_k (k = 0, 1, 2, …, 2n − 1) can be calculated as follows:where y_l is the zero point of Chebyshev polynomial of N-th order:

Usually, a_k will be accurate enough when N = 20. Here, for r > 0, more concise combinatorial notations are adopted for the coefficients of a_k and a_n+k (we have corrected minor clerical errors in the original formula which can be easily checked by readers).

2.2. Convergence and Accuracy

Assume L (s) is to be known and equation (5) is true, and its maximum relative error α_e can be found to satisfy

Note that we have not chosen a specific exponential approximation method at this point. Since Laplace transform is linear, we have

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 haveTake the absolute value of equation (14) and pay attention to the property of absolute value, and equation (13) could be founded.
According to Theorem 1, 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) =  (−ln [(1 + cost)/2])  Lip_Mα and t  [0, π], then there is | (x) −  (x)| ≤ 4π^αM/n^α uniformly on x  [0, +∞).

The proof of Theorem 2 can be found in Section 3 in reference [14].

Theorems 1 and 2 can be used for convergence analysis or error estimation. For example, when |L (s)| ≥ M_f > 0,  ≠ 0 and the condition for Theorem 2 is true; according to Theorem 1,

Therefore, as long as the property of function L (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, get L (s) in equation (1) by solving the differential equation of the original problem; secondly, get (p) by equation (4); thirdly, get a_k by equations (10) and (11); then get C_r 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 T_in flows into the fracture at a constant velocity of 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 and taking time t as the target variable, the solutions corresponding to differential equations of these two models, after the Laplace transforms to , are [16]where L₁ corresponds to the single fracture model and L₂ corresponds to the parallel fracture model. The meanings and values of these parameters in this example are as follows: T_in—initial temperature of water, 40°C; T₀—initial temperature of rock, 120°C; b—half of the fracture width, 0.1 mm; —flow velocity of water, 0.1 m/s; λ_r—thermal conductivity of rock, 3 W/(m·K); ρ_r—density of rock, 2700 kg/m³; c_r—specific heat capacity of rock, 800 J/(kg·K); —density of water, 1000 kg/m³; —specific heat capacity of water, 4200 J/(kg·K); α_r = λ_r/(ρ_rc_r), 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/(m²·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 (x, t):where H is the Heaviside step function and α and β are as follows (in order to obtain higher accuracy, equations (16) and (17), respectively, adopt and as independent variables in exponential approximation):

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:where

The calculation results of the three methods are compared in Figures 1(a), 1(b), 2(a), and 2(b).

(a)

(b)

Figure 1 

(a) Single fracture model (h = 10 W/(m²·K)). (b) Single fracture model (h = 100 W/(m²·K)).

(a)

(b)

Figure 2 

(a) Parallel fracture model (D = 3 m). (b) Parallel fracture model (D = 30 m).

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 if L (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.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Shuang Luo applied this method to a specific function and provided the computational code. Fu-yao Zhao generalized this method as a general method and wrote the manuscript.

Acknowledgments

We sincerely appreciate the work in reference [14] for the inspiration it has given us although we never had any chance to communicate with its authors. If the method in this paper really means something, we would like to call it LZW method (W represents the first author of reference [14]). We also sincerely appreciate Associate Professor Yiwei Li who made a lot of modifications to the original manuscript and gave us many valuable suggestions. This study was financially supported by PhD Scientific Research Initial Funding of Changzhi Medical College (no. BS202026).