High-Performance Parallel Solver for Integral Equations of Electromagnetics Based on Galerkin Method

Abstract

A new parallel solver for the volumetric integral equations (IE) of electrodynamics is presented. The solver is based on the Galerkin method, which ensures convergent numerical solution. The main features include: (i) memory usage eight times lower compared with analogous IE-based algorithms, without additional restrictions on the background media; (ii) accurate and stable method to compute matrix coefficients corresponding to the IE; and (iii) high degree of parallelism. The solver’s computational efficiency is demonstrated on a problem of magnetotelluric sounding of media with large conductivity contrast, revealing good agreement with results obtained using the second-order finite-element method. Due to the effective approach to parallelization and distributed data storage, the program exhibits perfect scalability on different hardware platforms.

Appendices

Appendix A: Green’s Tensor

Following the notations from (Dmitriev et al. 2002), the electrical $\widehat{G}^E$ and magnetic $\widehat{G}^H$ tensors of layered media can be written as

$$\begin{aligned} \begin{aligned} \widehat{G}^E&={\widehat{G}+\mathop {\mathrm {grad}}\nolimits \left( \frac{\mu _0}{k^2}\mathop {\mathrm {div}}\nolimits \frac{\widehat{G}}{\mu _0}\right) },\\ \widehat{G}^H&=\frac{1}{\mathrm {i}\omega \mu _0}\mathop {\mathrm {curl}}\nolimits \widehat{G},\\ k^2&=\mathrm {i}\omega \mu _0{\sigma }_{\text {b}}, \end{aligned} \end{aligned}$$

(21)

where

$$\begin{aligned} \widehat{G}(M,M_0)=\begin{pmatrix} G_1(M,M_0) &{} 0 &{} 0\\ 0 &{} G_1(M,M_0) &{} 0\\ \frac{\partial {g(M,M_0)}}{\partial {x}} &{} \frac{\partial {g(M,M_0)}}{\partial {y}} &{} G_2(M,M_0) \end{pmatrix} \end{aligned}$$

(22)

and

$$\begin{aligned}&G_1(M,M_0)=\frac{\mathrm {i}\omega \mu _0}{4\pi }\int \limits _0^\infty J_0(\lambda \rho )U_1(\lambda ,z,z_0)\lambda \mathrm{d}\lambda ,\nonumber \\&G_2(M,M_0)=\frac{\mathrm {i}\omega \mu _0}{4\pi }\int \limits _0^\infty J_0(\lambda \rho )U_\sigma (\lambda ,z,z_0)\lambda \mathrm{d}\lambda ,\nonumber \\&g(M,M_0)=-\frac{\mathrm {i}\omega \mu _0}{4\pi }\int \limits _0^\infty J_0(\lambda \rho )\left( \frac{\partial { }}{\partial {z_0}}U_\sigma (\lambda ,z,z_0)+\frac{\partial { }}{\partial {z}} U_1(\lambda ,z,z_0)\right) \frac{\mathrm{d}\lambda }{\lambda },\nonumber \\&M=M(x,y,z)\quad M_0=M_0(x_0,y_0,z_0),\quad \rho =\sqrt{(x-x_0)^2+(y-y_0)^2}. \end{aligned}$$

(23)

Here, $J_0$ is a zero-order Bessel function of the first kind, and functions $U_\gamma (\lambda ,z,z_0)$, $\gamma =1,\sigma$ are the fundamental functions of the layered media (Appendix B).

From Eqs. (21) and (22), one gets

