On the Combination of the Laplace Transform and the Integral Equation Method to the 3D Parabolic Initial Boundary Value Problem

1. Introduction

The boundary integral equation (BIE) method is a very powerful approach for the numerical solution of various boundary value problems (BVPs). The main advantage of BIE method consists in the dimensionality decrease of the given differential problem: the BVP is reduced to the BIE, where the unknown function is defined only on the domain boundary [1]. Clearly, the considered differential equation needs to have the fundamental solution and to be homogeneous. For the numerical solution of such BIE, effective numerical methods are developed, for example, projection methods [1].

In the case of non-stationary BVPs, there are additional difficulties caused by the presence of time as an independent variable. There are several ways to apply BIE to such BVPs. One approach involves a fundamental solution of the time-dependent differential equation. Then, by a direct or indirect BIE method, the initial BVP can be reduced to a time-boundary integral equation. The numerical solution of such a BIE is more difficult than in the stationary case. The most popular method for time-boundary integral equations is the convolution quadrature method suggested by Christian Lubich in the 1980s [2].

Another so-called two-steps methods consist of the semi-discretization of the given initial BVP with respect to the time variable. As a result, the set of stationary BVPs for elliptic equations is obtained. This time discretization can be achieved using approaches such as finite-difference approximations (e.g., the Rothe method [3]) or integral transforms (e.g., the Laguerre transform [4,5,6,7], the Laplace transform [8]). In the second step, which addresses the spatial variable, various techniques are available, including the BIE method. Two-step methods offer several advantages, such as dimension reduction and the avoidance of volume integrals. The finite-difference semi-discretization is the simplest approach, which gives the numerical solution in a fixed set of time moments. In the case of the integral transforms, we have an approximation for arbitrary time, but it is necessary to calculate the inverse transform numerically.

The Laguerre transform for a parabolic initial BVP leads to stationary BVPs for a recurrent sequence of elliptic equations, and to apply the BIE method, one needs to find the fundamental sequence [4,7]. The inverse Laguerre transform has the form of the Fourier-Laguerre series, and its summation is a complicated ill-posed problem, especially for long-time moments. In the case of the Laplace transform to a parabolic initial BVP, stationary BVPs for the Helmholtz-type equations with complex parameters can be obtained. The inverse Laplace transform is defined as the Bromwich integral on the complex plane, and there are several numerical methods for its calculation [9].

In this paper, we use the two-step approach based on the Laplace transform and the BIE method to solve parabolic BVPs in 3D domains. To calculate the inverse transform, the sinc-quadrature rule suggested in [10] is applied. Our contribution is to reduce computational costs by selecting optimal values of the quadrature formula parameters and applying an efficient method for numerically solving the resulting BIEs.

The outline of the present work is as follows. In Section 2, we apply the Laplace transform to the parabolic initial boundary value problem and describe the sinc-quadrature for the numerical inverse transform. Two ideas for decreasing computational cost are presented in subsections 2.1 and 2.2. In Subsection 2.1 it is shown that due to a certain symmetry of the sinc-quadrature nodes, the number of stationary problems can be reduced almost twice. In Subsection 2.2, it reflects how the choice of integration contour in the complex plane influences the precision of the sinc-quadrature. In Section 3, we apply the indirect BIE method to stationary elliptic problems. The unknown solution is presented in the form of the double-layer potential, and the BIE of the second kind is obtained. Taking into account that the boundary surface is diffeomorphic to the unit sphere we apply the Nyström method based on the Wienert’s quadrature rules. Section 4 presents numerical examples to clarify our approach and its optimization.

Before closing this section, we formulate the problem to be studied. Let $D\subset {\mathbb{R}}^3$ be a simply connected region with a smooth boundary $\Gamma$. It is necessary to find function $u\in C^{2,1}(D\times (0,\infty ))\cap C(\overline{D}\times [0,\infty ))$, which satisfies the heat equation

$${\displaystyle \frac{\partial u}{\partial t}}(x,t)=\Delta u(x,t),(x,t)\in D\times (0,\infty ),$$

(1)

the initial condition

$$u(x,0)=0,x\in D$$

(2)

and the Dirichlet boundary value condition

$$u(x,t)=g(x,t),(x,t)\in \Gamma \times (0,\infty ).$$

(3)

Assume that the given function g is bounded, continuous, and satisfies the compatibility condition $g(x,0)=0,x\in \Gamma$.

We consider surfaces $\Gamma$, diffeomorphic to the unit sphere

$${\mathbb{S}}^2=\left\{p(\theta ,\varphi )=(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta ),(\theta ,\varphi )\in [0,\pi ]\times [0,2\pi )\right\}$$

described by an analytic function $q:{\mathbb{S}}^2\to \Gamma$ with a non zero Jacobian J.

2. Time Semi-Discretization via Laplace Transform

The Laplace transform of the function $f\left(t\right)$ is given by

$${\mathcal{L}}_t\left(f\right)=F\left(s\right)={\int }_0^{\infty }{e}^{-st}f\left(t\right)dt,s=\sigma +i\tau ,s\in \mathbb{C}.$$

(4)

The integral (4) is convergent for $\mathrm{Re}\left(s\right)>a_0$, where $a_0$ is the order of growth of the function $f\left(t\right)$, and $F\left(s\right)$ is an analytic function.

For the known image F the original f can be reconstructed by using the inverse Laplace transform, described by the Bromwich integral

$${\mathcal{L}}_s^{-1}\left(F\right)=f\left(t\right)={\displaystyle \frac{1}{2\pi i}}{\int }_C{e}^{st}F\left(s\right)ds,$$

(5)

where C is the suitable integration contour (see [8,10,11]).

A popular strategy to use the Laplace transform for the heat problems is next:

1.

Apply the Laplace transform in time to the initial boundary value problem to obtain boundary value problems for the Helmholtz-type equations.

2.

Build an effective solver for stationary problems.

3.

Reconstruct time-domain solution via numerical inversion of the Laplace transform.

One approach to approximating the inverse Laplace transform was proposed in [10]. If F can be analytically continued to the set $\mathbb{C}\setminus {\Sigma }_{\delta }$, where

