Approximation of Uncoupled Quasi-Static Thermoelasticity Solutions Based on Gaussians

Abstract

A fast approximation method to three dimensional equations in quasi-static uncoupled thermoelasticity is proposed. We approximate the density via Gaussian approximating functions introduced in the method approximate approximations. In this way the action of the integral operators on such functions is presented in a simple analytical form. If the density has separated representation, the problem is reduced to the computation of one-dimensional integrals which admit efficient cubature procedures. The comparison of the numerical and exact solution shows that these formulas are accurate and provide the predicted approximation rate $2,4,6$ and $8$.

1 Introduction

The equations of thermoelasticity describe the elastic and the thermal behavior of elastic, heat conductive media, in particular the reciprocal actions between elastic stresses and temperature differences. We consider the classical thermoelastic system where the elastic part is the usual second-order one in the space variable. In the static uncoupled thermoelasticity, thermal effects on a body are restricted to strains due to a steady-state temperature distribution. Uncoupled quasi-static thermoelasticity can be employed when slowly varying thermal and mechanical loads are encountered and dissipative effected can be neglected. The equations are a coupling of the equations of elasticity and of the heat equation ([2, p.76], [3])

$$\begin{aligned} \mu \, \Delta \, {{\textbf{u}}}+(\mathbf \lambda +\mu ){} & {} \textrm{grad}\,\textrm{div}\,{{\textbf{u}}}-\gamma \, \textrm{grad}\,\, T+\rho F=0 \end{aligned}$$

(1.1)

$$\begin{aligned}{} & {} \frac{\partial T}{\partial t}-\kappa \Delta \, T=0 \end{aligned}$$

(1.2)

$$\begin{aligned}{} & {} T({{\textbf{x}}},0)=g({{\textbf{x}}}) \end{aligned}$$

(1.3)

for $({{\textbf{x}}},t)\in {{\mathbb {R}}}^3\times [0,\infty )$, together with the corresponding initial and boundary conditions. The set of quantities $\mu , \lambda ,\gamma ,\rho ,\kappa$ are positive and $3\lambda +2\mu >0$. We suppose that $g:{{\mathbb {R}}}^3\rightarrow {{\mathbb {R}}}$, $F=(F_1,F_2,F_3):{{\mathbb {R}}}^3\times [0,\infty )\rightarrow {{\mathbb {R}}}^3$ with $g, F_1(\cdot ,t),F_2(\cdot ,t),F_3(\cdot ,t)\in {\mathscr {S}}({{\mathbb {R}}}^3)$. Here ${\mathscr {S}}({{\mathbb {R}}}^3)$ denotes the Schwartz space of smooth functions whose derivatives (including the function itself) decay at infinity faster than any power. The function $T({{\textbf{x}}},t)$ is the temperature and the vector ${{\textbf{u}}}({{\textbf{x}}},t)=(u_1({{\textbf{x}}},t),u_2({{\textbf{x}}},t), u_3({{\textbf{x}}},t))$ is the thermoelastic displacement.

The problem of determining $T({{\textbf{x}}},t)$ is an independent of ${{\textbf{u}}}({{\textbf{x}}},t)$ problem. The Cauchy problem (1.2)–(1.3) can be solved by the Poisson integral

$$\begin{aligned} T({{\textbf{x}}},t)=({\mathscr {P}}g)({{\textbf{x}}},t)=\frac{1}{(4\pi \kappa t)^{3/2}}\int \limits _{{{\mathbb {R}}}^3} \textrm{e}^{-\frac{|{{\textbf{x}}}-{{\textbf{y}}}|^2}{4\kappa t}} g({{\textbf{y}}})d{{\textbf{y}}}. \end{aligned}$$

(1.4)

We get

$$\begin{aligned} \textrm{grad}\,\, T({{\textbf{x}}},t)=\frac{1}{(4\pi \kappa t)^{3/2}}\int \limits _{{{\mathbb {R}}}^3} \frac{-2({{\textbf{x}}}-{{\textbf{y}}})}{4\kappa t} \textrm{e}^{-\frac{|{{\textbf{x}}}-{{\textbf{y}}}|^2}{4\kappa t}} g({{\textbf{y}}})d{{\textbf{y}}}. \end{aligned}$$

Moreover

$$\begin{aligned} \lim _{|{{\textbf{x}}}|\rightarrow \infty } T({{\textbf{x}}},t)=\lim _{|{{\textbf{x}}}|\rightarrow \infty }|\textrm{grad}\,\, T({{\textbf{x}}},t)|=0,\quad \forall t>0. \end{aligned}$$

When the temperature field $T$ is known, the displacement field ${{\textbf{u}}}=(u_1,u_2,u_3)$ is obtained by solving (1.1) where the gradient of $T$ is treated as a body force. The displacement field ${{\textbf{u}}}=(u_1,u_2,u_3)$ with $\lim _{|{{\textbf{x}}}|\rightarrow \infty } |{{\textbf{u}}}({{\textbf{x}}},t)|=0$, $t>0$, can be represented by means of the Kelvin fundamental matrix $\{\Gamma _{k\ell }\}_{k,\ell =1,2,3}$ ([4, p.84])

$$\begin{aligned} \Gamma _{k\ell }({{\textbf{x}}})=\frac{\lambda ' \delta _{k\ell } }{8\pi |{{\textbf{x}}}|}+\frac{\mu '}{8\pi } \frac{x_k x_\ell }{|{{\textbf{x}}}|^3},\quad k,\ell =1,2,3 \end{aligned}$$

(1.5)

with

$$\begin{aligned} \lambda '=\frac{\lambda +3\mu }{\mu (\lambda +2\mu )},\qquad \mu '=\frac{\lambda +\mu }{\mu (\lambda +2\mu )}. \end{aligned}$$

Hence, we have

$$\begin{aligned} u_k({{\textbf{x}}},t)=\sum _{\ell =1}^3 \int \limits _{{{\mathbb {R}}}^3} \Gamma _{k\ell }({{\textbf{x}}}-{{\textbf{y}}})(\rho F_\ell ({{\textbf{y}}},t) -\gamma \frac{\partial }{\partial y_\ell } T({{\textbf{y}}},t))d{{\textbf{y}}},\quad k=1,2,3. \end{aligned}$$

We write

$$\begin{aligned} {{\textbf{u}}}({{\textbf{x}}},t)={{\textbf{u}}}^{(1)}({{\textbf{x}}},t)+{{\textbf{u}}}^{(2)}({{\textbf{x}}},t), \end{aligned}$$

where ${{\textbf{u}}}^{(1)}({{\textbf{x}}},t)$ is the solution of

$$\begin{aligned} \mu \, \Delta {{\textbf{u}}}^{(1)}+(\mathbf \lambda +\mu ) \textrm{grad}\,\textrm{div}\,{{\textbf{u}}}^{(1)}+\rho F=0, \quad \lim _{|{{\textbf{x}}}|\rightarrow \infty } |{{\textbf{u}}}^{(1)}({{\textbf{x}}},t)|=0 \end{aligned}$$

(1.6)

and ${{\textbf{u}}}^{(2)}({{\textbf{x}}},t)$ is the solution of

