Numerical Inverse Laplace Transform Methods for Advection-Diffusion Problems

Abstract

Partial differential equations arising in engineering and other sciences describe nature adequately in terms of symmetry properties. This article develops a numerical method based on the Laplace transform and the numerical inverse Laplace transform for numerical modeling of diffusion problems. This method transforms the time-dependent problem to a corresponding time-independent inhomogeneous problem by employing the Laplace transform. Then a local radial basis functions method is employed to solve the transformed problem in the Laplace domain. The main feature of the local radial basis functions method is the collocation on overlapping sub-domains of influence instead of on the whole domain, which reduces the size of the collocation matrix; hence, the problem of ill-conditioning in global radial basis functions is resolved. The Laplace transform is used in comparison with a finite difference technique to deal with the time derivative and avoid the effect of the time step on numerical stability and accuracy. However, using the Laplace transform sometimes leads to a solution in the Laplace domain that cannot be converted back into the real domain by analytic methods. Therefore, in such a case, the Laplace transform is inverted numerically. In this investigation, two inversion techniques are utilized; (i) the contour integration method, and (ii) the Stehfest method. Three test problems are used to evaluate the proposed numerical method. The numerical results demonstrate that the proposed method is computationally efficient and highly accurate.

1. Introduction

Advection-diffusion equations (ADEs) have an important role in many physical systems, i.e., vorticity, heat, energy, mass and fluid flow [1]. Symmetry is a fundamental property of nature and its phenomena; in physics, it ensures the existence of conservation laws, and, in mathematics, it means a transformation from one solution to another [2]. ADEs are able to adequately describe physical and biological processes and have symmetry properties [3]. Furthermore, many physical phenomena relating to the transformation of energy during advection and diffusion processes can be described by ADEs [4]. It is a $1\mathrm{D}$ PDE which demonstrates the advection and diffusion of mass, vorticity, heat and energy [5], etc. ADEs can be utilized to describe heat transfer in draining film [6], the dispersion of traces in porous media [7], the transfer of pollutants in the atmosphere [8], and water transmission in soil [9]. In this investigation, we consider an ADE of the form [10,11]:   

$$\frac{\partial C(x,t)}{\partial t}+\overline{\theta }\frac{\partial C(x,t)}{\partial x}-\overline{\alpha }\frac{{\partial }^{2}C(x,t)}{\partial x^2}=h(x,t),x\in \mathsf{\Omega },t>0,$$

(1)

with initial condition

$$C(x,0)=C_0\left(x\right),x\in \mathsf{\Omega },$$

and Dirichlet boundary conditions

$${\mathcal{L}}_{B}C(x,t)=g(x,t),x\in \partial \mathsf{\Omega },$$

(2)

where $\mathsf{\Omega }$ is the domain and $\partial \mathsf{\Omega }$ is its boundary, $\mathcal{L}=\frac{\partial }{\partial x}+\frac{{\partial }^2}{\partial x^2}$ and ${\mathcal{L}}_B$ are the linear differential and boundary operators; $\overline{\theta }$ and $\overline{\alpha }$ are the advection and diffusion coefficients, respectively, where $\overline{\theta }$ is an arbitrary constant and $\overline{\alpha }$ is a positive constant. Whereas, $C(x,t)$ is the concentration and $h(x,t)$, $g(x,t)$ are given smooth functions.

ADEs have been studied by many authors, and various approximate and analytic techniques have been developed for the solution of $1\mathrm{D}$ ADEs [12]. Kumar et al. [13,14] studied the analytic solution of $1\mathrm{D}$ ADEs with variable coefficients. Zoppou and Knight [15] investigated the analytic solution of an ADE with variable coefficients. In [16], the analytic solution of 2D ADE, arising in cytosolic calcium concentration distribution, was examined. Further information on other analytic methods for ADEs can be found in [17,18,19] and the references therein.

However, in many situations, the analytic solution of ADEs is difficult to obtain. Numerical techniques are then used to approximate the solution of ADEs. The numerical investigation of the solutions of ADEs has been considered by many researchers, and numerous methods have been proposed. These include the stable explicit finite difference method (FDM) [20], the implicit FDM [21], the compact FDM [10], the cubic trigonometric B-splines method [11], the finite element method (FEM) [22], and the boundary element method [23]. The numerical solutions of ADEs are difficult to obtain due to: (i) the first- and second-order derivatives with respect to space variables. Depending on the values of $\overline{\alpha }$ and $\overline{\theta }$, the equation becomes parabolic for diffusion-dominated and hyperbolic for advection-dominated processes. Traditionally used FDMs face problems when oscillations and smoothing of the wave front are introduced [24,25]; (ii) since all the methods depend on a mesh, it is necessary to construct a fine mesh for optimal results. However, mesh generation is the most time-consuming part of the solution process and, for complex geometries, problems can occur with implementation.

Due to the complexity of mesh generation, remarkable work has recently been undertaken on the development of mesh-free or mesh-less methods. Numerous mesh-free methods for approximating the solutions of PDEs in different fields of engineering and other fields have been developed. The major advantages of mesh-less methods are: (i) they do not require a mesh, which can be challenging in 3D-space cases; (ii) they are more suitable than mesh-based methods in cases of large deformations or moving discontinuities. One of the common features of all mesh-less methods is their ability to construct a functional interpolation or approximation entirely from information about a set of dispersed nodes between which there is no predetermined connectivity or relationship. Various mesh-less methods have been developed to date. These include, the Galerkin mesh-free methods [26], the boundary knot method [27], the singular boundary method [28], and the local point interpolation method [29].