$${\Sigma }_{\delta }=\{s\in \mathbb{C}:|arg(-s)|\le \delta ,0<\delta <{\displaystyle \frac{\pi }{2}}\}$$

and there exists $M>0$ such that

$$\left|F\left(s\right)\right|\le {\displaystyle \frac{M}{\left|s\right|}},s\in \mathbb{C}\setminus {\Sigma }_{\delta },$$

then to approximate the inverse Laplace transform of the function F, a quadrature formula is proposed based on the use of sinc-quadrature for integral (5) with a special integration contour (see Figure 1)

$$\gamma \left(\omega \right)=\lambda \left(1-sin(\alpha +i\omega )\right),\omega \in \mathbb{R}.$$

(6)

Here $\lambda >0$, $0<\alpha <{\displaystyle \frac{\pi }{2}}-\delta$ are arbitrary parameters that define the geometry of the contour (6).

Using contour (6) to parametrize integral (5), we obtain

$${\mathcal{L}}_s^{-1}\left(F\right)=f\left(t\right)=-{\displaystyle \frac{1}{2\pi i}}{\int }_{-\infty }^{\infty }{e}^{\gamma \left(\omega \right)t}F\left(\gamma \left(\omega \right)\right){\gamma }^{\prime }\left(\omega \right)d\omega .$$

(7)

Let $N_t>0,N_t\in \mathbb{N}$, $h_N=ln{N}_t/N_t$, ${\omega }_j=h_{N}j$,$j=-N_t,...,N_t$. Integral (7) can be approximated using the following quadrature formula [10]

$$f\left(t\right)\approx \left(T_{N_t}F\right)\left(t\right)={\displaystyle \frac{h_N}{2\pi i}}\sum _{j=-N_t}^{N_t}{e}^{\gamma \left({\omega }_j\right)t}F\left(\gamma \left({\omega }_j\right)\right){\gamma }^{\prime }\left({\omega }_j\right).$$

(8)

Let us denote $s_j=\gamma \left({\omega }_j\right)$ and ${\gamma }_j=h_N{\gamma }^{\prime }\left({\omega }_j\right)/2\pi i$, then we obtain

$$f\left(t\right)\approx \left(T_{N_t}F\right)\left(t\right)=\sum _{j=-N_t}^{N_t}{\gamma }_j{e}^{t{s}_j}F\left(s_j\right).$$

(9)

Note, when computing $\left(T_{N_t}F\right)\left(t\right)$ for different values of t, one can use the same set of values $F\left(s_j\right)$. The approximation error of (9) is shown to behave like $\mathcal{O}\left(e^{-c{N}_t/ln{N}_t}\right)$ and possess stability to the perturbations of $F\left(s_j\right)$. This is especially important when values $F\left(s_j\right)$ are computed numerically [10].

Since the solution of the non-stationary problem (1)–(3) u is bounded with respect to the time variable, i.e., its order of growth is equal to 0, the Laplace transform with respect to time can be applied to both parts of equation (1). Taking into account property ${\mathcal{L}}_t\left(f^{\prime }\left(t\right)\right)=sF\left(s\right)-f\left(0\right)$ and zero initial condition, we obtain the following equation for the Laplace image $U_s\left(x\right)={\mathcal{L}}_t\left(u(x,t)\right)$

$$\Delta U_s\left(x\right)-s{U}_s\left(x\right)=0,x\in D.$$

(10)

On the boundary of the domain, the function $U_s$ satisfies the following condition

$$U_s\left(x\right)=G_s\left(x\right),x\in \Gamma ,$$

(11)

where $G_s\left(x\right)={\mathcal{L}}_t\left(g(x,t)\right)$. Thus, for $U_s$ we get a boundary value problem (10)–(11) for the Helmholtz equation with a complex wavenumber ${\kappa }^2=-s$.

Applying the described approach for the inverse Laplace transform, in order to find an approximate solution of problem (1)–(3), it is necessary to compute

$$U_j\left(x\right)=U_{s_j}\left(x\right),j=-N_t,...,N_t,$$

that is, to solve a set of $2N_t+1$ problems (10)–(11) for $s=s_j$

$$\Delta U_j\left(x\right)-s_j{U}_j\left(x\right)=0,x\in D,$$

(12)

$$U_j\left(x\right)=G_j\left(x\right),x\in \Gamma .$$

(13)

Here $G_j\left(x\right)=G_{s_j}\left(x\right)$. It is important to emphasize that problems (12)–(13) are independent of each other, enabling their parallel solution.

In [8,11] it was shown that the image of the solution to the heat problem $U_s$, as a function of the complex argument s, can be analytically continued to the set $Z=\mathbb{C}\setminus (-\infty ,0]$, and there exists a constant $M>0$ such that

$$|U_s(·)|\le {\displaystyle \frac{M}{\left|s\right|}},s\in Z.$$

Thus, in our case, we can apply the approach from [10] and use the contour (6) for any $\lambda >0$ and $0<\alpha <{\displaystyle \frac{\pi }{2}}$.

Note, in order to solve problems (12)–(13) it is necessary to have boundary functions $G_j$, i.e., have the Laplace image $G_s$ of the original boundary condition g. If $G_s$ is not available in a closed form, it can be approximated using various techniques (see [12,13,14]). We leave the approximation of $G_s$ beyond the scope of the current article and will use examples of g with a known Laplace transform for the numerical experiments.

Recalling problems (12)–(13) are 3D stationary boundary value problems, it is easy to see that solving them numerically may pose a significant computational effort. The main motivation for this article was to suggest certain ideas for decreasing the amount of computational work, as described further.

2.1. Reducing Number of Stationary Problems

It is easy to notice that ${\omega }_{-j}=-{\omega }_j,j=1,...,N_t$.

Then

$$s_{-j}=\gamma \left({\omega }_j\right)=\lambda \left(1-sin(\alpha +i{\omega }_j)\right)=$$

$$=\lambda \left(1-sin(\alpha -i{\omega }_{N_t-j})\right)=\overline{\lambda (1-sin(\alpha +i{\omega }_{N_t-j})})=\overline{s_j},j=1,...,N_t.$$