$$\begin{aligned} \mu \, \Delta {{\textbf{u}}}^{(2)}+(\mathbf \lambda +\mu ) \textrm{grad}\,\textrm{div}\,{{\textbf{u}}}^{(2)}-\gamma \, \textrm{grad}\,\, T=0, \quad \lim _{|{{\textbf{x}}}|\rightarrow \infty } |{{\textbf{u}}}^{(2)}({{\textbf{x}}},t)|=0 \end{aligned}$$

(1.7)

with $T$ in (1.4).

The vectors ${{\textbf{u}}}^{(1)}=(u_1^{(1)},u_2^{(1)},u_3^{(1)})$ and ${{\textbf{u}}}^{(2)}=(u_1^{(2)},u_2^{(2)},u_3^{(2)})$ have the following integral representation by means of the Kelvin fundamental matrix

$$\begin{aligned}{} & {} u_k^{(1)}({{\textbf{x}}},t)=\rho \sum _{\ell =1}^3 \int \limits _{{{\mathbb {R}}}^3} \Gamma _{k\ell }({{\textbf{x}}}-{{\textbf{y}}}) F_\ell ({{\textbf{y}}},t) d{{\textbf{y}}},\quad k=1,2,3;\\{} & {} u_k^{(2)}({{\textbf{x}}},t)=-\gamma \sum _{\ell =1}^3 \int \limits _{{{\mathbb {R}}}^3} \Gamma _{k\ell }({{\textbf{x}}}-{{\textbf{y}}})\frac{\partial }{\partial y_\ell } T({{\textbf{y}}},t)d{{\textbf{y}}},\quad k=1,2,3. \end{aligned}$$

Fast formulas of high order for the approximation of ${{\textbf{u}}}^{(1)}$ were obtained in [9]. The goal of this paper is to derive semi-analytic cubature formulas for $({{\textbf{u}}}^{(2)},T)$ solutions to (1.7)–(1.2)–(1.3) of an arbitrary high-order which are fast and accurate by using the basis functions introduced in the theory approximate approximations ([11, 12]; see also [15] and the reference therein).

The approximate quasi-interpolant has the form

$$\begin{aligned} {\mathscr {M}}_{h,{\mathscr {D}}}{g}({{\textbf{x}}})={\mathscr {D}}^{-3/2} \sum _{{{\textbf{m}}}\in {{\mathbb {Z}}}^3}{g}(h{{\textbf{m}}})\eta \left( \frac{{{\textbf{x}}}-h{{\textbf{m}}}}{h\sqrt{{\mathscr {D}}}}\right) \end{aligned}$$

(1.8)

where h and ${\mathscr {D}}$ are positive parameters and $\eta$ is a smooth and rapidly decaying function which satisfies the moment conditions of order N

$$\begin{aligned} \int \limits _{{{\mathbb {R}}}^3}\eta ({{\textbf{x}}})\,{{\textbf{x}}}^{\alpha }d{{\textbf{x}}}=\delta _{0,\alpha },\quad 0\le |\alpha |<N. \end{aligned}$$

(1.9)

If we define the Fourier transform of $\eta$ as

$$\begin{aligned} {\mathscr {F}}\eta ({{\textbf{x}}})=\int \limits _{{{\mathbb {R}}}^3} \eta ({{\textbf{y}}}) \textrm{e}^{-2 i \pi \langle {{\textbf{x}}}, {{\textbf{y}}}\rangle }d{{\textbf{y}}}, \end{aligned}$$

then following [15, p.34] the approximate quasi-interpolant can be written in the form

$$\begin{aligned} {\mathscr {M}}_{h,{\mathscr {D}}}{g}({{\textbf{x}}})&= g({{\textbf{x}}}) +(-\sqrt{{\mathscr {D}}}h)^N g_N({{\textbf{x}}})\nonumber \\&\quad +\sum _{ |\varvec{\alpha }| = 0} ^{N-1} \frac{(\sqrt{{\mathscr {D}}}h)^{|\varvec{\alpha }|}}{\varvec{\alpha } !(2 \pi i)^{|{\varvec{\alpha }}|}} \, \partial ^{\varvec{\alpha }} g({{\textbf{x}}}) \, \rho _{\varvec{\alpha }} \big (\frac{{{\textbf{x}}}}{h},\eta ,{\mathscr {D}}\big ) \end{aligned}$$

(1.10)

with the function

$$\begin{aligned} g_N({{\textbf{x}}})= {{\mathscr {D}}}^{-3/2} \sum _{ |\varvec{\alpha }| = N} \frac{N}{\varvec{\alpha }!} \! \sum _{{{\textbf{m}}}\in {{{\mathbb {Z}}}}^3} \Big (\frac{ {{\textbf{x}}}-h {{\textbf{m}}}}{h\sqrt{{\mathscr {D}}}}\Big )^{\varvec{\alpha }} \! \eta \Big ( \frac{{{\textbf{x}}}\!-\!h{{\textbf{m}}}}{h\sqrt{{\mathscr {D}}}} \Big ) \int \limits _0 ^1 s^{N-1} \partial ^{\varvec{\alpha }} g (s {{\textbf{x}}}+ (1-s) h{{\textbf{m}}}) \, ds \end{aligned}$$

containing the remainder of the Taylor expansion of g. The functions

$$\begin{aligned} \rho _{\varvec{\alpha }} \big (\frac{{{\textbf{x}}}}{h},\eta ,{\mathscr {D}}\big ) =\sum _{\varvec{\nu } \in {{{\mathbb {Z}}}}^3\setminus \{0\}} \partial ^{\varvec{\alpha }}{\mathscr {F}}\eta (\sqrt{{\mathscr {D}}} {\varvec{\nu }}) \, \textrm{e}^{ \frac{2 \pi i}{h} \langle {{\textbf{x}}}, {\varvec{\nu }}\rangle } \end{aligned}$$

(1.11)

are rapidly oscillating multivariate trigonometric series and

$$\begin{aligned} | \rho _{\varvec{\alpha }}({{\textbf{x}}},\eta ,{\mathscr {D}})|\le \sum _{\varvec{\nu } \in {{{\mathbb {Z}}}}^3 \setminus \{0\}} \big | \partial ^{\varvec{\alpha }}{\mathscr {F}}\eta (\sqrt{{\mathscr {D}}} {\varvec{\nu }})\big | \end{aligned}$$

(1.12)

uniformly in ${{\textbf{x}}}$. Denoting

$$\begin{aligned} \varepsilon _k ({\mathscr {D}}) = \max _{|\varvec{\alpha }| = k} \sum _{\varvec{\nu } \in {{{\mathbb {Z}}}}^3 \setminus \{0\}} \big | \partial ^{\varvec{\alpha }}{\mathscr {F}}\eta (\sqrt{{\mathscr {D}}} {\varvec{\nu }})\big | \end{aligned}$$

we derive

$$\begin{aligned} \left| \sum _{ |\varvec{\alpha }| = 0} ^{N-1} \frac{(h \sqrt{{\mathscr {D}}})^{|\varvec{\alpha }|}}{\varvec{\alpha } ! (2 \pi i)^{|{\varvec{\alpha }}|}} \, \partial ^{\varvec{\alpha }} g({{\textbf{x}}}) \, \rho _{\varvec{\alpha }} \left( \frac{{{\textbf{x}}}}{h},\eta ,{\mathscr {D}}\right) \right| \le \sum _{ k = 0} ^{N-1}\varepsilon _k ({\mathscr {D}}) \frac{(\sqrt{{\mathscr {D}}}h)^{k}}{(2 \pi )^{k}} \sum _{ |\varvec{\alpha }| = k}|\partial ^{\varvec{\alpha }} g({{\textbf{x}}}) | \, . \end{aligned}$$