Another group of mesh-less methods that are based on radial basis functions (RBFs), represents one of the best tools for obtaining the numerical solutions to various real world problems [30,31,32,33,34]. The main characteristics of RBFs are their smoothness, spectral convergence and ease of implementation. RBF-based mesh-less methods are divided into two categories: global mesh-less methods (GRBFMs) and local mesh-less methods (LRBFMs). In GRBFMs, the solution is approximated using the linear combinations of RBFs; these methods involve full system matrices that are often ill-conditioned. Due to the involvement of full system matrices the approximation of large-scale problems becomes difficult for some RBF methods. The condition number of the system matrices increases as the number of interpolation nodes increases. Another factor is the use of infinitely smooth basis functions with arbitrary shape parameters [35]. The accuracy of the approximation and the conditioning of the system matrix depends on the value of the shape parameter. For GRBFMs, how to select an optimal shape parameter is still an open question [36].

In 2003, Shu et al. [37] proposed LRBFMs to solve $2\mathrm{D}$ Navier–Stokes equations. In LRBFMs, the collocation is made in overlapping sub-domains, which reduces the size of the system matrix, and then these small system matrices are solved for each node. Subsequently, LRBFMs have been applied to numerous problems in engineering and other science subjects, such as compressible flow problems [38], incompressible flow problems [39], arbitrarily shaped membranes [40], boundary value problems [41], and diffusion problems [42]. The main drawback of LRBFMs is that this method does not work for elliptic problems in a straightforward way. In addition, these LRBFMs are employed to approximate the solutions of time-dependent PDEs coupled with a time-stepping scheme. The drawback of the time-stepping technique is that it does not always lead to a stable solution. The finite difference time-stepping scheme is stable if the errors remain constant or decay during computation. In addition, in the finite difference time-stepping scheme, the accuracy is achieved at a very small time step and, hence, this scheme encounters an exponential increase in computing costs. To overcome this drawback of the time-stepping technique, the Laplace transform (LT) is used.

Some remarkable work has been published in the literature on the combination of LT with different space-discretization methods to avoid time-instability problems. Examples of this approach include the boundary particle method [43], Kansa’s method [44], the LRBFM on unit sphere method [45], the FDM [46], the spectral method [47], and the FEM [48,49]. In this study, we propose LRBFMs coupled with LT for the approximation of the solution of $1\mathrm{D}$ linear ADEs. For the numerical inversion, we use two well-known methods: (i) the contour integration method; and (ii) Stehfest’s method.

2. Implementation of the Method

In our numerical scheme, there are three main steps: (a) implementation of the Laplace transform, (b) implementation of a local RBF method, and (c) numerical inversion of the Laplace transform.

2.1. (a): Laplace Transform (LT)

The LT transforms the time-dependent advection diffusion equation from a time domain to an inhomogeneous time-independent problem to an LT domain.

Let $C(x,t)$, $t>0$ be a piecewise continuous, and let $\widehat{C}(x,s)$ be its LT, which is defined as

$$\widehat{C}(x,s)=\mathcal{L}\left\{C(x,t)\right\}={\int }_0^{\infty }{e}^{-st}C(x,t)dt,$$

(3)

where s denotes the LT parameter. The LT of $\frac{dC(x,t)}{dt}$ is given as

$$\mathcal{L}\left\{\frac{dC(x,t)}{dt}\right\}=s\widehat{C}(x,s)-C_0$$

(4)

Employing the LT to (1) and (2), we get

$$s\widehat{C}(x,s)-C_0-\mathcal{L}\widehat{C}(x,s)=\widehat{h}(x,s),$$

(5)

and

$$\mathcal{B}\widehat{C}(x,s)=\widehat{g}(x,s).$$

(6)

Simplifying (5) and (6), we get

$$(sI-\mathcal{L})\widehat{C}(x,s)=\widehat{G_1}(x,s),$$

(7)

$${\mathcal{L}}_B\widehat{C}(x,s)=\widehat{g}(x,s),$$

(8)

where

$$\widehat{G_1}(x,s)=C_0+\widehat{h}(x,s),$$

where I is the identity operator, $\mathcal{L}$ is the linear differential operator and ${\mathcal{L}}_B$ is the boundary differential operator. To solve the system (7) and (8), first, we discretize the operators $\mathcal{L}$ and ${\mathcal{L}}_B$ using LRBFM, then the system (7) and (8) is solved for each s in the LT domain. Finally, the solution of the problem can be obtained via inverse LT.

2.2. (b): Local RBF Method

Here, we propose a local RBF method for approximating the $1\mathrm{D}$ ADEs in Laplace space. This method provides a linear system that is sparse and well-conditioned. For the set $\{\widehat{C}\left(x_i\right):i=1,2,3,...,N\}$, where $\{x_i:i=1,2,3,...,N\}\subset \mathsf{\Omega }$, the local RBF approximation of $\widehat{C}\left(x\right)$ is of the form,

$$\widehat{C}\left(x_i\right)=\sum _{x_j\in {\mathsf{\Omega }}_i}{\alpha }_j^i\phi (\parallel x_i-x_j^i\parallel ),$$

(9)

where ${\mathit{\alpha }}^i={\left\{{\alpha }_j^i\right\}}_{j=1}^n$,  $\phi \left(r\right),r>0$ is the radial kernel. ${\mathsf{\Omega }}_i$ is the sub-domain which contains the node $x_i$ and its neighboring n nodes around it. Thus, we get N number linear systems of order $n\times n$ as