(14)

We use the fact that $\overline{sin\left(z\right)}=sin\left(\overline{z}\right)$ for any complex z.

Thus, quadrature nodes (14) are pairwise conjugate, except for the node $s_0$. This allows reducing the solution of the set of $2N_t+1$ stationary problems to $N_t+1$ problems.

We show that

$$U_{-j}\left(x\right)=\overline{U_j\left(x\right)},j=1,...,N_t.$$

Theorem 1. 

Let U be a solution of the problem

$$\Delta U\left(x\right)-sU\left(x\right)=0,x\in D,$$

(a)

$$U\left(x\right)=G\left(x\right),x\in \Gamma .$$

(b)

Then $\overline{U}$ is a solution of the problem

$$\Delta \overline{U\left(x\right)}-\overline{s}\overline{U\left(x\right)}=0,x\in D,$$

(c)

$$\overline{U\left(x\right)}=\overline{G\left(x\right)},x\in \Gamma .$$

(d)

Proof. 

Statement (d) follows directly from (b). Let us show that (c) follows from (a).

We denote

$$s=a+ib,U\left(x\right)=V\left(x\right)+iW\left(x\right),$$

$$\overline{s}=a-ib,\overline{U\left(x\right)}=V\left(x\right)-iW\left(x\right).$$

Then (a) can be written as

$$\Delta V\left(x\right)+i\Delta W\left(x\right)=aV\left(x\right)-bW\left(x\right)+i\left(bV\right(x)+aW(x\left)\right)$$

Thus

$$\Delta V\left(x\right)=aV\left(x\right)-bW\left(x\right),\Delta W\left(x\right)=bV\left(x\right)+aW\left(x\right).$$

Then

$$\Delta \overline{U\left(x\right)}=\Delta V\left(x\right)-i\Delta W\left(x\right)=$$

$$=aV\left(x\right)-bW\left(x\right)-ibV\left(x\right)-iaW\left(x\right)=$$

$$=(a-ib)V\left(x\right)-i(a-ib)W\left(x\right)=$$

$$=(a-ib)(V\left(x\right)-iW\left(x\right))=\overline{s}\overline{U\left(x\right)},$$

which proves statement (c). □

Corollary 1. 

Solutions of the problems(12)-(13)with indices $-j$ and j are complex conjugates

$$U_{-j}\left(x\right)=\overline{U_j\left(x\right)},j=1,...,N_t.$$

Proof. 

Using a well-known fact $F\left(\overline{s}\right)=\overline{F\left(s\right)}$ (for real valued $f\left(t\right)$) it is easy to see that the boundary conditions of the problems with indices $-j$ and j are complex conjugates. Since $s_{-j}=\overline{s_j}$, it follows from the Theorem 1 the solutions of the problems (12) - (13) with indices $-j$ and j are also complex conjugates. □

Thus, it is sufficient to solve the stationary problems for indices $j=0,\dots ,N_t$, and the solutions for the indices $j=-N_t,\dots ,-1$ can be obtained automatically from the Corollary 1.

2.2. Integration Contour Parameters Optimization

As mentioned earlier, the integration contour (6) depends on the parameters $\lambda >0$ and $0<\alpha <{\displaystyle \frac{\pi }{2}}$. The Figure 2 and Figure 3 show the influence of the parameters $\alpha$ and $\lambda$ on the shape of the contour and placement of the nodes for $N_t=4$

Since the approximate solution of the 3D stationary problems requires a large amount of computations, it makes sense to select the parameters $\alpha$, $\lambda$ in such a way as to reduce the expected error.

To find the parameters $\alpha$, $\lambda$ for which the error is minimized, we will define search intervals for the optimal values of $\alpha$, $\lambda$ and construct a uniform grid of test values for $\alpha ,\lambda$

$$\alpha \in [{\alpha }_0,{\alpha }_1]\subset (0;\pi /2),\lambda \in [{\lambda }_0,{\lambda }_1]\subset (0;\infty ),$$

$${\alpha }_{\nu }={\alpha }_0+\nu h_{\alpha },\nu =0,...,N_{\alpha },N_{\alpha }\in \mathbb{N},h_{\alpha }=({\alpha }_1-{\alpha }_0)/N_{\alpha },$$

(15)

$${\lambda }_{\mu }={\lambda }_0+\mu h_{\lambda },\mu =0,...,N_{\lambda },N_{\lambda }\in \mathbb{N},h_{\lambda }=({\lambda }_1-{\lambda }_0)/N_{\lambda }.$$

(16)

We fix certain values of $x,t,N_t$ and select a Laplace transform pair of test functions $u^T(x,t)$ and $U^T(x,s)$. It is natural to select $u^T$ to be similar to the behavior of the boundary condition g. Then, for each pair of values $({\alpha }_{\nu },{\lambda }_{\mu })$, we compute the absolute or relative errors $E_{\mathrm{abs}},E_{\mathrm{rel}}$ of the numerical Laplace transform inversion (9) for $U^T(x,s)$ and find the values $({\alpha }_{min},{\lambda }_{min})$ for which

$$E_{abs,min}=min\left\{E_{abs}\right|\left({\alpha }_{\nu ,{\lambda }_{\mu }}\right),\nu =0,...,N_{\alpha },\mu =0,...,N_{\lambda }\},$$

(17)

$$E_{rel,min}=min\left\{E_{rel}\right|\left({\alpha }_{\nu ,{\lambda }_{\mu }}\right),\nu =0,...,N_{\alpha },\mu =0,...,N_{\lambda }\}.$$

(18)

Obtained contour parameters $({\alpha }_{min},{\lambda }_{min})$ are then used to define quadrature nodes $s_j$ and solve $N_t+1$ stationary problems. We do not provide an explicit recipe to define $[{\alpha }_0,{\alpha }_1]$ and $[{\lambda }_0,{\lambda }_1]$. For $[{\alpha }_0,{\alpha }_1]$ it seems natural to define ${\alpha }_0$ close to 0 and ${\alpha }_1$ close to $\pi /2$ and thus "scan" most of the $(0;\pi /2)$ interval. For $[{\lambda }_0,{\lambda }_1]$ it is empirically observed that increasing ${\lambda }_1$ stops finding different ${\lambda }_{min}$ after certain values of ${\lambda }_1$.