Thus, at any point ${{\textbf{x}}}$ we have

$$\begin{aligned} |{g}({{\textbf{x}}})-{\mathscr {M}}_{h,{\mathscr {D}}}{g}({{\textbf{x}}})| \le c (\sqrt{{\mathscr {D}}}h)^N \Vert \nabla _N g\Vert _{L_\infty } + \sum _{k=0}^{N-1}\frac{ \varepsilon _k({\mathscr {D}}) }{(2\pi )^k}(\sqrt{{\mathscr {D}}}h)^k \left| \nabla _k g({{\textbf{x}}})\right| , \end{aligned}$$

(1.13)

where $\nabla _k g$ denotes the vector of partial derivatives $\{\partial ^{\alpha }g \}_{|{\alpha }=k}$. The second term in the right hand side of (1.13) is called the saturation error.

Since $\eta \in {\mathscr {S}}({{\mathbb {R}}}^3)$ implies $\varepsilon _k ({\mathscr {D}}) \rightarrow 0$ as ${\mathscr {D}}\rightarrow \infty$ a proper choice of the parameter ${\mathscr {D}}$ allows to make the terms $\varepsilon _k({\mathscr {D}})$ as small as necessary, for example less than the machine precision. Therefore, the quasi-interpolant ${\mathscr {M}}_{h,{\mathscr {D}}}g$ can behave in numerical computations like a converging approximation process. Similar estimates hold in integral norms.

Theorem 1.1

[15, p.42] Suppose that $\eta \in {\mathscr {S}}({{\mathbb {R}}}^3)$ satisfies the moment condition (1.9). Then for any $g\in W^L_p({{\mathbb {R}}}^3)$, $1\le p\le \infty$ and $L>3/p$, $L\ge N$, the quasi-interpolant (1.8) satisfies

$$\begin{aligned} ||g-{\mathscr {M}}_{h,{\mathscr {D}}}g||_{L_p}\le c_\eta (\sqrt{{\mathscr {D}}}h)^{N} ||\nabla _N g||_{L_p}+\sum _{k=0}^{N-1}\frac{\varepsilon _k({\mathscr {D}})}{(2\pi )^k} (\sqrt{{\mathscr {D}}}h)^k ||\nabla _k g||_{L_p} \end{aligned}$$

(1.14)

where the constant $c_\eta$ does not depend on $g$, $h$ and ${\mathscr {D}}$.

New classes of cubature formulas for important integral operators of mathematical physics by using approximate approximations were studied in [14]. They are based on replacing the density of the integral operator by its quasi-interpolant where the generating function $\eta$ is chosen such that the operator applied to it can be computed, analytically or at least efficiently. We choose as basis functions products of Gaussians and special polynomials. The use of the Gaussian functions for the numerical solution of the problems under consideration has the main advantage that the action of the integral operators on such functions may be presented in a simple analytical form.

By combining cubature formulas for volume potentials based on approximate approximations with the strategy of separated representations (cf., e.g. [1]), it is possible to derive a method for approximating volume potentials which is accurate and fast also in the multidimensional case and provides approximation formulas of high order. This procedure was applied successfully for the first time to the integration of the harmonic potential [5]. This approach was extended to the biharmonic [7], elastic and hydrodynamic [9] potentials, and to parabolic problems [6]. New approximation formulas for the solutions of nonstationary Stokes system were obtained in [8]. The static thermoelasticity was considered in [10]. Here we show that the fast method can be applied to uncoupled quasi-static thermoelasticity.

The outline of the paper is the following. In Sect. 2 we describe the fast formulas for the approximation of $T$ obtained in [6]. In Sect. 3 we use the approximants obtained in Sect. 2 to construct approximation formulas for ${{\textbf{u}}}^{(2)}$ and give error estimates. In Sect. 4 we provide results of numerical experiments, illustrating that our formulas are accurate and provide the predicted approximation rates $2$, $4$, $6$ and $8$.

2 Approximation of T

Cubature formulas for (1.4) are derived by replacing the density $g$ with the quasi-interpolant (1.8). Then

$$\begin{aligned} ({\mathscr {P}}{\mathscr {M}}_{h,{\mathscr {D}}}g)({{\textbf{x}}},t)=\frac{1}{{\mathscr {D}}^{3/2}}\sum _{{{\textbf{m}}}\in {{\mathbb {Z}}}^3} g(h{{\textbf{m}}}) ({\mathscr {P}}\eta )\left( \frac{{{\textbf{x}}}-h{{\textbf{m}}}}{h\sqrt{{\mathscr {D}}}},\frac{t}{h^2{\mathscr {D}}}\right) \end{aligned}$$

(2.1)

provides an approximation formula for $T({{\textbf{x}}},t)$.

The cubature error can be estimated by the following.

Theorem 2.1

[15, Theorem 6.1] Suppose that $\eta$ satisfies the moment condition (1.9). If the initial values of the parabolic problem (1.2)–(1.3) satisfy $g\in W^{N}_p({{\mathbb {R}}}^3)$, $1\le p\le \infty$, then the approximate solution (2.1) converges for any fixed $t>0$ with the order ${\mathscr {O}}(h^N)$ to the solution of the problem.

As basis functions in (1.8) we take the tensor products of univariate basis functions

$$\begin{aligned} {\eta _{2M}}({{\textbf{x}}})=\prod _{j=1}^{3} {{\widetilde{\eta }}}_{2M}(x_j); \quad {\widetilde{\eta }}_{2M}(x)=\frac{(-1)^{M-1}}{\sqrt{\pi }2^{2M-1} (M-1)!}\frac{H_{2M-1}(x) \textrm{e}^{-x^2}}{x} , \end{aligned}$$

(2.2)

where $H_k$ are the Hermite polynomials

$$\begin{aligned} H_k(x)=(-1)^k \textrm{e}^{x^2} \left( \frac{d}{dx}\right) ^k \textrm{e}^{-x^2}. \end{aligned}$$

Theorem 2.2

Let $M\ge 1$. The Poisson integral (1.4) applied to the generating functions $\eta _{2M}$ in (2.2) can be written as

$$\begin{aligned} ({\mathscr {P}}\eta _{2M})({{\textbf{x}}},t)=\prod _{j=1}^3 \frac{ {\mathscr {Q}}_M(x_j,4 \kappa t) }{\sqrt{1+4\kappa t}}\textrm{e}^{-|{{\textbf{x}}}|^2/(1+4\kappa t)} \end{aligned}$$

(2.3)

where ${\mathscr {Q}}_M(x,t)$ is a polynomial in $x$ of degree $2\,M-2$ whose coefficients depend on $t$, defined by

$$\begin{aligned} {\mathscr {Q}}_M(x,t)= \frac{1}{\sqrt{\pi }}\sum _{s=0}^{M-1} \frac{1}{(1+ t)^{s}} \frac{(-1)^s}{4^s s!}H_{2s}\left( \frac{x}{\sqrt{1+t}}\right) . \end{aligned}$$

(2.4)

Proof

We have