$$\begin{aligned} \widehat{G}^E=\left\{ \begin{array}{lll} G_1+\frac{1}{k^2}\frac{\partial ^2{ }}{\partial {x^2}}\left( G_1+\frac{\partial {g}}{\partial {z}}\right) &{} \frac{1}{k^2}\frac{\partial ^2{ }}{\partial {x}\partial {y}}\left( G_1+\frac{\partial {g}}{\partial {z}}\right) &{}\frac{1}{k^2}\frac{\partial ^2{G_2}}{\partial {x}\partial {z}}\\ \frac{1}{k^2}\frac{\partial ^2{ }}{\partial {x}\partial {y}}\left( G_1+\frac{\partial {g}}{\partial {z}}\right) &{} G_1+\frac{1}{k^2}\frac{\partial ^2{ }}{\partial {y^2}}\left( G_1+\frac{\partial {g}}{\partial {z}}\right) &{} \frac{1}{k^2}\frac{\partial ^2{G_2}}{\partial {y}\partial {z}}\\ \frac{\partial {g}}{\partial {x}}+\frac{1}{k^2}\frac{\partial ^2{ }}{\partial {x}\partial {z}}\left( G_1+ \frac{\partial {g}}{\partial {z}}\right) &{}\frac{\partial {g}}{\partial {y}}+\frac{1}{k^2}\frac{\partial ^2{ }}{\partial {y}\partial {z}}\left( G_1+\frac{\partial {g}}{\partial {z}} \right) &{} G_2+\frac{1}{k^2}\frac{\partial ^2{G_2}}{\partial {z^2}} \end{array}\right\} . \end{aligned}$$

(24)

Note that $G^E_{xz}(M,M_0)=-G^E_{zx}(M_0,M)$ and $G^E_{yz}(M,M_0)=-G^E_{zy}(M_0,M)$, according to Lorentz reciprocity.

Let $\mathrm {\Omega }_n=S_n \times [z_n^1,z_n^2]$, $\mathrm {\Omega }_m=S_m \times [z_m^1,z_m^2]$, where $S_n, S_m$ are horizontal rectangular domains, and let ${\sigma }_{\text {b}}={\sigma }_{\text {b}}^n$ inside $[z_n^1, z_n^2]$, ${\sigma }_{\text {b}}={\sigma }_{\text {b}}^m$ inside $[z_m^1, z_m^2]$; That is, the subdomains do not intersect the boundaries of the layers. Taking into account Eqs. (12), (21), and (24), one can see that $\widehat{B}_n^m$ is expressed in terms of double volumetric integrals with weak integrable singularity, so the order of integration can be changed. Then, using Eqs.  (23) and (24), one obtains

$$\begin{aligned} \widehat{B}_n^m=\frac{1}{4\pi }\left\{ \begin{aligned}&I_1^{n,m}+I_{xx}^{n,m}&I_{xy}^{n,m}&I_x^{n,m}\\&I_{xy}^{n,m}&I_1^{n,m}+I_{yy}^{n,m}&I_y^{n,m}\\&- I_x^{m,n}&- I_y^{m,n}&I_2^{n,m}+I_{zz}^{n,m}\\ \end{aligned}\right\} , \end{aligned}$$

(25)

where

$$\begin{aligned} I_1^{n,m}&=\int \limits _{S_n}\int \limits _{S_m}\left( \int \limits _0^{\infty }J_0(\lambda \rho )V_1^{n,m}(\lambda )\lambda \mathrm{d}\lambda \right) \mathrm{d}x_0\mathrm{d}y_0\mathrm{d}x\mathrm{d}y,\nonumber \\ I_2^{n,m}&=\int \limits _{S_n}\int \limits _{S_m}\left( \int \limits _0^{\infty }J_0(\lambda \rho )V_2^{n,m}(\lambda )\lambda \mathrm{d}\lambda \right) \mathrm{d}x_0\mathrm{d}y_0\mathrm{d}x\mathrm{d}y,\nonumber \\ I_{\alpha \beta }^{n,m}&=\int \limits _{S_n}\frac{\partial ^2{ }}{\partial {\alpha }\partial {\beta }}\left\{ \int \limits _{S_m}\left( \int \limits _0^{\infty }J_0(\lambda \rho )\left[ V_1^{n,m}(\lambda )+V_3^{n,m}(\lambda )\right] \frac{\mathrm{d}\lambda }{\lambda }\right) \mathrm{d}x_0\mathrm{d}y_0\right\} \mathrm{d}x\mathrm{d}y,\nonumber \\ I_{\alpha }^{n,m}&=\int \limits _{S_n}\frac{\partial {}}{\partial {\alpha }}\left\{ \int \limits _{S_m}\left( \int \limits _0^{\infty }J_0(\lambda \rho )V_4^{n,m}(\lambda )\lambda \mathrm{d}\lambda \right) \mathrm{d}x_0\mathrm{d}y_0\right\} \mathrm{d}x\mathrm{d}y,\nonumber \\ I_{zz}^{n,m}&=\int \limits _{S_n}\int \limits _{S_m}\left( \int \limits _0^{\infty }J_0(\lambda \rho )V_5^{n,m}(\lambda )\lambda \mathrm{d}\lambda \right) \mathrm{d}x_0\mathrm{d}y_0\mathrm{d}x\mathrm{d}y. \end{aligned}$$