3. Stationary Boundary Value Problems Solver

In this section we consider the numerical solution of the stationary problems (12)–(13). We will apply the BIE method with later application of the Nyström method based on the quadrature rules for surface integrals, proposed by Wienert [15].

For brevity, we rewrite problems (12)–(13) as

$$\Delta U_s\left(x\right)-s{U}_s\left(x\right)=0,x\in D,s\in \{s_j,j=0,...,N_t\},$$

(19)

$$U_s\left(x\right)=G_s\left(x\right),x\in \Gamma .$$

(20)

The fundamental solution of the equation (19) has the form

$${\Phi }_s(x,y)={\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{e^{-k|x-y|}}{|x-y|}},x,y\in {\mathbb{R}}^3,x\ne y,$$

(21)

$$k=\sqrt{s},Re\left(\kappa \right)>0.$$

Since $s\notin (-\infty ,0]$, it is known that under suitable assumptions on the boundary $\Gamma$ and for sufficiently smooth boundary data $G_s$, the solution of the Dirichlet problem exists and is unique, see [16] and references therein. The solution of (19) can be written in the form of a double layer potential

$$U_s\left(x\right)=\underset{\Gamma }{\int }\phi \left(y\right){\Phi }_{\nu ,s}(x,y)ds\left(y\right),x\in \Gamma ,$$

(22)

where ${\Phi }_{\nu ,s}(x,y)={\displaystyle \frac{\partial {\Phi }_s(x,y)}{\partial \nu \left(y\right)}}$, $\phi \left(y\right)$ is the potential density and $\nu \left(y\right)$ is the unit outward normal vector to $\Gamma$ at the point y.

Potential (22) is a solution of the problem (12)–(13) if the density $\phi$ is a solution of the Fredholm integral equation of the second kind

$$-{\displaystyle \frac{1}{2}}\phi \left(x\right)+\underset{\Gamma }{\int }\phi \left(y\right){\Phi }_{\nu ,s}(x,y)ds\left(y\right)=G_s\left(x\right),x\in \Gamma .$$

(23)

For any $G_s\in C(\Gamma )$ equation (23) has a unique solution $\phi$ in $C(\Gamma )$ [16].

Since $\Gamma =\left\{x=q\left(\widehat{x}\right),\widehat{x}\in {\mathbb{S}}^2)\right\}$, we can obtain the parametrized integral equation on ${\mathbb{S}}^2$

$$-{\displaystyle \frac{1}{2}}\psi \left(\widehat{x}\right)+\underset{{\mathbb{S}}^2}{\int }\psi \left(\widehat{y}\right)K(\widehat{x},\widehat{y})ds\left(\widehat{y}\right)={\widehat{G}}_s\left(\widehat{x}\right),\widehat{x}\in {\mathbb{S}}^2,$$

(24)

where we denoted $\psi \left(\widehat{x}\right)=\phi \left(q\left(\widehat{x}\right)\right),\widehat{G_s}\left(\widehat{x}\right)=G_s\left(q\left(\widehat{x}\right)\right),\widehat{x}\in {\mathbb{S}}^2$ and

$${\displaystyle K(\widehat{x},\widehat{y})={\Phi }_{\nu ,s}(q\left(\widehat{x}\right),q\left(\widehat{y}\right))J\left(\widehat{y}\right),\widehat{x},\widehat{y}\in {\mathbb{S}}^2,\widehat{x}\ne \widehat{y}.}$$

The function K is a weakly singular integral kernel that can be rewritten in the form

$$K(\widehat{x},\widehat{y})={\displaystyle \frac{M(\widehat{x},\widehat{y})}{|\widehat{x}-\widehat{y}|}},\widehat{x},\widehat{y}\in {\mathbb{S}}^2,\widehat{x}\ne \widehat{y},$$

where

$$M(\widehat{x},\widehat{y})={\displaystyle \frac{e^{-k|q\left(\widehat{x}\right)-q\left(\widehat{y}\right)|}}{4\pi }}{\displaystyle \frac{(q\left(\widehat{x}\right)-q\left(\widehat{y}\right),\nu \left(q\left(\widehat{y}\right)\right))}{|q\left(\widehat{x}\right)-q\left(\widehat{y}\right){|}^2}}{\displaystyle \frac{|\widehat{x}-\widehat{y}|J\left(\widehat{y}\right)}{|q\left(\widehat{x}\right)-q\left(\widehat{y}\right)|}}\left(k|q\left(\widehat{x}\right)-q\left(\widehat{y}\right)|+1\right).$$

Note, due to the analyticality of q well-posedness of (23) also applies to (24). In order to discretize (24), we consider quadrature rules, proposed by Wienert [15]. For a given space discretization parameter $N\in \mathbb{N}$ the following values are defined

$${\widehat{p}}_{\beta ,\mu }=p\left({\vartheta }_{\beta },{\phi }_{\mu }\right),{\vartheta }_{\beta }=arccos{\xi }_{\beta },{\phi }_{\mu }={\displaystyle \frac{\pi }{N}}\mu ,$$

(25)

$$w_{\beta }^{\left(1\right)}={\displaystyle \frac{2\pi }{N(1-{\xi }_{\beta }^2){\left(P_N^{\prime }\left({\xi }_{\beta }\right)\right)}^2}},\beta =1,\dots ,N,\mu =0,\dots ,2N-1.$$

where ${\xi }_1<\dots <{\xi }_N$ are the zeros of the Legendre polynomials $P_N$ [17].

For a given function $f\in C({\mathbb{S}}^2\setminus \left\{(0,0,\pm 1)\right\})$ approximation $A_{N}f$ is defined as

$$\left(A_{N}f\right)\left(\widehat{x}\right)=\sum _{\beta =1}^N\sum _{\mu =0}^{2N-1}{w}_{\beta }^{\left(1\right)}f\left({\widehat{p}}_{\beta ,\mu }\right)a_{\beta ,\mu }\left(\widehat{x}\right),$$