$$\begin{aligned} \left( {\mathscr {P}}\left( \prod _{j=1}^3 \widetilde{\eta }_{2M}\right) \right) ({{\textbf{x}}}, t)=\prod _{j=1}^3 \frac{1}{\sqrt{4\pi \kappa t}} \int \limits _{-\infty }^\infty \textrm{e}^{-(x_j-y_j)^2/(4\kappa t)} \widetilde{\eta }_{2M}(y_j)\, dy_j. \end{aligned}$$

Using the representation ([15, p.55])

$$\begin{aligned} {{\widetilde{\eta }}}_{2M}(x)=A\left( \frac{d}{dx}\right) \textrm{e}^{-x^2}, \quad A\left( \frac{d}{dx}\right) =\frac{1}{\sqrt{\pi }} \sum _{s=0}^{M-1}\frac{(-1)^s}{s! 4^s}\frac{d^{2s}}{dx^{2s}} \end{aligned}$$

(2.5)

and the relation

$$\begin{aligned} \int \limits _{-\infty }^\infty \textrm{e}^{-\frac{(x-y)^2}{a}} \textrm{e}^{-\frac{(y-z)^2}{b}} dy= \left( \frac{\pi a b}{a+b}\right) ^{1/2} e^{-\frac{(x-z)^2}{a+b}},\quad a>0,\quad b>0, \end{aligned}$$

(2.6)

we get

$$\begin{aligned} \int \limits _{-\infty }^\infty \textrm{e}^{-\frac{(x-y)^2}{4\kappa t}} \widetilde{\eta }_{2M}(y)\, dy =A\left( \frac{d}{dx}\right) \int \limits _{-\infty }^\infty \textrm{e}^{-\frac{(x-y)^2}{4\kappa t}}\textrm{e}^{-y^2}\, dy =\left( \frac{4\pi \kappa t}{1+4\kappa t}\right) ^{1/2}A\left( \frac{d}{dx}\right) \textrm{e}^{-\frac{x^2}{1+4\kappa t}}. \end{aligned}$$

By direct computation, the polynomials ${\mathscr {Q}}_M$ satisfy

$$\begin{aligned} A\left( \frac{d}{dx}\right) \textrm{e}^{-\frac{x^2}{1+t}}={\mathscr {Q}}_M(x,t)\textrm{e}^{-\frac{x^2}{1+t}}. \end{aligned}$$

(2.7)

Formula (2.3) easily follows. $\square$

Using formula (2.3), we can specify the high order approximation $T^{(M)}_{h,{\mathscr {D}}}({{\textbf{x}}},t):=({\mathscr {P}}{\mathscr {M}}_{h,{\mathscr {D}}}g)({{\textbf{x}}},t)$ as follows

$$\begin{aligned} T^{(M)}_{h,{\mathscr {D}}}({{\textbf{x}}},t)=\sum _{{{\textbf{m}}}\in {{\mathbb {Z}}}^3} g(h{{\textbf{m}}}) \frac{\textrm{e}^{-\frac{|{{\textbf{x}}}-h{{\textbf{m}}}|^2}{h^2{\mathscr {D}}(1+4\kappa t/(h^2{\mathscr {D}}))}}}{({\mathscr {D}}(1+4\kappa t/(h^2{\mathscr {D}})))^{3/2}}\prod _{j=1}^3 {\mathscr {Q}}_M\left( \frac{x_j-hm_j}{h\sqrt{{\mathscr {D}}}},\frac{4\kappa t}{h^2{\mathscr {D}}}\right) \end{aligned}$$

(2.8)

for the generating function $\eta _{2M}$ defined in (2.2). This is a semi-analytic cubature formula for (1.4) with the error ${\mathscr {O}}(h^{2M})$.

From (2.8), at the points $(h{{\textbf{s}}},t)$, ${{\textbf{s}}}=(s_1,s_2,s_3)\in {{\mathbb {Z}}}^3$,

$$\begin{aligned} T^{(M)}_{h,{\mathscr {D}}}(h{{\textbf{s}}},t)=\sum _{{{\textbf{m}}}\in {{\mathbb {Z}}}^3} g(h{{\textbf{m}}}) \frac{\textrm{e}^{-\frac{|{{\textbf{s}}}-{{\textbf{m}}}|^2}{{\mathscr {D}}(1+4\kappa t/(h^2{\mathscr {D}}))}}}{({\mathscr {D}}(1+4\kappa t/(h^2{\mathscr {D}})))^{3/2}}\prod _{j=1}^3 {\mathscr {Q}}_M\left( \frac{s_j-m_j}{\sqrt{{\mathscr {D}}}},\frac{4\kappa t}{h^2{\mathscr {D}}}\right) . \end{aligned}$$

(2.9)

Remark 2.3

The polynomials ${\mathscr {Q}}_M$ for $M=1,2,3,4$ are given by

$$\begin{aligned} {\mathscr {Q}}_1(x,t)= & {} 1/\sqrt{\pi },\\ {\mathscr {Q}}_2(x,t)= & {} \frac{1}{\sqrt{\pi }}\left( -\frac{x^2}{(t+1)^2}+\frac{1}{2 (t+1)}+1\right) ,\\ {\mathscr {Q}}_3(x,t)= & {} {\mathscr {Q}}_2(x,t)+\frac{1}{\sqrt{\pi }}\left( \frac{x^4}{2 (t+1)^4}-\frac{3 x^2}{2 (t+1)^3}+\frac{3}{8 (t+1)^2}\right) ,\\ {\mathscr {Q}}_4(x,t)= & {} {\mathscr {Q}}_3(x,t)+ \frac{1}{\sqrt{\pi }}\left( -\frac{x^6}{6 (t+1)^6}+\frac{5 x^4}{4 (t+1)^5}-\frac{15 x^2}{8 (t+1)^4}+\frac{5}{16 (t+1)^3} \right) . \end{aligned}$$

The approximation formulas (2.9) are very efficient if $g$ has a separated representation, i.e. for a given accuracy $\varepsilon$ it can be represented as the sum of products of vectors in dimension $1$

$$\begin{aligned} g({{\textbf{x}}})= \sum _{\ell =1}^L \prod _{r=1}^3g_r^{(\ell )}(x_r)+{\mathscr {O}}(\varepsilon ). \end{aligned}$$

(2.10)

Then $T(h{{\textbf{s}}},t)$ can be approximated by the sum of products of one-dimensional sums

$$\begin{aligned} T^{(M)}_{h,{\mathscr {D}}} g(h{{\textbf{s}}},t)=\sum _{\ell =1}^L \prod _{r=1}^3 S_r^{(\ell )} \left( s_r,{4\kappa t} \right) \end{aligned}$$

where

$$\begin{aligned} S_r^{(\ell )}\left( s,t\right) =\sum _{m\in {{\mathbb {Z}}}} g_r^{(\ell )}(hm) \frac{\textrm{e}^{-\frac{(s-m)^2}{{\mathscr {D}}(1+t)}}}{({\mathscr {D}}(1+t))^{1/2}} {\mathscr {Q}}_M\left( \frac{s-m}{\sqrt{{\mathscr {D}}}}, \frac{t}{h^2{\mathscr {D}}}\right) , \quad r=1,2,3. \end{aligned}$$

3 Approximation of ${{\textbf{u}}}^{(2)}{}$

In this section we propose formulas for the approximation of