$$\left(\begin{array}{c}\widehat{C}\left(x_1^i\right)\\ \widehat{C}\left(x_2^i\right)\\ .\\ .\\ .\\ \widehat{C}\left(x_n^i\right)\end{array}\right)=\left(\begin{array}{cccccc}{\phi }_{11}& {\phi }_{12}& .& .& .& {\phi }_{1n}\\ {\phi }_{21}& {\phi }_{22}& .& .& .& {\phi }_{2n}\\ .& .& & .& & .\\ .& .& & .& & .\\ .& .& & .& & .\\ {\phi }_{n1}& {\phi }_{n2}& .& .& .& {\phi }_{nn}\end{array}\right)\left(\begin{array}{c}{\alpha }_1^i\\ {\alpha }_2^i\\ .\\ .\\ .\\ {\alpha }_n^i\end{array}\right),i=1,2,3,...,N,$$

(10)

which can be written as,

$${\widehat{\mathbf{C}}}^i={\mathsf{\Phi }}^i{\mathit{\alpha }}^i,i=1,2,3,...,N,$$

(11)

the matrix ${\mathsf{\Phi }}^i$ has entries $d_{kj}^i=\phi (\parallel x_k^i-x_j^i\parallel ),\mathrm{where}{x}_k^i,x_j^i\in {\mathsf{\Omega }}_i$; each of the above systems is then solved for the ${\mathit{\alpha }}^i={\left\{{\alpha }_j^i\right\}}_{j=1}^n$. The index i denotes that the nodes belong to the local domain ${\mathsf{\Omega }}_i$ of each node $x_i.$ In addition, for $\mathcal{L}$ we have,

$$\mathcal{L}\widehat{C}\left(x_i\right)=\sum _{x_j\in {\mathsf{\Omega }}_i}{\mathit{\alpha }}_j^i\mathcal{L}\phi (\parallel x_i-x_j^i\parallel ),$$

(12)

From Equation (12), we have

$$\mathcal{L}\widehat{C}\left(x_i\right)={\mathbf{\nu }}^i·{\mathit{\alpha }}^i,$$

(13)

where ${\mathit{\alpha }}^i$ is a column vector of order n and ${\mathbf{\nu }}^i$ is a row vector of order n having elements of the form

$${\mathbf{\nu }}^i=\mathcal{L}\phi (\parallel x_i-x_j^i\parallel ),x_j^i\in {\mathsf{\Omega }}_i,$$

(14)

From Equation (11), we have

$${\mathit{\alpha }}^i={\left({\mathsf{\Phi }}^i\right)}^{-1}{\widehat{\mathbf{C}}}^i,$$

(15)

From (13) and (15), we have

$$\begin{array}{cc}\hfill \mathcal{L}\widehat{C}\left(x_i\right)=& {\mathbf{\nu }}^i{\left({\mathsf{\Phi }}^i\right)}^{-1}{\widehat{\mathbf{C}}}^i\hfill \\ \hfill =& {\Theta }^i{\widehat{\mathbf{C}}}^i,\hfill \end{array}$$

(16)

where,

$${\Theta }^i={\mathbf{\nu }}^i{\left({\mathsf{\Phi }}^i\right)}^{-1},$$

(17)

thus, the operator $\mathcal{L}$ has an approximation via local radial basis functions at each node $x_i$ as

$$\mathcal{L}\widehat{C}\equiv \mathbf{D}\widehat{\mathit{C}},$$

(18)

where D denotes the sparse matrix of order N which approximates the linear differential operator $\mathcal{L}$. A similar procedure can be applied to the boundary operator ${\mathcal{L}}_B$ as

$${\mathcal{L}}_B\widehat{C}\equiv \mathbf{B}\widehat{\mathit{C}},$$

(19)

substituting Equations (18) and (19) in Equations (7) and (8), we get

$$(sI-\mathbf{D})\widehat{C}(x,s)=\widehat{G_1}(x,s),$$

(20)

$$\mathbf{B}\widehat{C}(x,s)=\widehat{g}(x,s),$$

(21)

Solving the system (20) and (21) for each point s, we will obtain the approximate solution $\widehat{C}(x,s).$ Finally, we obtain the solution $C(x,t)$ of the problem defined in Equations (1) and (2) using the numerical inverse LT.

2.3. (c): Inverse Laplace Transform

This section is about numerical inverse LT methods. In this investigation, we use two numerical inverse LT methods: (i) the contour integration method; and (ii) Stehfest’s method.

(i) Contour Integration Method

In the contour integration method, the solution $C(x,t)$ is obtained via the inverse LT given as

$$C(x,t)=\frac{1}{2\pi \iota }{\int }_{\zeta -\iota \infty }^{\zeta +\iota \infty }{e}^{st}\widehat{C}(x,s)ds=\frac{1}{2\pi \iota }{\int }_{\mathsf{\Gamma }}{e}^{st}\widehat{C}(x,s)ds,t>0,$$

(22)

where $\mathsf{\Gamma }$ is an appropriately selected contour of integration in a complex plane. The solution of the problem (1) and (2) mainly depends on $\mathsf{\Gamma }$. For optimal solution of (22), an optimal contour parabolic or hyperbolic may be utilized. Various contour integration methods are developed for the approximation of the inverse Laplace transform. However, in this study, we utilize two recently proposed contours in [47,48]. The parabolic contour is given as

$$s\left(\upsilon \right)=\eta {(1+\iota \upsilon )}^2\left({\mathsf{\Gamma }}_P\right),$$

(23)

where $\eta$ is an unknown parameter. In [48], the authors proposed a contour given as

$$s\left(\upsilon \right)=\chi +\epsilon -\chi sin(\varsigma -\iota \upsilon ),\mathrm{for}-\infty <\upsilon <\infty ,\left({\mathsf{\Gamma }}_h\right)$$

(24)