(26)

Here $\alpha =x,y$, $\beta =x,y$, and

$$\begin{aligned}&V_1^{n,m}(\lambda )=\mathrm {i}\omega \mu _0\int \limits _{z_n^1}^{z_n^2}\int \limits _{z_m^1}^{z_m^2}U_1(\lambda ,z,z_*)\mathrm{d}z_*\mathrm{d}z,\nonumber \\&V_2^{n,m}(\lambda )=\mathrm {i}\omega \mu _0\int \limits _{z_n^1}^{z_n^2}\int \limits _{z_m^1}^{z_m^2}U_{\sigma }(\lambda ,z,z_*)\mathrm{d}z_*\mathrm{d}z,\nonumber \\&V_3^{n,m}(\lambda )=-\frac{1}{{\sigma }_{\text {b}}^n}\int \limits _{z_n^1}^{z_n^2}\frac{\partial { }}{\partial {z}} \left( \int \limits _{z_m^1}^{z_m^2}\left[ \frac{\partial { }}{\partial {z_*}}U_{\sigma }(\lambda ,z,z_*)\right] \mathrm{d}z_*\right) \mathrm{d}z,\nonumber \\&V_4^{n,m}(\lambda )=\frac{1}{{\sigma }_{\text {b}}^n}\int \limits _{z_n^1}^{z_n^2}\frac{\partial { }}{\partial {z}}\left( \int \limits _{z_m^1}^{z_m^2}U_{\sigma }(\lambda ,z,z_*)\mathrm{d}z_*\right) \mathrm{d}z,\nonumber \\&V_5^{n,m}(\lambda )=\frac{1}{{\sigma }_{\text {b}}^n}\int \limits _{z_n^1}^{z_n^2}\frac{\partial ^2{ }}{\partial {z^2}}\left( \int \limits _{z_m^1}^{z_m^2}U_{\sigma }(\lambda ,z,z_*)\mathrm{d}z_*\right) \mathrm{d}z. \end{aligned}$$

(27)

Therefore, to obtain the coefficients of $\hat{B}_n^m$, one needs computational methods to find “horizontal” integrals (26) and “vertical” integrals (27). These methods are presented in the following appendices.

Appendix B: Vertical Integration

The integrals in formulas (27) are expressed in terms of the so-called fundamental function of layered media (Dmitriev et al. 2002). Consider media with $N_\mathrm{lay}-1$ homogeneous layers with complex conductivities ${\sigma }_n$, $n=1\ldots N_\mathrm{lay}-1$, the upper halfspace (air, the zeroth layer) with complex conductivity ${\sigma }_0$, and the lower halfspace (the $N_\mathrm{lay}$-th layer) with conductivity ${\sigma }_{N_\mathrm{lay}}$. Note that, in EM sounding problems, typically $\mathop {\mathbf {Re}}{{\sigma }_0} \le 10^{-9}$.