$$\begin{aligned} u_k^{(2)}({{\textbf{x}}},t)=-\gamma \sum _{\ell =1}^3 \int \limits _{{{\mathbb {R}}}^3} \Gamma _{k\ell }({{\textbf{x}}}-{{\textbf{y}}})\frac{\partial }{\partial y_\ell } T({{\textbf{y}}},t)d{{\textbf{y}}},\quad k=1,2,3 \end{aligned}$$

where $T$ is given in (1.4).

Integrating by parts and using the relation

$$\begin{aligned} \frac{\partial }{\partial y_\ell } \Gamma _{k\ell }({{\textbf{x}}}-{{\textbf{y}}})=-\frac{\partial }{\partial x_\ell } \Gamma _{k\ell }({{\textbf{x}}}-{{\textbf{y}}}) \end{aligned}$$

we get

$$\begin{aligned} u_k^{(2)}({{\textbf{x}}},t) =-\gamma \sum _{\ell =1}^3 \int \limits _{{{\mathbb {R}}}^3}\frac{\partial }{\partial x_\ell } \Gamma _{k\ell }({{\textbf{x}}}-{{\textbf{y}}}) T({{\textbf{y}}},t)d{{\textbf{y}}},\quad k=1,2,3. \end{aligned}$$

From the relation [4, p.84]

$$\begin{aligned} \sum _{\ell =1}^3 \frac{\partial }{\partial x_\ell } \Gamma _{k\ell }({{\textbf{x}}})=\frac{1}{4\pi (\lambda +2\mu )}\frac{\partial }{\partial x_k}\frac{1}{|{{\textbf{x}}}|} \end{aligned}$$

we obtain

$$\begin{aligned} u_k^{(2)}({{\textbf{x}}},t)=-\frac{c_{\gamma ,\lambda +2\mu }}{4\pi }\frac{\partial }{\partial x_k}\int \limits _{{{\mathbb {R}}}^3}\frac{T({{\textbf{y}}},t)}{|{{\textbf{x}}}-{{\textbf{y}}}|}d{{\textbf{y}}},\quad k=1,2,3 \end{aligned}$$

(3.1)

where we set

$$\begin{aligned} c_{\gamma ,\lambda +2\mu }=\frac{\gamma }{\lambda +2\mu }. \end{aligned}$$

Since $T({{\textbf{y}}},t)=({\mathscr {P}}g)({{\textbf{y}}},t)$ we can also write

$$\begin{aligned} u_k^{(2)}({{\textbf{x}}},t)=-c_{\gamma ,\lambda +2\mu }\frac{\partial }{\partial x_k}{\mathscr {L}}(({\mathscr {P}}g)(\cdot ,t))({{\textbf{x}}}),\quad k=1,2,3 \end{aligned}$$

(3.2)

where we denote by ${\mathscr {L}}$ the harmonic potential

$$\begin{aligned} {\mathscr {L}}(g)({{\textbf{x}}})=\frac{1}{4\pi } \int \limits _{{{\mathbb {R}}}^3} \frac{g({{\textbf{y}}})}{|{{\textbf{x}}}-{{\textbf{y}}}|}d{{\textbf{y}}}. \end{aligned}$$

(3.3)

We use the representation (1.4) to get

$$\begin{aligned} u_k^{(2)}({{\textbf{x}}},t)=-\frac{c_{\gamma ,\lambda +2\mu }}{4\pi }\frac{1}{(4\pi \kappa t)^{3/2}} \frac{\partial }{\partial x_k}\int \limits _{{{\mathbb {R}}}^3}\frac{d{{\textbf{y}}}}{|{{\textbf{x}}}-{{\textbf{y}}}|}\int \limits _{{{\mathbb {R}}}^3} \textrm{e}^{-\frac{|{{\textbf{y}}}-{{\textbf{z}}}|^2}{4\kappa t}} g({{\textbf{z}}})d{{\textbf{z}}},\quad k=1,2,3 \end{aligned}$$

and we change the order of integration

$$\begin{aligned}{} & {} u_k^{(2)}({{\textbf{x}}},t)=-\frac{c_{\gamma ,\lambda +2\mu }}{4\pi }\frac{1}{(4\pi \kappa t)^{3/2}}\frac{\partial }{\partial x_k}\int \limits _{{{\mathbb {R}}}^3}g({{\textbf{z}}})d{{\textbf{z}}}\int \limits _{{{\mathbb {R}}}^3}\frac{ \textrm{e}^{-\frac{|{{\textbf{y}}}-{{\textbf{z}}}|^2}{4\kappa t}}}{|{{\textbf{x}}}-{{\textbf{y}}}|} d{{\textbf{y}}}\\{} & {} \quad =-\frac{c_{\gamma ,\lambda +2\mu }}{4\pi }\frac{1}{\pi ^{3/2}}\frac{1}{(4\kappa t)^{1/2}}\frac{\partial }{\partial x_k}\int \limits _{{{\mathbb {R}}}^3}g({{\textbf{z}}})d{{\textbf{z}}}\int \limits _{{{\mathbb {R}}}^3}\frac{\textrm{e}^{-|{{\textbf{w}}}|^2}}{\left| \frac{{{\textbf{x}}}-{{\textbf{z}}}}{\sqrt{4\kappa t}}-{{\textbf{w}}}\right| }d{{\textbf{w}}},\quad k=1,2,3. \end{aligned}$$

Then

$$\begin{aligned} u_k^{(2)}({{\textbf{x}}},t)=-\frac{c_{\gamma ,\lambda +2\mu }}{\pi ^{3/2}}\frac{1}{(4\kappa t)^{1/2}}\frac{\partial }{\partial x_k}\int \limits _{{{\mathbb {R}}}^3}{\mathscr {L}}(\textrm{e}^{-|\cdot |^2})(\frac{{{\textbf{x}}}-{{\textbf{z}}}}{\sqrt{4\kappa t}})g({{\textbf{z}}})d{{\textbf{z}}},\> k=1,2,3. \end{aligned}$$

(3.4)

We use the representation ([15, p.128])

$$\begin{aligned} {\mathscr {L}}(\textrm{e}^{-|\cdot |^2})({{\textbf{x}}})=\frac{1}{4} \int \limits _0^\infty \frac{\textrm{e}^{-\frac{|{{\textbf{x}}}|^2}{1+\tau }}}{(1+\tau )^{3/2}}d\tau \end{aligned}$$

to get

$$\begin{aligned} u_k^{(2)}({{\textbf{x}}},t)=-\frac{c_{\gamma ,\lambda +2\mu }}{4}\frac{1}{\pi ^{3/2}}\frac{1}{(4\kappa t)^{1/2}}\frac{\partial }{\partial x_k} \int \limits _0^\infty \frac{d\tau }{(1+\tau )^{3/2}} \int \limits _{{{\mathbb {R}}}^3} \textrm{e}^{-\frac{|\frac{{{\textbf{x}}}-{{\textbf{z}}}}{\sqrt{4\kappa t}}|^2}{1+\tau }} g({{\textbf{z}}})d{{\textbf{z}}}. \end{aligned}$$

(3.5)

Now we replace $g$ in (3.5) by the approximate quasi-interpolant (1.8) and we set