where $\chi >0,$ $\frac{\pi }{2}<{\alpha }_1<\pi ,$ $0<\varsigma <{\alpha }_1-\frac{\pi }{2}$, $\epsilon \ge 0.$

Using Equations (23) and (24) in Equation (22), we have

$$C(x,t)=\frac{1}{2\pi \iota }{\int }_{-\infty }^{\infty }{e}^{s\left(\upsilon \right)t}\widehat{C}(x,s\left(\upsilon \right))s^{\prime }\left(\upsilon \right)d\upsilon$$

(25)

Equation (25) is solved via the trapezoidal rule, with step k as

$$C_{App}(x,t)=\frac{1}{2\pi i}\sum _{j=-M_p}^{M_p}k{e}^{s\left({\upsilon }_j\right)t}\widehat{C}(x,s\left({\upsilon }_j\right))\acute{s}\left({\upsilon }_j\right),{\upsilon }_j=jk.$$

(26)

2.4. Error Analysis

2.4.1. Parabolic Contour

Here, we discuss the error analysis of our scheme. We consider the integral

$$I={\int }_{-\infty }^{\infty }G\left(\vartheta \right)d\vartheta ,$$

and its infinite and finite approximations

$$I_k=k\sum _{k=-\infty }^{\infty }G\left(jk\right),I_{k;M_p}=k\sum _{j=-M_p}^{M_p}G\left(jk\right).$$

The quantities $E_T=|I_{k;M_p}-I_k|$ and $E_d=|I-I_k|$ are called the truncation and discretization errors. The proof of the discretization error of the proposed numerical scheme is based on the following theorem:

Theorem 1

([47]). Let $s=\xi +\iota \rho$, where $\xi ,\rho \in \mathbb{R}.$ Let $G\left(s\right)$ be analytic in the strip $-a_1<\rho <a_2,$ for some $a_1,a_2>0,$ with $G\left(s\right)\to 0,$ as $\left|s\right|\to \infty$. Let for $Q_{+},Q_{-},$ $G\left(s\right)$ satisfies

$${\int }_{-\infty }^{\infty }\left|G(\xi +\iota r)\right|d\xi \le Q_{+},$$

(27)

$${\int }_{-\infty }^{\infty }\left|G(\xi -\iota z)\right|d\xi \le Q_{-},$$

(28)

$\forall 0<r<a_1,0<z<a_2.$ Then

$$|I_k-I|\le \frac{M_p(+)}{exp\left(\frac{2\pi a_1}{k}\right)-1}+\frac{M_p(-)}{exp\left(\frac{2\pi a_2}{k}\right)-1}.$$

Remark 1.

In the literature, G is typically a real valued function, $a_1=a_2,Q_{+}=Q_{-}=Q,$ say, then we obtain an estimate given as

$$E_d=|I_k-I|\le \frac{2Q}{exp\left(\frac{2\pi a_1}{k}\right)-1}.$$

For the best approximation of $E_d$ in [47], the authors obtained the best parameters. They obtained the optimal values of $\eta$ and k by asymptotically equalizing the discretization and the truncation errors. The optimal values of the parameters given by

$$k=\frac{1}{M_p{(1+8\omega )}^{1/2}},\varrho =\frac{\pi M_p}{4T{(8\omega +1)}^{1/2}},\omega =\frac{T}{t_0},$$

thus, we have the following estimate

$${\mathsf{\Gamma }}_p\left(est\right)=|C_{App}(x,t)-C(x,t)|=O\left(e^{-\left(2\pi M_p\frac{1}{\sqrt{8\omega +1}}\right)}\right),$$

(29)

For more details, readers are referred to [47]. In our numerical experiments, $M_p$ is defined as $m=\frac{2M_p}{k}$, where k and m correspond to the quadrature step and the nodes, respectively.

2.4.2. Hyperbolic Contour

For the contour ${\mathsf{\Gamma }}_h$, the authors in [48] obtained the optimal parameters, with the error estimate given as

$${\mathsf{\Gamma }}_h\left(est\right)=|C_{App}(x,t)-C(x,t)|=O\left(\ell \left({\sigma }_r{M}_h\right)e^{-\varpi M_h}\right),$$

(30)

where $\varpi =\frac{\overline{r}}{b}(1-\vartheta )$,  $\overline{r}=2r\pi$,  $0<\vartheta <1$, $\tau =t_0/T$, $\ell \left(x\right)=max\left(log\right(1/x),1)$, ${\sigma }_r=(\tau sin(\theta -r))\frac{\vartheta \overline{r}}{\mathcal{U}}$, where $\mathcal{U}$ satisfies $cosh\left(\mathcal{U}\right)(\vartheta \tau sin(\theta \left)\right)=1$, $k=\mathcal{U}/M_h,$ $0<t_0\le t\le T$, and $\chi =\vartheta \overline{r}{M}_h/\left(\mathcal{U}T\right).$

2.4.3. (ii) Stehfest’s Method

In this section, we utilize Stehfest’s method for inverting the LT numerically. Stehfest’s method was developed in the 1960s. It is one of the most popular methods for numerical inversion of the LT. This method is fast and gives good results for smooth functions. The solution $C(x,t)$ can be approximated via Stehfest’s method as

$$C_{App}(x,t)=\frac{ln2}{t}\sum _{i=1}^{M_s}{w}_i\widehat{C}(x,s),$$

(31)

where $w_i$ are given by

$$w_i={(-1)}^{\frac{M_s}{2}+i}\sum _{h=\lfloor \frac{i+1}{2}\rfloor }^{min(i,\frac{M_s}{2})}\frac{h^{\frac{M_s}{2}}(2h!)}{(\frac{M_s}{2}-h)!h!(h-1)!(i-h)!(2h-i)!}.$$