(26)

where ${\displaystyle a_{\beta ,\mu }\left(\widehat{x}\right)=\sum _{n=0}^{N-1}\frac{2n+1}{4\pi }{P}_n(\widehat{x}·{\widehat{p}}_{\beta ,\mu })}$ and by $\widehat{x}·{\widehat{p}}_{\beta ,\mu }$ we denote a scalar product of two vectors. Note, the poles are excluded from the continuity requirements, since values of the parametrized functions $f\circ p$ at the poles may depend on the direction of approach (i.e. specific value of $\phi$) and may be not continous at the poles.

For the non-singular integrands, the following quadrature rule is suggested

$${\int }_{{\mathbb{S}}^2}f\left(\widehat{y}\right)ds\left(\widehat{y}\right)\approx {\int }_{{\mathbb{S}}^2}\left(A_{N}f\right)\left(\widehat{y}\right)ds\left(\widehat{y}\right)=\sum _{\beta =1}^N\sum _{\mu =0}^{2N-1}{w}_{\beta }^{\left(1\right)}f\left({\widehat{p}}_{\beta ,\mu }\right).$$

(27)

For the weakly singular integrands, the following quadrature rule can be used

$${\int }_{{\mathbb{S}}^2}{\displaystyle \frac{f\left(\widehat{y}\right)}{|{\widehat{n}}_p-\widehat{y}|}}ds\left(\widehat{y}\right)\approx {\int }_{{\mathbb{S}}^2}{\displaystyle \frac{\left(A_{N}f\right)\left(\widehat{y}\right)}{|{\widehat{n}}_p-\widehat{y}|}}ds\left(\widehat{y}\right)=\sum _{\beta =1}^N\sum _{\mu =0}^{2N-1}{w}_{\beta }^{\left(2\right)}f\left({\widehat{p}}_{\beta ,\mu }\right),$$

(28)

where ${\displaystyle w_{\beta }^{\left(2\right)}=w_{\beta }^{\left(1\right)}\sum _{n=0}^{N-1}{P}_n\left({\xi }_{\beta }\right)}$ and ${\widehat{n}}_p=(0,0,1).$

Both quadratures are obtained by approximation of the regular part of the integrand via approximation $A_N$ and then using exact integration. According to results in [15], these quadrature rules have super-algebraic or even exponential convergence order, depending on the smoothness of f.

By simple substitution, the quadrature rule (28) can be extended to a more general case

$${\int }_{{\mathbb{S}}^2}{\displaystyle \frac{f\left(\widehat{y}\right)}{|\widehat{x}-\widehat{y}|}}ds\left(\widehat{y}\right)={\int }_{{\mathbb{S}}^2}{\displaystyle \frac{f\left(T_{\widehat{x}}^{-1}\widehat{y}\right)}{|{\widehat{n}}_p-\widehat{y}|}}ds\left(\widehat{y}\right)\approx \sum _{\beta =1}^N\sum _{\mu =0}^{2N-1}{w}_{\beta }^{\left(2\right)}f\left(T_{\widehat{x}}^{-1}{\widehat{p}}_{\beta ,\mu }\right),\widehat{x}\in {\mathbb{S}}^2,$$

(29)

where $T_{\widehat{x}}$ is usually a rotation, such that $T_{\widehat{x}}\widehat{x}={\widehat{n}}_p$, see [15].

Applying (29) to the integral in (24), for $\tilde{\psi }\approx \psi$ we get an approximation equation

$$-{\displaystyle \frac{1}{2}}\tilde{\psi }\left(\widehat{x}\right)+\sum _{\beta =1}^N\sum _{\mu =0}^{2N-1}{w}_{\beta }^{\left(2\right)}\tilde{\psi }\left(T_{\widehat{x}}^{-1}{\widehat{p}}_{\beta ,\mu }\right)M(\widehat{x},T_{\widehat{x}}^{-1}{\widehat{p}}_{\beta ,\mu })={\widehat{G}}_s\left(\widehat{x}\right),\widehat{x}\in {\mathbb{S}}^2.$$

(30)

We observe that (30) contains values of the density $\tilde{\psi }$ in the rotated nodes $T_{\widehat{x}}^{-1}{\widehat{p}}_{\beta ,\mu }$. In order to be able to construct a system of linear equations, we replace $\tilde{\psi }\left(T_{\widehat{x}}^{-1}{\widehat{p}}_{\beta ,\mu }\right)$ with its approximation by $A_N$ (26)

$$\tilde{\psi }\left(T_{\widehat{x}}^{-1}{\widehat{p}}_{\beta ,\mu }\right)\approx \left(A_N\tilde{\psi }\right)\left(T_{\widehat{x}}^{-1}{\widehat{p}}_{\beta ,\mu }\right)=\sum _{{\beta }^{\prime }=1}^N\sum _{{\mu }^{\prime }=0}^{2N-1}{w}_{{\beta }^{\prime }}^{\left(1\right)}\tilde{\psi }\left({\widehat{p}}_{{\beta }^{\prime },{\mu }^{\prime }}\right)a_{{\beta }^{\prime },{\mu }^{\prime }}\left(T_{\widehat{x}}^{-1}{\widehat{p}}_{\beta ,\mu }\right).$$

(31)

Substituting (31) back into (30), we get

$$-{\displaystyle \frac{1}{2}}\tilde{\psi }\left(\widehat{x}\right)+\sum _{{\beta }^{\prime }=1}^N\sum _{{\mu }^{\prime }=0}^{2N-1}\tilde{\psi }\left({\widehat{p}}_{{\beta }^{\prime },{\mu }^{\prime }}\right)w_{{\beta }^{\prime },{\mu }^{\prime }}^{\left(3\right)}\left(\widehat{x}\right)={\widehat{G}}_s\left(\widehat{x}\right),\widehat{x}\in {\mathbb{S}}^2,$$

(32)