$$\begin{aligned}{} & {} ({\mathscr {N}}_{h,{\mathscr {D}}} g)_k({{\textbf{x}}},t):=\nonumber \\{} & {} -\frac{c_{\gamma ,\lambda +2\mu }}{4\pi ^{3/2}}\frac{h^3}{(4\kappa t)^{1/2}} \sum _{{{\textbf{m}}}\in {{\mathbb {Z}}}^3} g(h{{\textbf{m}}}) \frac{\partial }{\partial x_k} \int \limits _0^\infty \frac{d\tau }{(1+\tau )^{3/2}} \int \limits _{{{\mathbb {R}}}^3} \textrm{e}^{- \frac{|r_{{\textbf{m}}}({{\textbf{x}}})-{{\textbf{w}}}|^2}{(1+\tau )\frac{4\kappa t}{h^2{\mathscr {D}}}}} \eta \left( {{\textbf{w}}}\right) d{{\textbf{w}}}\end{aligned}$$

(3.6)

with

$$\begin{aligned} r_{{\textbf{m}}}({{\textbf{x}}})=\frac{{{\textbf{x}}}-h{{\textbf{m}}}}{h\sqrt{{\mathscr {D}}}}. \end{aligned}$$

In the next theorem we estimate the error of the cubature formula ${\mathscr {N}}_{h,{\mathscr {D}}} g$.

Theorem 3.1

Suppose that $\eta$ satisfies the moment condition (1.9). Let $1<p<3$, $q=3p/(3-p)$, and let $g\in W_p^L({{\mathbb {R}}}^3)$ with $L>3/p$, $L\ge N$. Then there exist two constants $c$ and $C$ such that, for any fixed $t>0$

$$\begin{aligned} ||({\mathscr {N}}_{h,{\mathscr {D}}} g)(\cdot ,t))-{{\textbf{u}}}^{(2)}(\cdot ,t)||_{L_q}\le c (h\sqrt{{\mathscr {D}}})^N||\nabla _Ng||_{L_p}+ C\,h^N. \end{aligned}$$

(3.7)

The constant $c$ does not depend on $g$, $h$ and ${\mathscr {D}}$ and $C{ isindependentof}h$.

Proof

Since $({\mathscr {N}}_{h,{\mathscr {D}}} g)({{\textbf{x}}},t)=-c_{\gamma ,\lambda +2\mu }\nabla {\mathscr {L}}(({\mathscr {P}}{\mathscr {M}}_{h,{\mathscr {D}}}g)(\cdot ,t))({{\textbf{x}}})$ and ${\textbf {u}}^{(2)}({{\textbf{x}}},t)=-c_{\gamma ,\lambda +2\mu }\nabla {\mathscr {L}}(({\mathscr {P}}g)(\cdot ,t))({{\textbf{x}}})$, we have to estimate the difference

$$\begin{aligned} ||\nabla {\mathscr {L}}(({\mathscr {P}}g)(\cdot ,t))-\nabla {\mathscr {L}}(({\mathscr {P}}{\mathscr {M}}_{h,{\mathscr {D}}}g)(\cdot ,t))||_{L_q}. \end{aligned}$$

(3.8)

Since

$$\begin{aligned} \nabla {\mathscr {L}}(({\mathscr {P}}g)(\cdot ,t))-\nabla {\mathscr {L}}(({\mathscr {P}}{\mathscr {M}}_{h,{\mathscr {D}}}g)(\cdot ,t))=\nabla {\mathscr {L}}({\mathscr {P}}(g- {\mathscr {M}}_{h,{\mathscr {D}}}g)(\cdot ,t)), \end{aligned}$$

the norm $||\nabla u||_{L_q}$ is equivalent to the norm $||(-\Delta )^{1/2} u||_{L_q}$ ([13, p.458]) and ${\mathscr {L}}$ is the inverse of the Laplacian, we obtain

$$\begin{aligned} ||(-\Delta )^{1/2}({\mathscr {L}}({\mathscr {P}}g)(\cdot ,t)-{\mathscr {L}}(({\mathscr {P}}{\mathscr {M}}_{h,{\mathscr {D}}}g))(\cdot ,t))||_{L_q}\le B_{pq} ||{\mathscr {P}}(g- {\mathscr {M}}_{h,{\mathscr {D}}}g)(\cdot ,t))||_{L_p} \end{aligned}$$

where $B_{pq}$ denotes the norm of the bounded mapping $(-\Delta )^{-1/2}:L_p\rightarrow L_q$ [17, Theorem V.1]. From [15, (6.14)], [16, (2.68)]) we see that

$$\begin{aligned} ||{\mathscr {P}}(g- {\mathscr {M}}_{h,{\mathscr {D}}}g)(\cdot ,t))||_{L_p}\le ||g- {\mathscr {M}}_{h,{\mathscr {D}}}g||_{L_p} \end{aligned}$$

(3.9)

for any $t>0$ and $p\ge 1$. In addition, the saturation error converges to zero with the order ${\mathscr {O}}(h^N)$. In [15, Paragraph 6.2.1] the inequality

$$\begin{aligned} \Big | \frac{1}{(4\pi \kappa t)^{3/2}}&\int \limits _{{{\mathbb {R}}}^3} \textrm{e}^{-\frac{|{{\textbf{x}}}-{{\textbf{y}}}|^2}{4\kappa t}} \, \partial ^{\varvec{\alpha }} g({{\textbf{y}}}) \, \rho _{\varvec{\alpha }} \big (\frac{{{\textbf{y}}}}{h},\eta ,{\mathscr {D}}\big ) d{{\textbf{y}}}\Big | \\&\le \frac{c_{\varvec{\alpha }} h^{N-{|\varvec{\alpha }|}}}{(4\pi \kappa t)^{(N-{|\varvec{\alpha }|)/2}}}\sum _{\varvec{\nu } \in {{{\mathbb {Z}}}}^3\setminus \{0\}} |\partial ^{\varvec{\alpha }}{\mathscr {F}}\eta (\sqrt{{\mathscr {D}}} {\varvec{\nu }})| |\varvec{\nu }|^{|\varvec{\alpha }|-N} \end{aligned}$$

is proved with a constant $c_{\varvec{\alpha }}{} { dependingon}g{ and}t$. This shows that

$$\begin{aligned} \Big \Vert \frac{1}{(4\pi \kappa t)^{3/2}} \sum _{ |\varvec{\alpha }| = 0} ^{N-1} \frac{(\sqrt{{\mathscr {D}}}h)^{|\varvec{\alpha }|}}{\varvec{\alpha } !(2 \pi i)^{|{\varvec{\alpha }}|}} \int \limits _{{{\mathbb {R}}}^3} \textrm{e}^{-\frac{|{{\textbf{x}}}-{{\textbf{y}}}|^2}{4\kappa t}} \, \partial ^{\varvec{\alpha }} g({{\textbf{y}}}) \, \rho _{\varvec{\alpha }} \big (\frac{{{\textbf{y}}}}{h},\eta ,{\mathscr {D}}\big ) d{{\textbf{y}}}\Big \Vert _{L_p} \le c\, h^N \, . \end{aligned}$$

Hence, by Theorem 1.1 the assertion follows. $\square$

We assume the basis function (2.2). Keeping in mind (2.5) and (2.7) we have, for $b>0$