(32)

Here, $M_s$, referred to as Stehfest’s number, should be even. Solving (20) and (21) for the parameters $s=\frac{ln2}{t}i,i=1,2,...,M_s$ in the Laplace domain. Using (31), the original solution can be obtained.

The Gaver–Stehfest algorithm has some attractive properties: (i) the approximations of any function using this algorithm are linear in values of the transform; (ii) it requires the values of the transform function for real s only; (iii) the computation of the coefficients is very easy; and (iv) the approximations using this algorithm are exact for constant functions. In the literature, use of this algorithm is reported by many authors [50,51], where it has been demonstrated that the method converges very rapidly for non-oscillatory functions.

2.4.4. Error Analysis

In this section, we present an error analysis of our proposed numerical scheme. In step (a), we use the LT which is analytic and in this step no error occurs. Errors in our method occur in step (b) and (c). In step (b), we use the local RBs method, which has an error estimate $O\left({\sigma }^{\frac{1}{\epsilon Q_d}}\right),0<\sigma <1,$ where $Q_d$ and $\epsilon$ denote the fill distance and the shape parameter [52]. In step (c), Stehfest’s numerical inversion for the LT is used. An error analysis of this method is described in [53,54], in which the authors report the effect of parameters on the numerical efficiency and accuracy. They concluded that “If $\beta$ significant digits are required, suppose $M_s=\lceil 2.2\beta \rceil$ positive integer. Set the precision of the system at $\upsilon =\lceil 1.1M\rceil$. Then given the transform $\widehat{C}(x,s)$ and the argument t, calculate $C_{App}(x,t)$ in (31)”. On the basis of these conclusions, the following error estimation is obtained:

Remark 2

([43]). If the error of the input data is ${10}^{-\upsilon +1}\le \frac{\parallel \widehat{C}-C_{App}\parallel }{\parallel C_{App}\parallel }\le {10}^{-\upsilon }$, with $M_S$ a positive even integer via $\upsilon =\lceil 1.1M_s\rceil$. Then the error is ${10}^{-\beta +1}\le \frac{\parallel \widehat{C}-C_{App}\parallel }{\parallel C_{App}\parallel }\le {10}^{-\beta }$, where $M_s=\lceil 1.1\beta \rceil .$

The key steps of the proposed numerical scheme are presented in the following Algorithm 1:

Algorithm 1: The necessary steps of the proposed numerical scheme

1: Input: The computational domain, the time derivative in [0, 1], the final time, the initial shape parameter and the other parameter of the given model, the inhomogeneous function, and initial boundary conditions. 2: Step 1: Reduce the problem defined in Equations (1) and (2) to an equivalent time-independent problem defined in Equations (7) and (8) by employing the Laplace transform. 3: Step 2: Discretize the linear differential operator $\mathcal{L}$ and boundary operator ${\mathcal{L}}_B$ using Equations (18) and (19). 4: Step 3: Solve the system of Equations (20) and (21) in parallel for each point s. 5: Step 4: Compute the approximate solution using (26) or (31). 6: Output: The approximate solution is $C_{App}(x,t)$.

3. Stability

In this section, we discuss the stability of the system (20) and (21) by writing it in discrete form as

$$\mathbf{G}\widehat{\mathbf{C}}=\mathbf{p},$$

(33)

where $\mathbf{G}$ is obtained using the local radial basis functions method. The stability constant for the discrete system (33) is defined as

$$\sigma =\underset{\widehat{\mathbf{C}}\ne 0}{sup}\frac{\parallel \widehat{\mathbf{C}}\parallel }{\parallel \mathbf{G}\widehat{\mathbf{C}}\parallel },$$

(34)

where the value of the constant $\sigma$ is finite for any discrete norm $\parallel .\parallel$ defined on $R^N$. Equation (34) implies

$${\parallel \mathbf{G}\parallel }^{-1}\le \frac{\parallel \widehat{\mathbf{C}}\parallel }{\parallel \mathbf{G}\widehat{\mathbf{C}}\parallel }\le \sigma .$$

(35)

For the pseudoinverse ${\mathbf{G}}^{†}$ of $\mathbf{G}$, we have

$$\parallel {\mathbf{G}}^{†}\parallel =\underset{\mathbf{v}\ne 0}{sup}\frac{\parallel {\mathbf{G}}^{†}\mathbf{v}\parallel }{\parallel \mathbf{v}\parallel }.$$

(36)

Therefore,

$$\parallel {\mathbf{G}}^{†}\parallel \ge \underset{\mathbf{v}=\mathbf{G}\widehat{\mathbf{C}}\ne 0}{sup}\frac{\parallel {\mathbf{G}}^{†}\mathbf{G}\widehat{\mathbf{C}}\parallel }{\parallel \mathbf{G}\widehat{\mathbf{C}}\parallel }=\underset{\widehat{\mathbf{C}}\ne 0}{sup}\frac{\parallel \widehat{\mathbf{C}}\parallel }{\parallel \mathbf{G}\widehat{\mathbf{C}}\parallel }=\sigma .$$

(37)

Equations (35) and (37) give the upper and lower bounds for $\sigma$. For the system (33), the pseudoinverse may be difficult to calculate, but it ensures the stability. We can utilize the command condest in MATLAB for estimating $\parallel {\mathbf{G}}^{-1}{\parallel }_{\infty }$, thus

$$\sigma =\frac{condest\left({\mathbf{G}}^{\prime }\right)}{{\parallel \mathbf{G}\parallel }_{\infty }}$$

(38)

4. Numerical Examples