where ${\displaystyle w_{{\beta }^{\prime },{\mu }^{\prime }}^{\left(3\right)}\left(\widehat{x}\right)=w_{{\beta }^{\prime }}^{\left(1\right)}\sum _{\beta =1}^N\sum _{\mu =0}^{2N-1}{w}_{\beta }^{\left(2\right)}M(\widehat{x},T_{\widehat{x}}^{-1}{\widehat{p}}_{\beta ,\mu })a_{{\beta }^{\prime },{\mu }^{\prime }}\left(T_{\widehat{x}}^{-1}{\widehat{p}}_{\beta ,\mu }\right).}$

Collocating the equation (32) in the nodes ${\widehat{p}}_{{\beta }_2,{\mu }_2},{\beta }_2=1,\dots ,N,{\mu }_2=0,\dots ,2N-1$, we get a $2N^2\times 2N^2$ system of linear equations for the unknown values $\tilde{\psi }\left({\widehat{p}}_{{\beta }_2,{\mu }_2}\right)$

$$-{\displaystyle \frac{1}{2}}\tilde{\psi }\left({\widehat{p}}_{{\beta }_2,{\mu }_2}\right)+\sum _{{\beta }^{\prime }=1}^N\sum _{{\mu }^{\prime }=0}^{2N-1}\tilde{\psi }\left({\widehat{p}}_{{\beta }^{\prime },{\mu }^{\prime }}\right)w_{{\beta }^{\prime },{\mu }^{\prime }}^{\left(3\right)}\left({\widehat{p}}_{{\beta }_2,{\mu }_2}\right)={\widehat{G}}_s\left({\widehat{p}}_{{\beta }_2,{\mu }_2}\right).$$

(33)

After solving (33), the approximate solutions of problems (19)–(20) for the parameter $s=\{s_j,j=0,\dots ,N_t\}$ can be found by applying the quadrature rule (27) to (22)

$$U_{j,N}\left(x\right)=\sum _{{\beta }_2=1}^N\sum _{{\mu }_2=0}^{2N-1}{w}_{{\beta }_2}^{\left(1\right)}\widetilde{{\psi }_j}\left({\widehat{p}}_{{\beta }_2,{\mu }_2}\right){\Phi }_{\nu ,s_j}(x,q\left({\widehat{p}}_{{\beta }_2,{\mu }_2}\right))J\left({\widehat{p}}_{{\beta }_2,{\mu }_2}\right),x\in D,j=0,\dots ,N_t,$$

(34)

where $\widetilde{{\psi }_j}\left({\widehat{p}}_{{\beta }_2,{\mu }_2}\right)=\tilde{\psi }\left({\widehat{p}}_{{\beta }_2,{\mu }_2}\right)$ for the parameter value $s=s_j$.

Having solved a set of problems (19)–(20), we can construct the approximate solution of the original non-stationary problem (1)–(3)

$$u(x,t)\approx u_{N_t,N}(x,t)=\sum _{j=0}^{N_t}{\gamma }_j{e}^{t{s}_j}{U}_{j,N}\left(x\right)+\sum _{j=-N_t}^{-1}{\gamma }_j{e}^{t{s}_j}\overline{U_{-j,N}\left(x\right)}.$$

(35)

As mentioned, the error rate of the numerical inversion of the Laplace transform behaves like $\mathcal{O}\left(e^{-c{N}_t/ln{N}_t}\right)$ and is stable to perturbations of $U_j$ values. In our case these perturbations are created by the fact that $U_j$ are approximated by $U_{j,N}$, which in practice exhibits super-algebraic convergence rate for the sufficiently smooth surfaces and boundary conditions. As result, when $N_t$ and N are selected in the balanced way, the overall error rate of the original non-stationary problem is super-algebraic, which is shown by the numerical experiments.

4. Numerical Experiments

We will consider the following examples of the regions $D_k,k=1,2$ and their boundaries ${\Gamma }_k$ to perform numerical experiments (see Figure 4)

$${\Gamma }_1=\left\{r_1(\theta ,\phi )\left(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta \right),(\theta ,\varphi )\in [0,\pi ]\times [0,2\pi )\right\},$$

$$r_1(\theta ,\phi )=A_1\left(0.6+\sqrt{4.25+2cos3\theta }\right),A_1\in \mathbb{R}>0,$$

$${\Gamma }_2=\left\{r_2(\theta ,\phi )\left(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta \right),(\theta ,\varphi )\in [0,\pi ]\times [0,2\pi )\right\},$$

$$r_2(\theta ,\phi )=A_2\sqrt{0.8+0.2(cos(2\phi )-1)(cos(4\theta )-1)},A_2\in \mathbb{R}>0.$$

Note, when defining specific surfaces $\Gamma$ using mapping $q:{\mathbb{S}}^2\to \Gamma$, it is possible to use spherical or Cartesian coordinates to describe points on the unit sphere. For the mentioned surfaces ${\Gamma }_1$ and ${\Gamma }_2$ we used spherical coordinates.

4.1. Inverse Laplace Transform

Here we test the numerical inversion of the Laplace transform and suggested optimizations. Let us consider the fundamental solution of the heat equation (1)

$$\tilde{G}(x,y,t)={\displaystyle \frac{1}{\sqrt{4\pi t^3}}}{e}^{{\displaystyle \frac{{|x-y|}^2}{4t}}},x,y\in R^3,x\ne y,t>0,$$

for which the Laplace image is a fundamental solution of the Helmholtz equation (10)

$${\mathcal{L}}_t\left(\tilde{G}(x,y,t)\right)={\int }_0^{\infty }{e}^{-st}\tilde{G}(x,y,t)dt={\Phi }_s(x,y),$$

$${\mathcal{L}}_s^{-1}\left({\Phi }_s(x,y)\right)=\tilde{G}(x,y,t).$$

(36)

For a given source point $y^{*}\notin D_1,$ function

$$w(x,t)=\tilde{G}(x,y^{*},t),x\in D_1,y^{*}\notin D_1,$$

is an exact solution of the heat equation (1) and its Laplace transform $W_s\left(x\right)={\Phi }_s(x,y^{*})$ is an exact solution of the Hemholtz equation (10). Clearly, the Theorem 1 holds true in this case, so the values of $W_s$ can be computed only at the nodes $s_j,j=0,...,N_t$.