$$\begin{aligned} \int \limits _{{{\mathbb {R}}}^3} \textrm{e}^{- \frac{|{{\textbf{z}}}-{{\textbf{w}}}|^2}{b}} \eta \left( {{\textbf{w}}}\right) d{{\textbf{w}}}= & {} \prod _{j=1}^3 A\left( \frac{d}{dz_j}\right) \int _{{\mathbb {R}}}\textrm{e}^{-w^2} \textrm{e}^{-(z_j-w^2)/b}dw\\= & {} \left( \frac{\pi b}{1+b}\right) ^{3/2} \prod _{j=1}^3 A\left( \frac{d}{dz_j}\right) \textrm{e}^{-z_j^2/(1+b)}= \left( \frac{\pi b}{1+b}\right) ^{3/2} \prod _{j=1}^3 {\mathscr {Q}}_M(z_j,b)\textrm{e}^{-z_j^2/(1+b)}. \end{aligned}$$

Substituting in (3.6) we obtain, for $k=1,2,3$,

$$\begin{aligned} ({\mathscr {N}}^{(M)}_{h,{\mathscr {D}}} g)_k({{\textbf{x}}},t)= & {} -c_{\gamma ,\lambda +2\mu }{ {\kappa \, t}\, h^3}{}\sum _{{{\textbf{m}}}\in {{\mathbb {Z}}}^3} g(h{{\textbf{m}}}) \nonumber \\{} & {} \times \int \limits _0^\infty \frac{\partial }{\partial x_k} \prod _{j=1}^3{\mathscr {Q}}_M\left( \frac{x_j-hm_j}{h\sqrt{{\mathscr {D}}}},(1+\tau ) \frac{4\kappa t}{h^2{\mathscr {D}}}\right) \frac{\textrm{e}^{-\frac{|{{\textbf{x}}}-h{{\textbf{m}}}|^2}{h^2{\mathscr {D}}+(1+\tau ) 4\kappa t}} }{(h^2 {\mathscr {D}}+(1+\tau ) 4\kappa t)^{3/2}}d\tau \nonumber \\= & {} c_{\gamma ,\lambda +2\mu }\frac{{\kappa \, t}}{h\,{\mathscr {D}}^{2}}\sum _{{{\textbf{m}}}\in {{\mathbb {Z}}}^3} g(h{{\textbf{m}}})\int \limits _0^\infty \prod _{\begin{array}{c} j=1 j\ne k \end{array}}^3{\mathscr {Q}}_M\left( \frac{x_j-hm_j}{h\sqrt{{\mathscr {D}}}},(1+\tau ) \frac{4\kappa t}{h^2{\mathscr {D}}}\right) \nonumber \\{} & {} \times {\mathscr {R}}_M\left( \frac{x_k-hm_k}{h\sqrt{{\mathscr {D}}}},(1+\tau ) \frac{4\kappa t}{h^2{\mathscr {D}}}\right) \textrm{e}^{-\frac{|{{\textbf{x}}}-h{{\textbf{m}}}|^2}{h^2{\mathscr {D}}+(1+\tau ) 4\kappa t}} \frac{ (h\sqrt{{\mathscr {D}}})^3}{(h^2 {\mathscr {D}}+(1+\tau ) 4\kappa t)^{3/2}}d\tau , \end{aligned}$$

(3.10)

where

$$\begin{aligned} {\mathscr {R}}_M(x,\Lambda ) =\frac{2x}{1+\Lambda }{\mathscr {Q}}_M(x,\Lambda )+2\sqrt{1+\Lambda }{\mathscr {A}}_M(x,\Lambda ) \end{aligned}$$

with

$$\begin{aligned} {\mathscr {A}}_M(x,\Lambda )=\frac{2}{\sqrt{\pi }}\sum _{s=1}^{M-1}\frac{1}{(1+\Lambda )^{s-1/2}}\frac{(-1)^{s-1}}{(s-1)! 4^s}H_{2s-1}\left( \frac{x}{\sqrt{1+\Lambda }}\right) . \end{aligned}$$

For example, for $M=1$ we get the following formula suitable for fast computation

$$\begin{aligned}{} & {} ({\mathscr {N}}^{(1)}_{h,{\mathscr {D}}} g)_k({{\textbf{x}}},t)\\{} & {} \quad = \frac{2}{\pi ^{3/2}}c_{\gamma ,\lambda +2\mu }\, {\kappa \, t}\, h^4\sqrt{{\mathscr {D}}}\sum _{{{\textbf{m}}}\in {{\mathbb {Z}}}^3} g(h{{\textbf{m}}}) \frac{x_k-hm_k}{h\sqrt{{\mathscr {D}}}} \int \limits _0^\infty \frac{\textrm{e}^{-\frac{|{{\textbf{x}}}-h{{\textbf{m}}}|^2}{h^2{\mathscr {D}}+(1+\tau ) 4\kappa t}} }{(h^2 {\mathscr {D}}+(1+\tau ) 4\kappa t)^{5/2}}d\tau . \end{aligned}$$

4 Implementation and Numerical Experiments

In this section we provide numerical experiments for the approximation of ${{\textbf{u}}}^{(2)}$ and $T$ by means of (3.10) and (2.8), respectively.

The quadrature of the one-dimensional integrals which appears in $({\mathscr {N}}^{(M)}_{h,{\mathscr {D}}} g)_k$, $k=1,2,3$, with certain quadrature weights $\omega _p$ and nodes $\tau _p$ leads to the approximation formulas at the point of a uniform grid $\{h{{\textbf{s}}}\}$

$$\begin{aligned}{} & {} u_k^{(2)}(h{{\textbf{s}}},t)\approx ({\mathscr {N}}^{(M)}_{h,{\mathscr {D}}} g)_k(h{{\textbf{s}}},t)\\{} & {} \quad = c_{\gamma ,\lambda +2\mu }\frac{{\kappa \, t}}{h{{\mathscr {D}}^2}}\sum _{{{\textbf{m}}}\in {{\mathbb {Z}}}^3} g(h{{\textbf{m}}})\sum _p \omega _p \prod _{\begin{array}{c} j=1 j\ne k \end{array}}^3{\mathscr {Q}}_M\left( \frac{s_j-m_j}{\sqrt{{\mathscr {D}}}},(1+\tau _p) \frac{4\kappa t}{h^2{\mathscr {D}}}\right) \\{} & {} \qquad \times {\mathscr {R}}_M\left( \frac{s_k-m_k}{\sqrt{{\mathscr {D}}}},(1+\tau _p) \frac{4\kappa t}{h^2{\mathscr {D}}}\right) \textrm{e}^{-\frac{|{{\textbf{s}}}-{{\textbf{m}}}|^2}{{\mathscr {D}}+(1+\tau _p) 4\kappa t/h^2}} \frac{ (h\sqrt{{\mathscr {D}}})^3}{(h^2 {\mathscr {D}}+(1+\tau _p) 4\kappa t)^{3/2}}. \end{aligned}$$

Table 1 Exact values $u_1^{(2)}$ in (4.2) at some grid points ${{\textbf{x}}}=(x,x,x)$ and $t=1$, approximated values using ${\mathscr {N}}_{0.025,2}^{(3)}$ and relative errors

Table 2 Absolute error and rate of convergence for $u_1^{(2)}$ in (4.2) at ${{\textbf{x}}}=(1,0,0)$ and $t=1$ using ${\mathscr {N}}_{h,2}^{(M)}$

Table 3 Absolute error and rate of convergence for $u_1^{(2)}$ in (4.2) at ${{\textbf{x}}}=(0.8,0.8,0.8)$ and $t=2$ using ${\mathscr {N}}_{h,2}^{(M)}$

Table 4 Exact values $T$ in (4.1) at some grid points ${{\textbf{x}}}=(x,x,x)$ and $t=1$, approximated values using $T^{(3)}_{0.025,2}$ and relative errors