The function $U_\gamma (z,z_*,\lambda )$ is defined as a unique solution of the problem

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\partial {^2 }}{\partial {z^2}}U_\gamma (z,z_*,\lambda ) -\eta _0^2U_\gamma (z,z_*,\lambda )=0,z<d_1, z \not =z_*,\\&\frac{\partial {^2 }}{\partial {z^2}}U_\gamma (z,z_*,\lambda ) -\eta _n^2U_\gamma (z,z_*,\lambda )=0, d_n<z<d_{n+1}, z \not =z_*,\\&\frac{\partial {^2 }}{\partial {z^2}}U_\gamma (z,z_*,\lambda ) -\eta _{N_\mathrm{lay}}^2U_\gamma (z,z_*,\lambda )=0, z >d_{N_\mathrm{lay}}, z \not =z_*,\\&\left[ U_\gamma (z,z_*,\lambda )\right] _{z=d_n}=0,\\&\left[ \frac{1}{\gamma }\frac{\partial { }}{\partial {z}}U_\gamma (z,z_*,\lambda )\right] _{z=d_n}=0,\\&\left[ U_\gamma (z,z_*,\lambda )\right] _{z=z_*}=0,\\&\left[ \frac{\partial { }}{\partial {z}}U_\gamma (z,z_*,\lambda )\right] _{z=z_*}=-2,\\&\left| U_\gamma (z,z_*,\lambda )\right| \rightarrow 0 \text { as } z \rightarrow \pm \infty ,\\&\eta _m^2=\lambda ^2-k^2_m,\quad k^2_m=\mathrm {i}\omega \mu _0{\sigma }_m,\\&m=0\ldots N_\mathrm{lay},\quad n=1\ldots N_\mathrm{lay}-1,\quad 0< \lambda < \infty . \end{aligned} \right. \end{aligned}$$

(28)

The following procedure is performed to obtain an explicit expression for $U_\gamma$ that allows analytical integration. Let $l_0=0$, $l_{N_\mathrm{lay}}=0$, $l_n=d_{n+1}-d_n$, $n=1\ldots N_\mathrm{lay}-1$. Define $p_{m}^\gamma$, $q_{m}^\gamma$, ${m=0\ldots N_\mathrm{lay}}$ by the recurrent expressions