To test the effect of the $\alpha ,\lambda$ selection, we choose some random values of ${\alpha }_0,{\lambda }_0$ and compare the absolute error $E_{abs}=|T_{N_t}\left(W_s\right)-w|$ for the approximate computation of the inverse Laplace transform (36) to the absolute error $E_{abs,min}$ for the optimal parameters ${\alpha }_{min},{\lambda }_{min}$, obtained via optimization process (17).

Table 1 and Table 2 show the comparison of $E_{abs,min}$ and $E_{abs}$ for the different values of $x,y^{*},t$.

Obtained results support expected error rates. Note, the tested $({\alpha }_{min},{\lambda }_{min})$ search routine is computationally fast (involves $N_{\alpha }\times N_{\lambda }$ Laplace inversions) and is negligible compared to the computational effort of solving the stationary 3D problem. Comparing specific values in Table 1 and Table 2 one could expect significant reduction in necessary $N_t$, i.e. number of stationary problems to solve.

4.2. Stationary Problem

In this section we test the numerical solution of the stationary problems (10)–(11) using BIE method, described in section 3. As a sample boundary condition we choose $G_s$ as the narrowing of the fundamental solution ${\Phi }_s(·,y^{*})$ onto $\Gamma$ with a source point $y^{*}$ outside the region $D_k$.

In this case the exact solution of the problem (10)–(11) is

$$U_s^{*}\left(x\right)={\Phi }_s(x,y^{*}),x\in D_k,y^{*}\notin D_k,k=1,2.$$

Let $y^{*}=(0,0,5)\notin D_k,k=1,2$. To measure the accuracy of the numerical approximation, we use the following discrete $L_2$ error

$$E_{L_2}={\displaystyle \frac{1}{\tilde{N}}}{\left(\sum _{m=1}^{\tilde{N}}\sum _{j=1}^{\tilde{N}}{\left|U_{s_k}^{*}\left(x_{m,j}\right)-U_{k,N}\left(x_{m,j}\right)\right|}^2\right)}^{1/2},$$

where $U_{k,N}$ is the approximate solution obtained by (34) and $\tilde{N}=10$. The test points $x_{m,j}$ are uniformly distributed on a diminished artificial surface located within the solution domain, according to the following rule

$$x_{m,j}=0.5r_k({\theta }_m,{\varphi }_j)p({\theta }_m,{\varphi }_j),{\theta }_m={\displaystyle \frac{\pi }{\tilde{N}+1}}m,{\varphi }_j={\displaystyle \frac{2\pi }{\tilde{N}}}j,m,j=1,\dots ,\tilde{N},k=1,2.$$

Tables below show the error $E_{L_2}$ for the two test surfaces, different equation parameter values s and space discretization parameter N.

Table 3. Discrete $L_2$ error for the case $D=D_1$.

N	Nodes	${\mathit{s}}_1=1$	${\mathit{s}}_2=0.5-3\mathit{i}$	${\mathit{s}}_3=2+5\mathit{i}$
4	32	4.35e-06	2.41e-06	5.43e-08
8	128	6.12e-07	3.51e-07	3.26e-09
16	512	2.56e-09	6.80e-10	2.27e-11
32	2048	7.32e-12	3.66e-12	4.24e-12

Table 3. Discrete $L_2$ error for the case $D=D_1$.

N	Nodes	${\mathit{s}}_1=1$	${\mathit{s}}_2=0.5-3\mathit{i}$	${\mathit{s}}_3=2+5\mathit{i}$
4	32	4.35e-06	2.41e-06	5.43e-08
8	128	6.12e-07	3.51e-07	3.26e-09
16	512	2.56e-09	6.80e-10	2.27e-11
32	2048	7.32e-12	3.66e-12	4.24e-12

Table 4. Discrete $L_2$ error for the case $D=D_2$.

N	Nodes	${\mathit{s}}_1=1$	${\mathit{s}}_2=0.5-3\mathit{i}$	${\mathit{s}}_3=2+5\mathit{i}$
4	32	3.75e-06	8.92e-08	1.37e-08
8	128	4.74e-07	1.27e-08	2.43e-09
16	512	6.31e-09	5.83e-09	1.29e-10
32	2048	7.14e-11	4.29e-11	3.92e-11

Table 4. Discrete $L_2$ error for the case $D=D_2$.

N	Nodes	${\mathit{s}}_1=1$	${\mathit{s}}_2=0.5-3\mathit{i}$	${\mathit{s}}_3=2+5\mathit{i}$
4	32	3.75e-06	8.92e-08	1.37e-08
8	128	4.74e-07	1.27e-08	2.43e-09
16	512	6.31e-09	5.83e-09	1.29e-10
32	2048	7.14e-11	4.29e-11	3.92e-11

The obtained results support error rates, provided by Wienert [15].

4.3. Non-Stationary Problem

In this subsection we test the numerical solution of the non-stationary problem (1)–(3). The first example shows a case with an exactly known solution. The second example shows a case where exact solution is not known, but the boundary condition (3) has a Laplace transform available in closed form. For all examples, as suggested in section 2, we solve only $N_t+1$ stationary problems to provide approximate values $U_{j,N}\left(x\right)$ for the numerical inversion of the Laplace transform. We also apply parameters $(\alpha ,\lambda )$ selection technique, described in section 2 and tested in subsection 4.1.

4.3.1. Example with an Exactly Known Solution

As a sample boundary condition (20) we will choose narrowing of the fundamental solution (4.1) onto ${\Gamma }_1$

$$f(x,t)=\tilde{G}(x,y^{*},t),x\in {\Gamma }_1,y^{*}\notin D_1.$$

In this case the exact solution of the problem (1) - (3) is

$$u^{*}(x,t)=\tilde{G}(x,y^{*},t),x\in D_1,y^{*}\notin D_1.$$

Let $y^{*}=(0,0,5),x_0=(0.1,0.1,0.2)$. Table 5 shows the absolute error $E_{abs}^{\left(t\right)}=|u^{*}(x_0,t)-u_{N_t,N}(x_0,t)|$ of the approximate solution $u_{N_t,N}(x,t)$ for the different values of t and discretization parameters $N_t,N$.