Table 5 Absolute error and rate of convergence for $T$ in (4.1) at ${{\textbf{x}}}=(1,0,0)$ and $t=1$ using $T^{(M)}_{h,2}$

Table 6 Absolute error and rate of convergence for $T$ in (4.1) at ${{\textbf{x}}}=(0.8,0.8,0.8)$ and $t=2$ using $T^{(M)}_{h,2}$

The approximation formulas $({\mathscr {N}}^{(M)}_{h,{\mathscr {D}}} g)_k,k=1,2,3$ are very efficient if $g$ has a separated representation (2.10). Then an approximate value of $u_k^{(2)}(h{{\textbf{s}}},t)$ can be approximated using only one-dimensional operations as follows

$$\begin{aligned} u_k^{(2)}(h{{\textbf{s}}},t)\approx & {} ({\mathscr {N}}^{(M)}_{h,{\mathscr {D}}} g)_k(h{{\textbf{s}}},t)\\\approx & {} \frac{\gamma }{\lambda +2\mu } \frac{{\kappa }\,t}{h{{\mathscr {D}}^2}}\sum _{\ell =1}^L \sum _p \omega _p R_k^{(\ell )}(s_k,(1+\tau _p) \frac{4\kappa t}{h^2{\mathscr {D}}}) \prod _{\begin{array}{c} j=1 j\ne k \end{array}}^3T_j^{(\ell )}\left( s_j,(1+\tau _p) \frac{4\kappa t}{h^2{\mathscr {D}}}\right) \end{aligned}$$

with the one-dimensional convolutions

$$\begin{aligned} T_r^{(\ell )}(s,\Lambda )= & {} \Lambda ^{-1}\sum _{m\in {{\mathbb {R}}}} g_r^{(\ell )}(hm) {\mathscr {Q}}_M\left( \frac{s-m}{\sqrt{{\mathscr {D}}}},\Lambda \right) \textrm{e}^{-\frac{(s-m)^2}{{\Lambda {\mathscr {D}}}}},\\ R_r^{(\ell )}(s,\Lambda )= & {} \Lambda ^{-1}\sum _{m\in {{\mathbb {R}}}} g_r^{(\ell )}(hm) {\mathscr {R}}_M\left( \frac{s-m}{\sqrt{{\mathscr {D}}}},\Lambda \right) \textrm{e}^{-\frac{(s-m)^2}{{\Lambda {\mathscr {D}}}}} . \end{aligned}$$

We provide results of some experiments which show the accuracy and the convergence order of the method. We compute the solution of (1.7),(1.2),(1.3) with $g({{\textbf{x}}})=\textrm{e}^{-|{{\textbf{x}}}|^2}$. The exact solution of (1.2)–(1.3) is given by

$$\begin{aligned} T({{\textbf{x}}},t)={\mathscr {P}}(\textrm{e}^{-|\cdot |^2})({{\textbf{x}}},t)=\frac{\textrm{e}^{-|{{\textbf{x}}}|^2/(1+4 \kappa t)}}{(1+4 \kappa t)^{3/2}} \end{aligned}$$

(4.1)

and, by using (3.2) and

$$\begin{aligned} {\mathscr {L}}(\textrm{e}^{-|\cdot |^2})({{\textbf{x}}})=\frac{\sqrt{\pi }}{4|{{\textbf{x}}}|}\textrm{erf}\,(|{{\textbf{x}}}|), \end{aligned}$$

we get

$$\begin{aligned} u_k^{(2)}({{\textbf{x}}},t)=-c_{\gamma ,\lambda +2\mu }\frac{\sqrt{\pi }}{4}\frac{\partial }{\partial x_k} \frac{\textrm{erf}\,\left( \frac{|{{\textbf{x}}}|}{\sqrt{1+4\kappa t}}\right) }{|{{\textbf{x}}}|} \end{aligned}$$

(4.2)

$$\begin{aligned}= \frac{c_{\gamma ,\lambda +2\mu }}{4}\frac{x_k}{|{{\textbf{x}}}|^2} \left( \sqrt{\pi } \frac{\textrm{erf}\,\left( \frac{|{{\textbf{x}}}|}{\sqrt{1+4 \kappa t}}\right) }{|{{\textbf{x}}}|}-\frac{\textrm{e}^{-|{{\textbf{x}}}|^2/(1+4 \kappa t)}}{\sqrt{1+4\kappa t}} \right) ,\quad k=1,2,3. \end{aligned}$$

We assume $\kappa =1$ and the parameters $\gamma ,\lambda ,\mu$ such that $c_{\gamma ,\lambda +2\mu }=1$.

Following [18] the one-dimensional integrals in (3.10) are transformed to integrals over ${{\mathbb {R}}}$ with integrands decaying doubly exponentially by making the substitutions

$$\begin{aligned} t=\textrm{e}^\xi ,\quad \xi =\alpha (\sigma +\textrm{e}^\sigma ),\quad \sigma =\beta (u-\textrm{e}^{-u}) \end{aligned}$$

(4.3)

with certain positive constants $\alpha ,\beta$, and the computation is based on the classical trapezoidal rule. Then the tensor product structure of the integrands allows the efficient computation of ${\mathscr {N}}_{h,{\mathscr {D}}}^{(M)}g$.

In Table 1 we compare the exact values $u_1^{(2)}$ in (4.2) and the approximates values $({\mathscr {N}}_{0.025,2}^{(3)}(\textrm{e}^{-|\cdot |^2}))_1$ at some grid points ${{\textbf{x}}}=(x,x,x)$ and $t=1$. In Tables 2 and 3 we report on the absolute errors and approximate rates for the computation of $u_1^{(2)}$ at ${{\textbf{x}}}=(1,0,0)$, $t=1$ and ${{\textbf{x}}}=(0.8,0.8,0.8)$, $t=2$, respectively. The approximate values are computed by the formulas $({\mathscr {N}}_{h,2}^{(M)}(\textrm{e}^{-|\cdot |^2}))_1$ for $M=1,2,3,4$ and uniform grids size $h=0.1\times 2^{-s}$, $s=0, \ldots ,4$. The convergence rate is calculated as

$$\begin{aligned} (\log |u_1^{(2)}-{\mathscr {N}}_{2h,2}^{(M)}(\textrm{e}^{-|\cdot |^2})_1|-\log |u_1^{(2)}-{\mathscr {N}}_{h,2}^{(M)}(\textrm{e}^{-|\cdot |^2})_1|)/\log 2. \end{aligned}$$

We have chosen ${\alpha }=6$, $\beta =5$ in the transformation (4.3) and $\tau =0.003$ with $600$ terms in the trapezoidal rule. The numerical results confirm the $h^{2M}$ convergence of the approximating formula when $M=1,2,3,4$. For small $h$, the 8th-order formula has reached the machine precision.

In the next tables we report on numerical experiments for the approximation of $T$ in (4.1) by means of (2.8).

In Table 4 we compare the values of the exact solution and the approximate solution at some points. The approximations in Table 4 have been computed on a uniform grid with step size $h = 0.025$ and $N = 6$.

In Tables 5 and 6 we show that formula (2.8) approximates the exact solution with the predicted approximate orders $h^{2M}$ with $M=1,2,3,4$. For small $h$, the 6th-order and 8th-order formulas have reached the machine precision.