To validate the proposed numerical methods, we consider three test problems. The MQ RBFs defined by $\phi =\sqrt{r^2+{\vartheta }^2}$ are used in this paper. For the optimal shape parameter, the uncertainty principal due to [55] is utilized. We performed our experiments in MATLAB R2018a on a Windows 10 (64 bit) PC equipped with an Intel(R) Core(TM) i5-3317U CPU @ 1.70 GHz and with 4 GB of RAM. The performance of our method is tested using the $L_{\infty }$ error norm which is denoted and defined by

$$Erro{r}_{max}={\parallel C(x_i,t)-C_{App}(x_i,t)\parallel }_{\infty },$$

where $C_{App}(x,t)$ and $C(x,t)$ are the numerical and analytical solutions, respectively.

4.1. Problem 1

Here, we consider a $1\mathrm{D}$ diffusion equation of the form

$$\frac{\partial C(x,t)}{\partial t}-\frac{{\partial }^{2}C(x,t)}{\partial x^2}=exp\left(x\right)(2t-t^2),t>0,0<x<1,$$

(39)

with the initial and boundary data

$$C(x,0)=0,$$

and

$$C(0,t)=t^2,$$

$$C(1,t)=exp\left(1\right)t^2.$$

The analytical solution is

$$C(x,t)=exp\left(x\right)t^2.$$

The results obtained along the parabolic contour are shown in

Table 1. The $Erro{r}_{max}$ using the proposed scheme along the path ${\mathsf{\Gamma }}_p,$ at $t=1$.

${\mathit{M}}_{\mathit{p}}$	N	n	${\mathit{Error}}_{\mathit{max}}$	$\mathit{\vartheta }$	C. Time (s)
30	70	10	1.18 $\times {10}^{-2}$	6.4	0.353447
40			2.55 $\times {10}^{-4}$	6.4	0.416162
50			2.46 $\times {10}^{-4}$	6.4	0.538895
60	20	8	3.60 $\times {10}^{-5}$	1.4	0.271687
	30		5.01 $\times {10}^{-5}$	2.1	0.324185
	40		6.04 $\times {10}^{-5}$	2.9	0.386367

, and along the hyperbolic contour are shown in

Table 2. The $Erro{r}_{max}$ using the proposed scheme along the path ${\mathsf{\Gamma }}_h$ at $t=1$.

${\mathit{M}}_{\mathit{h}}$	N	n	${\mathit{Error}}_{\mathit{max}}$	$\mathit{\vartheta }$	C. Time (s)
30	20	10	5.49 $\times {10}^{-4}$	1.7	0.202415
40			1.67 $\times {10}^{-4}$	1.7	0.225109
50			1.45 $\times {10}^{-4}$	1.7	0.308366
60			1.53 $\times {10}^{-4}$	1.7	0.436627
70	30	12	7.50 $\times {10}^{-5}$	3.0	0.843163
	50		7.42 $\times {10}^{-5}$	5.1	2.020373
	80		9.54 $\times {10}^{-5}$	8.3	4.864506

. Similarly, the results obtained using Stehfest’s method are presented in

Table 3. The $Erro{r}_{max}$ using the proposed scheme at $t=1$.

${\mathit{M}}_{\mathit{s}}$	N	n	${\mathit{Error}}_{\mathit{max}}$	$\mathit{\vartheta }$	CPU (s)
8	70	10	6.22 $\times {10}^{-2}$	6.4	0.346949
10			9.56 $\times {10}^{-4}$	6.4	0.334023
12			7.99 $\times {10}^{-4}$	6.4	0.281686
14			2.16 $\times {10}^{-4}$	6.4	0.288764
16	50	12	1.61 $\times {10}^{-4}$	5.1	0.225202
	60		1.81 $\times {10}^{-4}$	6.2	0.244202
	70		5.04 $\times {10}^{-4}$	7.3	0.336965
	80		3.47 $\times {10}^{-4}$	8.3	0.346112
	90		6.01 $\times {10}^{-4}$	9.4	0.418864

. For the contours, we obtained the quadrature nodes using the MATLAB command $\upsilon =-M_h:k:M_h$. In our experiments, different parameter values are used which are $\varsigma =0.15410,\chi =2,c=0.1,r=0.13870,\vartheta =0.10$, and $t\in [0.5,5].$ Various nodes in the local and global domains are selected for the numerical experiments. Figure 1a–d depicts the numerical and exact solutions of problem 1, computed using the hyperbolic contour with parameter values $M_h=50,N=30,n=12.$ Figure 2a,b displays the exact and numerical space-time solutions of problem 1, computed along the hyperbolic contour with parameter values $M_h=50,N=40,n=12.$ Figure 3a shows a comparison between the errors of problem 1 obtained using the three inversion techniques. It is observed that the performance of the two contour integration methods is better than for Stehfest’s method. Figure 3b shows the plot of error versus the nodes $M_p$ along the parabolic contour. Similarly Figure 4a shows the plot of error versus $M_h$ using the hyperbolic contour. A comparison between the observed and theoretical error is shown in Figure 4b; a good agreement between the error and the error estimate is observed. Figure 4c shows a plot of error versus $M_s$ using Stehfest’s method. In Figure 5a, the plot shows error versus t using the hyperbolic path. Figure 5b depicts a plot of error versus t using Stehfest’s method. From all the figures and tables, a very high accuracy and stability is evident. However, it is observed that the results along the hyperbolic path are more accurate than that of the results along the parabolic path and using Stehfest’s method. The results demonstrate that the method is capable of approximating the solution of this type of equation efficiently.

4.2. Problem 2

Here, a $1\mathrm{D}$ equation of the form [10] is considered