$$\begin{aligned} \begin{aligned}&p_0^\gamma =0;&q_{N_\mathrm{lay}}^\gamma =0;\\&p_{1}^\gamma =\frac{1-\alpha _0^\gamma \frac{\eta _0}{\eta _{1}}}{1+\alpha _0^\gamma \frac{\eta _0}{\eta _{1}}};\quad&q_{N_\mathrm{lay}-1}^\gamma =\frac{1-\beta _{N_\mathrm{lay}}^\gamma \frac{\eta _{N_\mathrm{lay}}}{\eta _{N_\mathrm{lay}-1}}}{1+\beta _{N_\mathrm{lay}}^\gamma \frac{\eta _{N_\mathrm{lay}}}{\eta _{N_\mathrm{lay}-1}}};\\&p_{m+1}^\gamma =\frac{1+\alpha _m^\gamma \frac{\eta _m}{\eta _{m+1}}\frac{p_{m}^\gamma \mathrm {e}^{-2\eta _ml_m}-1}{p_{m}^\gamma \mathrm {e}^{-2\eta _ml_m}+1}}{1-\alpha _m^\gamma \frac{\eta _m}{\eta _{m+1}}\frac{p_{m}^\gamma \mathrm {e}^{-2\eta _ml_m}-1}{p_{m}^\gamma \mathrm {e}^{-2\eta _ml_m}+1}}, \quad&q_{m-1}^\gamma =\frac{1+\beta _m^\gamma \frac{\eta _m}{\eta _{m-1}}\frac{q_{m}^\gamma \mathrm {e}^{-2\eta _ml_m}-1}{q_{m}^\gamma \mathrm {e}^{-2\eta _ml_m}+1}}{1-\beta _m^\gamma \frac{\eta _m}{\eta _{m-1}}\frac{q_{m}^\gamma \mathrm {e}^{-2\eta _ml_m}-1}{q_{m}^\gamma \mathrm {e}^{-2\eta _ml_m}+1}},\\& m\not =N_\mathrm{lay}; \qquad& m\not =0;\\&\alpha _m^\gamma =\left\{ \begin{aligned} 1,\quad \gamma =1;\\ \frac{{\sigma }_{m+1}}{{\sigma }_{m}},\quad \gamma ={\sigma }; \end{aligned} \right.&\beta _m^\gamma =\left\{ \begin{aligned} 1,\quad \gamma =1;\\ \frac{{\sigma }_{m-1}}{{\sigma }_{m}},\quad \gamma ={\sigma }. \end{aligned} \right. \end{aligned} \end{aligned}$$

(29)

Let $d_0=d_1$, $d_{N_{lay+1}}=d_{N_\mathrm{lay}}$ and let points $z_r$, $z_s$ belong to layers r and s, respectively, $0 \le r, s \le N_\mathrm{lay}$. Then, using expressions (29), one gets

$$\begin{aligned} U_\gamma (z_r,z_s,\lambda )= {\left\{ \begin{array}{ll} &{}A_{r,s}^\gamma \left( p_r^\gamma \mathrm {e}^{2\eta _rd_r}\mathrm {e}^{-\eta _rz_r}+\mathrm {e}^{\eta _r z_r}\right) \left( \mathrm {e}^{-\eta _sz_s}+q_s^\gamma \mathrm {e}^{-2\eta _sd_{s+1}}\mathrm {e}^{\eta _s z_s}\right) \\ &{}\text {for}\quad z_r\le z_s;\\ &{}A_{s,r}^\gamma \left( p_s^\gamma \mathrm {e}^{2\eta _sd_s}\mathrm {e}^{-\eta _sz_s}+\mathrm {e}^{\eta _s z_s}\right) \left( \mathrm {e}^{-\eta _rz_r}+q_r^\gamma \mathrm {e}^{-2\eta _rd_{r+1}}\mathrm {e}^{\eta _r z_r}\right) \\ &{}\text {for}\quad z_r>z_s, \end{array}\right. } \end{aligned}$$

(30)

where

$$\begin{aligned} \begin{aligned} A_{r,s}^\gamma&=Q_r^\gamma \times Q_{r+1}^\gamma \times \cdots \times Q_{s-1}^\gamma A_{s,s}^\gamma ,\quad \text {for} \quad r < s,\\ Q_m^\gamma&=\frac{1+p_{m+1}^\gamma }{1+p_m^\gamma \mathrm {e}^{-2\eta _m l_m }}\mathrm {e}^{(\eta _{m+1}-\eta _m)d_{m+1}},\quad \text {for}\quad m=1\ldots N_\mathrm{lay}-1,\\ A_{n,n}^\gamma&=\frac{1}{\eta _n\left( 1-p_n^\gamma q_n^\gamma \mathrm {e}^{-2\eta _nl_n}\right) },\quad \text { for} \quad r=s=n, n=0\ldots N_\mathrm{lay},\\ A_{r,s}^1&=A_{s,r}^1, \quad A_{r,s}^{{\sigma }}=\frac{{\sigma }_r}{{\sigma }_s}A_{s,r}^{{\sigma }},\quad \text { for} \quad r>s. \end{aligned} \end{aligned}$$

(31)

To check Eq. (30), one can explicitly substitute Eq. (30) into Eq. (28), taking into account expressions (29) and (31).

In view of Eq. (30), one can see that integrals in formulas (27) (i.e., the integrals over $U_\gamma (z,z_*,\lambda )$ and its partial derivatives) can be obtained analytically with respect to z, $z_*$ over any domains that do not intersect the layer boundaries. However, rounding errors arising in addition and multiplication of very small or large quantities make the formula (30) impractical for $\lambda \gg 1$. Instead, the following formula is used:

$$\begin{aligned} U_\gamma (z_r,z_s,\lambda )= {\left\{ \begin{array}{ll} &{}A^\gamma _{r,s}\left( p_r^\gamma \mathrm {e}^{-(\eta _rz_r+\eta _sz_s-2\eta _rd_r)}+q_s^\gamma \mathrm {e}^{-(2\eta _s d_{s+1}-(\eta _rz_r+\eta _sz_s))}\right. +\\ &{}\qquad \left. \mathrm {e}^{-(\eta _sz_s-\eta _rz_r)}+p_r^\gamma q_s^\gamma \mathrm {e}^{-(2(\eta _sd_{s+1}-\eta _rd_r)-(\eta _sz_s-\eta _rz_r))}\right) \\ &{}\text {for}\quad z_r \le z_s;\\ &{}A^\gamma _{s,r}\left( p_s^\gamma \mathrm {e}^{-(\eta _sz_s+\eta _rz_r-2\eta _sd_s)}+q_r^\gamma \mathrm {e}^{-(2\eta _r d_{r+1}-(\eta _sz_s+\eta _rz_r))}\right. +\\ &{}\qquad \left. \mathrm {e}^{-(\eta _rz_r-\eta _sz_s)}+p_s^\gamma q_r^\gamma \mathrm {e}^{-(2(\eta _rd_{r+1}-\eta _sd_s)-(\eta _rz_r-\eta _sz_s))}\right) \\ &{}\text {for}\quad z_r>z_s. \end{array}\right. } \end{aligned}$$

(32)

Formula (32) overcomes the aforementioned problem, since the real parts of all the exponents’ powers are negative. The consequent calculations provide accurate and robust results for any $0< \lambda <\infty$.

Consider $N_z$ subdomains in the discretization in vertical direction. To obtain the matrix $\widehat{B}_n^m$ for the system (15), one needs to compute $O\left( N_z^2\right)$ complex exponents in Eq. (32). An algorithm requiring only $O\left( N_z\right)$ complex exponent calculations has been developed to speed up this integration procedure.

Let $z_0<z_1< \cdots < z_{N_z}$. Suppose that the intervals $[z_l,z_{l+1}]$, $l=0\ldots N_z-1$ do not intersect the layers’ boundaries. For $i,j=0\ldots N_z-1$, $0 \le \alpha +\beta \le 2$, one needs to calculate $W_{i,j}^{\alpha ,\beta }(\gamma )=\int \nolimits _{z_i}^{z_{i+1}}\frac{\partial {^\alpha }}{\partial {z^\alpha }}\left( \int \nolimits _{z_j}^{z_{j+1}}\left[ \frac{\partial {^\beta }}{\partial {z_*^\beta }}\, U_\gamma (z,z_{*},\lambda )\right] \mathrm{d}z_{*}\right) \mathrm{d}z$.

Let $r_l$ be an index of the layer containing $[z_l,z_{l+1}]$, $l=0 \ldots N_{z}-1$. Then, using Eq. (30), one obtains for $z_i < z_j$

$$\begin{aligned} \begin{aligned} W_{i,j}^{\alpha ,\beta }(\gamma )=&\int \limits _{z_i}^{z_{i+1}}\left( \frac{\partial {^\alpha }}{\partial {z^\alpha }} \int \limits _{z_j}^{z_{j+1}}\left[ \frac{\partial {^\beta }}{\partial {z_*^\beta }} U_\gamma (z,z_*,\lambda )\right] \mathrm{d}z_*\right) \mathrm{d}z=\\&H_\gamma ^\alpha (z_i,z_{i+1})\prod \limits _{l={i+1}}^{l={j}}\Theta ^\gamma _l \int \limits _{z_j}^{z_{j+1}}\left( \frac{\partial {^\beta }}{\partial {z_*^\beta }} U\gamma (z_j,z_*,\lambda )\right) \mathrm{d}z_*,\\ \end{aligned} \end{aligned}$$

(33)

where

$$\begin{aligned} \begin{aligned} H_\gamma ^\alpha (z_i,z_{i+1})&=\frac{\int \nolimits _{z_i}^{z_{i+1}}\left( \frac{\partial {^\alpha }}{\partial {z^\alpha }}\left[ p_{r_i}^\gamma \mathrm {e}^{-\eta _{r_i}\left( z+z_{i+1}-2d_{r_i}\right) }+\mathrm {e}^{-\eta _{r_i}\left( z_{i+1}-z\right) }\right] \right) \mathrm{d}z}{p_{r_i}^\gamma \mathrm {e}^{-2\eta _{r_i}\left( z_{i+1}-d_{r_i}\right) }+1},\\ \\ \Theta _l^\gamma&=\varkappa _l^\gamma \frac{ p_{r_{l}}^\gamma \mathrm {e}^{-\eta _{r_l}\left( (z_l+z_{l+1})-2d_{r_l}\right) }+\mathrm {e}^{-\eta _{r_l}( z_{l+1}-z_l)}}{p_{r_{l}}^\gamma \mathrm {e}^{2\eta _{r_{l}}(d_{r_{l}}-z_{l+1})}+1},\\ \varkappa _l^\gamma&=\left\{ \begin{aligned} 1,\quad r_l=r_{l+1}\\ Q_{r_{l+1}}^\gamma ,\quad r_l \not =r_{l+1}\\ \end{aligned} \right\} .\\ \end{aligned} \end{aligned}$$

(34)

All the exponents in Eq. (34) vanish as $\lambda \rightarrow \infty$, so the corresponding computations do not depend on round-off errors due to machine precision. The formulas for $z_i>z_j$ are similar. The integrals $W_{ii}^{\alpha ,\beta }$ are computed analytically using formulas (32). Since $\Theta _l^\gamma$ depends only on $l=1\ldots N_{z}$ and $\gamma ={1,{\sigma }}$, one only needs to calculate $O\left( N_z\right)$ complex exponents using factorization from expressions (33), (34).

Appendix C: Horizontal Integration

The integrals (26) are a particular case of the integral

$$\begin{aligned} \begin{aligned}&I_{\alpha ,\beta }=\int \limits _{S_n}\frac{\partial ^{{\alpha +\beta }}}{\partial {x^\alpha }\partial {y^\beta }}\left\{ \int \limits _{S_m}\left[ \int \limits _0^\infty J_0(\rho \lambda )f(\lambda )\mathrm{d}\lambda \right] \mathrm{d}S_m\right\} \mathrm{d}S_n,\\&\rho =\sqrt{(x-x_0)^2+(y-y_0)^2}, \quad 0 \le \alpha +\beta \le 2, \end{aligned} \end{aligned}$$

(35)

where $f(\lambda )$ is some easily computed function, and $S_n=[x_n,x_n+h_x]\times [y_n,y_n+h_y]$, $S_m=[x_m,x_m+h_x]\times [y_m,y_m+h_y]$ are rectangular domains with similar sizes.

The key feature of the proposed method is transformation of the integrals (35) to one-dimensional convolution integrals. Taking for simplicity $\alpha =\beta =0$, one has

$$\begin{aligned} I_{0,0}=F(R;p,q,\varphi )&=\int \limits _0^\infty K(R\lambda ;p,q,\varphi ) f(\lambda )\frac{\mathrm{d}\lambda }{\lambda ^4},\nonumber \\ K(R\lambda ;p,q,\varphi )&=\lambda ^4\int \limits _{S_m}\int \limits _{S_n} J_0(\rho \lambda ) \mathrm{d}S_m\mathrm{d}S_n\nonumber \\&=\lambda ^4\int \limits _{x_n}^{x_n+h_x}\int \limits _{y_n}^{y_n+h_y}\int \limits _{x_m}^{x_m+h_x}\int \limits _{y_m}^{y_m+h_y}J_0(\rho \lambda )\mathrm{d}x_0\mathrm{d}y_0\mathrm{d}x\mathrm{d}y\nonumber \\&=\int \limits _{0}^{R\lambda p}\int \limits _{0}^{R\lambda q}\int \limits _{R\lambda \left( \cos \varphi -\frac{p}{2}\right) }^{R\lambda \left( \cos \varphi +\frac{p}{2}\right) }\int \limits _{R\lambda \left( \sin \varphi -\frac{q}{2}\right) }^{R\lambda \left( \sin \varphi +\frac{q}{2}\right) }J_0(\tau )d\tilde{x} d\tilde{x}_0 d\tilde{y} d\tilde{y}_0, \end{aligned}$$

(36)

where

$$\begin{aligned}&\tau =\sqrt{(\tilde{x}-\tilde{x}_0)^2+(\tilde{y}-\tilde{y}_0)^2};\\&R=\sqrt{\left( x_n-x_m-\frac{h_x}{2}\right) ^2+\left( y_n-y_m-\frac{h_y}{2}\right) ^2};\\&p=\frac{h_x}{R} \quad q=\frac{h_y}{R} \quad \varphi =\arctan \frac{y_n-y_m-\frac{h_y}{2}}{x_n-x_m-\frac{h_x}{2}}. \end{aligned}$$

Let $\lambda =\mathrm {e}^{-t}$, $R=\mathrm {e}^{s}$

$$\begin{aligned} \begin{aligned} F\left( \mathrm {e}^{s};p,q,\varphi \right) =\mathrm {e}^{3s}\int \limits _{-\infty }^{\infty }\Phi (s-t;p,q,\varphi )f\left( \mathrm {e}^{-t}\right) \mathrm{d}t,\\ \Phi (s-t;p,q,\varphi )=K\left( \mathrm {e}^{s-t};p,q,\varphi \right) \mathrm {e}^{-3(s-t)}. \end{aligned} \end{aligned}$$

(37)

For fixed $p,q,\varphi$ the integral in Eq. (37) is the convolution integral with kernel $\Phi$. Note that, for different values of $\alpha$ and $\beta$, the kernels can be obtained similarly.

The main advantage of using the convolution integrals is that their computation does not require explicit calculation of the kernel $\Phi$. Consider the input function v(t) and the output function u(s) such that

$$\begin{aligned} u(s)=\int \limits _{-\infty }^{+\infty }\Phi (s-t)v(t)\mathrm{d}t. \end{aligned}$$

(38)

For some $N=2M$, l, $0<\xi <0.5l$, $k=-M\ldots M-1$, define

$$\begin{aligned} W_s=(-1)^s\frac{1}{N}\sum \limits _{n=-M}^{M-1}\left\{ \frac{\sum _{m=-M}^{M-1}(-1)^m u(ml-\xi )\mathrm {e}^{-2\mathrm {i}\pi \frac{mn}{N}} }{\sum _{m=-M}^{M-1}(-1)^m v(ml+0.5l)\mathrm {e}^{-2\mathrm {i}\pi \frac{mn}{N} } }\right\} \mathrm {e}^{2\mathrm {i}\pi \frac{sn}{N}}. \end{aligned}$$

(39)

Then, for any g(t), one gets

$$\begin{aligned} \int \limits _{-\infty }^{+\infty }\Phi (s-t)g(t)\mathrm{d}t \approx \sum \limits _{s=N_1}^{N_2} W_s g(sl+\xi ), \quad -M\le N_1 < N_2 \le M-1. \end{aligned}$$

(40)

The tradeoff between the accuracy of approximation (40) and the computational time is achieved by the particular selection of M, $N_1$, $N_2$, l, $\xi$ and functions u and v.

From Eqs. (36) and (40), the approximation formulas for integral (35) can be obtained as

$$\begin{aligned} \begin{aligned}&I_{\alpha ,\beta } \approx R^{(3-\alpha -\beta )} \sum \limits _{m=N_1}^{m=N_2}W_m^{\alpha ,\beta }(p,q,\varphi ,\alpha ,\beta )f\left( \frac{\lambda _m}{R}\right) ,\\&\lambda _m=\mathrm {e}^{ml+\xi }. \end{aligned} \end{aligned}$$

(41)

For the given input function $v(t)=8\mathrm {e}^{-t^2}\left( t^5-4t^3+2t\right)$, the output functions for different kernels can be expressed analytically by Gaussian and error functions. Inspired by (Anderson 1979), the parameters used are $l=0.2$, $\xi =0.0964$, $M=512$, $N_1=-250$, $N_2=200$. In computational experiments, these parameters provided appropriate accuracy in calculation of $\widehat{B}_n^m$.