Table 5. Absolute error $E_{abs}^{\left(t\right)}$, $\Gamma ={\Gamma }_1$.

${\mathit{N}}_{\mathit{t}}$	N	$\mathit{t}=2$	$\mathit{t}=2.2$	$\mathit{t}=2.5$
2	4	2.136812e-05	2.462481e-05	2.837609e-05
	8	1.946759e-06	2.188658e-06	2.331827e-06
	16	2.325619e-08	2.773015e-10	9.600292e-08
	32	1.540043e-08	8.910118e-09	8.651947e-08
4	4	2.137816e-05	2.461644e-05	2.844138e-05
	8	1.962045e-06	2.179749e-06	2.418113e-06
	16	7.838391e-09	8.531891e-09	9.433593e-09
	32	1.884036e-11	5.333509e-11	1.013018e-10
8	4	2.137815e-05	2.461638e-05	2.844134e-05
	8	1.962032e-06	2.179657e-06	2.418072e-06
	16	7.851940e-09	8.624379e-09	9.475018e-09
	32	5.286807e-12	8.816986e-12	1.190821e-11

Table 5. Absolute error $E_{abs}^{\left(t\right)}$, $\Gamma ={\Gamma }_1$.

${\mathit{N}}_{\mathit{t}}$	N	$\mathit{t}=2$	$\mathit{t}=2.2$	$\mathit{t}=2.5$
2	4	2.136812e-05	2.462481e-05	2.837609e-05
	8	1.946759e-06	2.188658e-06	2.331827e-06
	16	2.325619e-08	2.773015e-10	9.600292e-08
	32	1.540043e-08	8.910118e-09	8.651947e-08
4	4	2.137816e-05	2.461644e-05	2.844138e-05
	8	1.962045e-06	2.179749e-06	2.418113e-06
	16	7.838391e-09	8.531891e-09	9.433593e-09
	32	1.884036e-11	5.333509e-11	1.013018e-10
8	4	2.137815e-05	2.461638e-05	2.844134e-05
	8	1.962032e-06	2.179657e-06	2.418072e-06
	16	7.851940e-09	8.624379e-09	9.475018e-09
	32	5.286807e-12	8.816986e-12	1.190821e-11

To verify that obtained error rates agree with error rates of numerical Laplace inversion and stationary problems solution, next table highlights decimal exponents of absolute errors for the same point $x_0=(0.1,0.1,0.2)$ and $t=2$.

Table 6. Error rates comparison.

Laplace Inversion		Stationary		Non-Stationary
$N_t$	$lo{g}_{10}\left(E_{abs}\right)$	N	$lo{g}_{10}\left(E_{abs}\right)$	$N_t,N$	$lo{g}_{10}\left(E_{abs}\right)$
2	-8	16	-10	2, 16	-8
4	-11	32	-12	4, 32	-11
8	-15	32	-12	8, 32	-12

Table 6. Error rates comparison.

Laplace Inversion		Stationary		Non-Stationary
$N_t$	$lo{g}_{10}\left(E_{abs}\right)$	N	$lo{g}_{10}\left(E_{abs}\right)$	$N_t,N$	$lo{g}_{10}\left(E_{abs}\right)$
2	-8	16	-10	2, 16	-8
4	-11	32	-12	4, 32	-11
8	-15	32	-12	8, 32	-12

It is easy to see full discretization of the non-stationary problem results in error rate, defined by the worse error of Laplace inversion and stationary problem solution, which is expected. It also indicates a balanced selection of $N_t$ and N may provide best overall result.

4.3.2. Example Without Exactly Known Solution

Let us consider numerical solution of the non-stationary problem (1)–(3) with the following boundary condition

$$g(x,t)=t^2{e}^{-t},(x,t)\in {\Gamma }_1\times (0,\infty ).$$

(37)

For the stationary problems boundary condition we will use closed form of the Laplace transform of $g(x,t)$

$$G_s\left(x\right)={\mathcal{L}}_t\left(f(x,t)\right)={\displaystyle \frac{2}{{(s+1)}^3}},x\in {\Gamma }_1,s\in \mathbb{C},s\ne -1.$$

Table 7 shows values of the approximate solution of (1)–(3) for different combinations of $N_t,N$, different points $x\in D_1$ and time points t.

For each combination of $x,t$ we observe increasing number of same decimal digits as $N_t,N$ grow.

Let us consider region $D_1$ with scaling parameter of the boundary ${\Gamma }_1$ set to $A_1=6$. We will calculate and plot numerical solution of the non-stationary problem (1)–(3) with the boundary condition (37) over time interval $t\in (0;10]$ in two points $x_1=(0,0,0)$ and $x_2=(2.5,0,0)$. Note, two test points are intentionally selected so $x_1$ is placed "deeper" within the region and point $x_2$ is closer to the boundary. Parameter $A_1$ is intentionally selected to significantly scale up the region in order to observe time delay in propagation of the boundary condition behaviour inside the region. To produce the plots, we have chosen time step $\Delta t=0.2$ and numerically solved problem (1)–(3) 50 times for each test point using discretization parameters $N_t=4,N=8$. The entire computation process fits into 2-3 minutes using an average-level PC. This and previous numerical examples were performed using MATLAB.

As can be seen from the Figure 5, the peak value at $x_1$ is observed at a later time point compared to the peak at $x_2$, which agrees with the expected behaviour.

5. Conclusions

In this article, an effective combination of the Laplace transform and boundary integral equation method for the numerical solution of 3D initial boundary value problem for the heat equation was proposed and tested. Using observed and proven symmetry in the set of boundary value problems, computation work is reduced by almost a factor of 2. Additional optimization of the integration contour parameters was shown to further reduce the error. As mentioned in the article, due to the independence of the boudanry value problems, the proposed approach is also suitable for parallel computing.

In future work, the presented approach is planned to be tested with other methods for solving boundary value problems, as well as different equations or types of boundary conditions. Additionaly, further research of the optimization techniques for the numerical Laplace transform inversion is planned.