$$\frac{\partial C(x,t)}{\partial t}+\overline{\theta }\frac{\partial C(x,t)}{\partial x}-\overline{\alpha }\frac{{\partial }^{2}C(x,t)}{\partial x^2}=0,\mathrm{where}0<x<1,t>0,$$

(40)

with the initial and boundary data

$$C(x,0)=\overline{a}exp(-\overline{c}x),$$

and

$$C(0,t)=\overline{a}exp\left(\left(\overline{b}t\right)\right),$$

$$C(1,t)=\overline{a}exp\left((\overline{b}t-\overline{c})\right),$$

the analytical solution of the problem is

$$C(x,t)=\overline{a}exp\left((\overline{b}t-\overline{c}\xi )\right),$$

where

$$\overline{c}=\frac{-\overline{\beta }+\sqrt{{\overline{\beta }}^2+4\overline{\alpha }\overline{b}}}{2\overline{\alpha }}.$$

The ratio of the two constants $\overline{\theta }$ and $\overline{\alpha }$ is denoted by $P_e$, called the P$\acute{e}$clet number, i.e., $P_e=\left|\frac{\overline{\theta }}{\overline{\alpha }}\right|$. For a large value of $P_e$, the convection term is dominant and, for a small value of $P_e$, the diffusion term is dominant [10]. The results using the proposed numerical scheme along the parabolic contour are shown in

Table 4. The $Erro{r}_{max}$ using the proposed scheme along the path ${\mathsf{\Gamma }}_p$ at $t=1$.

	${\mathit{M}}_{\mathit{p}}$	N	n	${\mathit{Error}}_{\mathit{max}}$	$\mathit{\vartheta }$	C. Time (s)
	40	90	10	1.12 $\times {10}^{-4}$	8.3	0.450731
	50			3.54 $\times {10}^{-5}$	8.3	0.671841
	60			3.53 $\times {10}^{-5}$	8.3	1.206206
	70			3.53 $\times {10}^{-5}$	8.3	1.150397
	80	60	10	7.61 $\times {10}^{-5}$	5.5	0.839738
		70		3.61 $\times {10}^{-5}$	6.4	0.997628
		100		1.07 $\times {10}^{-5}$	9.2	1.612720
[56]				4.41 $\times {10}^{-3}$

, and the results obtained using Stehfest’s method are presented in

Table 5. The $Erro{r}_{max}$ using the proposed scheme at $t=1$.

	${\mathit{M}}_{\mathit{s}}$	N	n	${\mathit{Error}}_{\mathit{max}}$	$\mathit{\vartheta }$	CPU (s)
	8	70	10	1.10 $\times {10}^{-3}$	6.4	0.317848
	10			1.08 $\times {10}^{-4}$	6.4	0.332286
	12			3.25 $\times {10}^{-5}$	6.4	0.323966
	14			3.64 $\times {10}^{-5}$	6.4	0.331252
	16	50	12	1.53 $\times {10}^{-4}$	5.1	0.268543
		60		1.34 $\times {10}^{-4}$	6.2	0.321235
		70		8.73 $\times {10}^{-5}$	7.3	0.365170
		80		7.73 $\times {10}^{-5}$	8.3	0.443169
		90		4.96 $\times {10}^{-5}$	9.4	0.533850
		100		2.72 $\times {10}^{-5}$	9.4	0.558243
[56]				4.41 $\times {10}^{-3}$

. The same set of parameter values for the contours is used in this problem. Figure 6a,b depicts the numerical solutions of problem 2 computed using the parabolic contour at different time levels, with $P_e=10$ and $P_e=20$ and other parameter values $M_p=70,N=30,n=10.$ Figure 7a shows a comparison between the errors of problem 2 obtained using the three inversion techniques with various values of $P_e$. It is observed that the performance of the two contour integration methods is better than for Stehfest’s method. Figure 7b shows the plot of error versus the nodes $M_p$ along the parabolic contour. Figure 8a shows the plot of error versus $M_h$ using the hyperbolic contour. A comparison between the observed and theoretical error is shown in Figure 8b; a good agreement between error and error estimate is observed. Figure 8c shows the plot of error versus $M_s$ using Stehfest’s method. In Figure 9a, the plot shows error versus t using the hyperbolic path. Figure 9b shows a plot of error versus t using Stehfest’s method. From all the figures and tables, a very high accuracy and stability is evident. However, it is observed that the results along the hyperbolic path are more accurate than those for the results along the parabolic path and using Stehfest’s method. The results demonstrate that the method is capable of approximating the solution of this type of equation efficiently.

4.3. Problem 3

Here, a $1\mathrm{D}$ equation of the form [10] is considered

$$\frac{\partial C(x,t)}{\partial t}+\overline{\theta }\frac{\partial C(x,t)}{\partial x}-\overline{\alpha }\frac{{\partial }^{2}C(x,t)}{\partial x^2}=0,0<x<2,t>0,$$

(41)

with the initial and boundary data

$$C(x,0)=sin\left(x\right),$$

and

$$C(0,t)=exp(-\overline{\alpha }t)sin(-\overline{\theta }t),$$

$$C(2,t)=exp(-\overline{\alpha }t)sin(2-\overline{\theta }t),$$

the analytical solution is given as

$$C(x,t)=exp(-\overline{\alpha }t)sin(x-\overline{\theta }t).$$

The results obtained for various values of $P_e$ along the hyperbolic contour are shown in

Table 6. $Erro{r}_{max}$ of the proposed method for problem 3 using ${\mathsf{\Gamma }}_h$ at $t=1$ for various P$\acute{e}$clet numbers.

	${\mathit{M}}_{\mathit{h}}$	N	n	${\mathit{P}}_{\mathit{e}}=100$	${\mathit{P}}_{\mathit{e}}=1000$	${\mathit{P}}_{\mathit{e}}=\mathrm{10,000}$	${\mathit{P}}_{\mathit{e}}=\mathrm{20,000}$
	30	40	8	9.83 $\times {10}^{-2}$	9.83 $\times {10}^{-2}$	9.83 $\times {10}^{-2}$	9.83 $\times {10}^{-2}$
	40			2.20 $\times {10}^{-3}$	2.20 $\times {10}^{-3}$	2.20 $\times {10}^{-3}$	2.20 $\times {10}^{-3}$
	50			5.70 $\times {10}^{-3}$	5.70 $\times {10}^{-3}$	5.70 $\times {10}^{-3}$	5.70 $\times {10}^{-3}$
	60			4.37 $\times {10}^{-4}$	4.37 $\times {10}^{-4}$	4.37 $\times {10}^{-4}$	4.38 $\times {10}^{-4}$
	70			2.96 $\times {10}^{-4}$	2.96 $\times {10}^{-4}$	2.96 $\times {10}^{-4}$	2.96 $\times {10}^{-4}$
	90			1.79 $\times {10}^{-5}$	1.76 $\times {10}^{-5}$	1.79 $\times {10}^{-5}$	1.80 $\times {10}^{-5}$
	100			1.25 $\times {10}^{-5}$	1.24 $\times {10}^{-5}$	1.29 $\times {10}^{-5}$	1.29 $\times {10}^{-5}$
	90	20	10	1.67 $\times {10}^{-4}$	1.95 $\times {10}^{-4}$	2.02 $\times {10}^{-4}$	2.02 $\times {10}^{-4}$
		30		2.59 $\times {10}^{-5}$	4.09 $\times {10}^{-5}$	5.12 $\times {10}^{-5}$	5.20 $\times {10}^{-5}$
		40		4.38 $\times {10}^{-5}$	4.62 $\times {10}^{-5}$	5.16 $\times {10}^{-5}$	5.20 $\times {10}^{-5}$
		50		5.93 $\times {10}^{-5}$	6.40 $\times {10}^{-5}$	6.41 $\times {10}^{-5}$	6.40 $\times {10}^{-5}$
[10]					2.87 $\times {10}^{-7}$	2.68 $\times {10}^{-4}$	7.40 $\times {10}^{-4}$

. The results obtained for $P_e$ using Stehfest’s numerical inversion method are presented in

Table 7. $Erro{r}_{max}$ using the proposed scheme at $t=1$ and $P_e=100$.

${\mathit{M}}_{\mathit{s}}$	N	n	${\mathit{Error}}_{\mathit{max}}$	$\mathit{\vartheta }$	CPU (s)
8	70	10	1.60 $\times {10}^{-1}$	6.4	0.334523
10			6.29 $\times {10}^{-2}$	6.4	0.426634
12			8.80 $\times {10}^{-3}$	6.4	0.347591
14			2.40 $\times {10}^{-3}$	6.4	0.479109
16	50	12	2.20 $\times {10}^{-4}$	5.1	0.361524
	60		2.61 $\times {10}^{-4}$	6.2	0.322293
	70		4.20 $\times {10}^{-4}$	7.3	0.376275
	80		3.90 $\times {10}^{-4}$	8.3	0.447432
	90		3.44 $\times {10}^{-4}$	9.4	0.718522
	100		4.98 $\times {10}^{-4}$	9.4	0.672432

. Here, we use the same set of parameters used for problem 1 and problem 2. Figure 10a–f depict the numerical and exact solutions of problem 3, computed using Stehfest’s method, with $M_s=16,N=100,n=12,$ and $P_e=2$, at different time levels. Figure 11a,b shows the exact and numerical space-time solutions of problem 3 using Stehfest’s method at $t=1$ with $M_s=16,N=30,n=12,$ and $P_e=2.$ Figure 12a shows a comparison between the errors of problem 3 obtained using the three inversion techniques with various values of $P_e$. It is observed that the performance of the two contour integration methods is better than for Stehfest’s method. Figure 12b shows a plot of error versus the nodes $M_p$ along the parabolic contour. Figure 13a shows a plot of error versus $M_h$ using the hyperbolic contour. A comparison between the observed and theoretical error is shown in Figure 13b; a good agreement between the error and the error estimate is observed. Figure 13c shows the plot of error versus $M_s$ using Stehfest’s method. In Figure 14a, the plot shows error versus t using the hyperbolic path. Figure 14b depicts a plot of error versus t using Stehfest’s method. From all the figures and tables, a very high accuracy and stability is evident. However, it is observed that the results along the hyperbolic path are more accurate than the results along the parabolic path and using Stehfest’s method.

5. Conclusions

In this paper, a numerical method based on the Laplace transform and local radial basis functions was developed for numerical simulation diffusion equations. The method employed the Laplace transform and inverse Laplace transform methods to cope with the time derivative. This method first uses the Laplace transform to convert the given problem into an inhomogeneous problem in the Laplace domain. Then it utilizes the local radial basis functions for discretization of the spatial derivatives. Finally, it uses the numerical inverse Laplace transform methods to convert the solution obtained in the Laplace domain into a time domain. The present method is almost exact in time, having no stability restrictions, which is a common issue in finite-difference time-stepping methods. In these methods, for optimal accuracy, we need a very small time step. The local radial basis functions makes this method very interesting and easy to implement. An error analysis and the results of computation demonstrate that the present method is an efficient and competitive method for the solution of diffusion problems, with the benefits that it is a simple, mesh-less, easy to implement and highly accurate numerical method.