# A Tree-Based Scheme for the Numerically Stable Compression of Integral Equations

Giulio Cosentino, Clément Henry, Maxime Monin, Adrien Merlini, Francesco

Andriulli

► To cite this version:

Giulio Cosentino, Clément Henry, Maxime Monin, Adrien Merlini, Francesco Andriulli. A Tree-Based Scheme for the Numerically Stable Compression of Integral Equations. 2024 IEEE International Symposium on Antennas and Propagation and INC/USNC-URSI Radio Science Meeting (AP-S/INC-USNC-URSI), Jul 2024, Firenze, France. pp.503-504, 10.1109/AP-S/INC-USNC-URSI52054.2024.10686275. hal-04801504

HAL Id: hal-04801504

https://hal.science/hal-04801504v1

Submitted on 25 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Tree-Based Scheme for the Numerically Stable Compression of Integral Equations

Giulio Cosentino(1), Clément Henry(2), Maxime Y. Monin(1), Adrien Merlini(2), and Francesco P. Andriulli(1)

**Abstract**—When compressing the electric field integral equation with the adaptive cross approximation or related fast solvers, the resulting matrix can often be affected by numerical cancellations at low frequency in the computation of the scalar potential. This work theoretically studies the source of this instability and proposes a solution based on a tree scheme. Numerical results corroborate and confirm the theory and show the efficacy of the proposed solution.

## I. INTRODUCTION

The electric field integral equation (EFIE), solved by the boundary element method (BEM), becomes a very competitive solver for electromagnetic analyses when associated with a fast solver such as the multilevel adaptive cross approximation (MLACA) [1]. At low frequency, however, this formulation is affected by ill-conditioning and loss of accuracy in the solution process, which is often referred to as the low-frequency breakdown. An efficient and widely spread remedy to this limitation leverages the quasi-Helmholtz projectors that enable the proper rescaling of the components of the EFIE system matrix at low frequencies. A crucial step when manipulating these projectors is to enforce the divergence-free nature of the solenoidal functions employed in the quasi-Helmholtz decomposition, which is a challenging operation when using a standard compression scheme at the basis functions level of the EFIE. The reason behind this is the lack of orthogonality between the compressed hypersingular operator ($\mathbf{T}_{\Phi}$) of the EFIE and the solenoidal functions. A way to address this is to compress the kernel at quadrature points level and thereby leaving the divergence part of $\mathbf{T}_{\Phi}$ as uncompressed sparse matrices. In addition to proposing a stable compression of the EFIE at very low frequencies, we show in this contribution, with the adequate mathematical framework, that this method suffers from a loss of accuracy in the mapping between the compressed kernel and $\mathbf{T}_{\Phi}$. Moreover, we introduce a remedy based on a hierarchical tree that prevents the numerical loss of accuracy caused by the divergence operator and retains the low frequency stability of the preconditioned EFIE. Numerical results will corroborate the theory and validate the accuracy of this new scheme.

## II. BACKGROUND AND NOTATION

Let $E^i$ be a time-harmonic incident electric field impinging on a perfectly electrically conducting scatterer delimited by the simply connected and closed boundary $\Gamma \subset \mathbb{R}^3$ with outward pointing normal $\hat{n}$. The induced surface current density can be obtained by solving the EFIE, which reads

$\mathcal{T}\mathbf{J} = ik\mathcal{T}_A + \frac{1}{ik}\mathcal{T}_{\Phi} = -\hat{n} \times \mathbf{E}^i$ where $\mathcal{T}_A\mathbf{J}(\mathbf{r}) = \hat{n} \times \int_{\Gamma} G(\mathbf{r}, \mathbf{r}')\mathbf{J}(\mathbf{r}')dS'$ is the vector potential operator, $\mathcal{T}_{\Phi}\mathbf{J}(\mathbf{r}) = -\hat{n} \times \nabla_{\Gamma} \int_{\Gamma} G(\mathbf{r}, \mathbf{r}') \text{div}_{\Gamma} \mathbf{J}(\mathbf{r}')dS'$ is the scalar potential operator and $k$ is the wavenumber. Here, $G(\mathbf{r}, \mathbf{r}') = e^{ik\|\mathbf{r}-\mathbf{r}'\|}/(4\pi \|\mathbf{r}-\mathbf{r}'\|)$ is the Green's function, $\nabla_{\Gamma}$ and $\text{div}_{\Gamma}$ are the surface gradient and surface divergence, respectively.

The EFIE is discretized using a triangular approximation of $\Gamma$ with average edge length $h$ on which the Rao-Wilton-Glisson (RWG, [2]) basis functions $\{f_i\}_{i=1}^N$ are defined to expand the current as $\mathbf{J} \approx \sum_{i=1}^N [\mathbf{j}_i]f_i$. The EFIE is then tested with rotated RWGs $\{\hat{\mathbf{n}} \times f_i\}$, which yields the linear system $\mathbf{T}\mathbf{j} = (ik\mathbf{T}_A + \frac{1}{ik}\mathbf{T}_{\Phi})\mathbf{j} = -e^i$ where $[\mathbf{T}_A]_{ij} := \langle \hat{\mathbf{n}} \times f_i, \mathcal{T}_A f_j \rangle_{\Gamma}$, $[\mathbf{T}_{\Phi}]_{ij} := \langle \hat{\mathbf{n}} \times f_i, \mathcal{T}_{\Phi} f_j \rangle_{\Gamma}$, $[e^i]_i := \langle \hat{\mathbf{n}} \times f_i, \hat{\mathbf{n}} \times \mathbf{E}^i \rangle_{\Gamma}$ with $\langle \mathbf{f}, \mathbf{g} \rangle_{\Gamma} = \int_{\Gamma} \mathbf{f}(\mathbf{r}) \cdot \mathbf{g}(\mathbf{r})dS$. We also define the single layer potential operator $Sf(\mathbf{r}) = \int_{\Gamma} G(\mathbf{r}, \mathbf{r}')f(\mathbf{r}')dS'$ that can be used for an alternative definition of $[\mathbf{T}_{\Phi}]_{ij} = \langle \text{div}_{\Gamma} f_i, S \text{div}_{\Gamma} f_j \rangle_{\Gamma}$. The interaction between two RWG functions $f_i$ and $f_j$ consists in a weighted sum of interactions between the four triangular cells $T_i^{\pm}$, $T_j^{\pm}$ on which the RWGs are defined, which can be integrated by Gaussian quadrature with a set of points $\{\mathbf{r}_i\}_{i=1}^{N_{\text{pts}}}$.

## III. NUMERICAL CANCELLATION: THEORETICAL ANALYSIS

To compress $\mathbf{T}$ at quadrature points level, $\mathbf{T}_A$ and $\mathbf{T}_{\Phi}$ are decomposed as $\tilde{\mathbf{T}}_A = \mathbf{X}_A \mathbf{K} \mathbf{X}_A^{\top}$ and $\tilde{\mathbf{T}}_{\Phi} = \mathbf{X}_{\Phi} \mathbf{K} \mathbf{X}_{\Phi}^{\top}$, where the point-point interactions are stored as $[\mathbf{K}]_{ij} := G(\mathbf{r}_i, \mathbf{r}_j)$ and the sparse mapping matrices $\mathbf{X}_{\Phi}$ and $\mathbf{X}_A$ combine the Gaussian weights and cells coefficients. However, $\text{div}_{\Gamma} f_i|_{T_i^{\perp}}$ and $\text{div}_{\Gamma} f_i|_{T_i^{\parallel}}$ are constant and of opposite value and the above-mentioned construction of $\tilde{\mathbf{T}}_{\Phi}$ results in numerical cancellation. We will now prove the following result:

**Proposition 1.** Given two RWG basis functions $f_i$ and $f_j$ whose cells $T_i^{\pm}$, $T_j^{\pm}$ have aligned centroids, suppose that $[\mathbf{T}_{\Phi}]_{ij}$ is computed via Gaussian quadrature with one quadrature point per cell: $c_i^{\pm}, c_j^{\pm}$. Let $\Delta r = \|c_i^{\pm} - c_j^{\pm}\| = \|c_j^{\mp} - c_j^{\pm}\|$ and $R = \|c_j^{\pm} - c_i^{\pm}\|$, and suppose $\Delta r \ll R$ (see Fig. 1a). If the Green function is compressed at tolerance $\epsilon$, then the upper bound on the error on the resulting compressed element $[\mathbf{T}_{\Phi}]_{ij}$ is

$$ \frac{4\epsilon}{h^2} \left| k^2 + \frac{2ik}{R} - \frac{2}{R^2} \right|^{-1} \approx \frac{4\epsilon}{h^2} \max \left\{ \frac{2}{R^2}, k^2 \right\}^{-1} . \quad (1) $$

*Proof.* Let $q(R^{\pm\pm}) := G(c_i^{\pm}, c_j^{\pm}) = e^{ikR^{\pm\pm}}/(4\pi R^{\pm\pm})$ be the point-point interaction. Then, using $\text{div}_{\Gamma} f_i|_{T_i^{\perp}} =$

<footnote>(1) Politecnico di Torino, Turin, Italy</footnote>

<footnote>(2) IMT Atlantique, Brest, France</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: (a) The triangular support of a pair of interacting basis functions. (b) Relative error of the far interaction between two RWGs computed from the compression of the cells. (c) Relative error of a far block of the compressed $\mathbf{T}_{\Phi}$ for a sphere.</figcaption></figure>

$-\operatorname{div}_{\Gamma} f_i|_{T_i^-} > 0$ and the definition of $[\mathbf{T}_{\Phi}]_{ij}$, we expand the ratio $[\mathbf{T}_{\Phi}]_{ij}/q(R^{++})$ as $h \to 0$

$$ 
\begin{aligned}
\frac{[\mathbf{T}_{\Phi}]_{ij}}{q(R^{++})} &= \frac{4\pi R}{e^{ikR}} (q(R^{++}) - q(R^{+-}) - q(R^{-+}) + q(R^{--})) \\
&= \frac{4\pi R}{e^{ikR}} \left( \frac{2e^{ikR}}{4\pi R} - \frac{e^{ik(R+\Delta r)}}{4\pi(R+\Delta r)} - \frac{e^{ik(R-\Delta r)}}{4\pi(R-\Delta r)} \right) \\
&= \left( k^2 + \frac{2ik}{R} - \frac{2}{R^2} \right) h^2 + \mathcal{O}(h^4),
\end{aligned}
\quad (2) $$

where we have used the fact that $\Delta r = \mathcal{O}(h)$. Now, assuming that $q$ is compressed with relative error $\epsilon$, then its absolute error is $|q(R^{\pm\pm}) - \tilde{q}(R^{\pm\pm})| = \epsilon|q(R^{\pm\pm})|$. Since the maximum error of a sum equals the sum of the maximum errors of the terms, then the error on $[\tilde{\mathbf{T}}_{\Phi}]_{ij}$ is less than $4\epsilon|q(R^{\pm\pm})|$. By dividing by $|[\tilde{\mathbf{T}}_{\Phi}]_{ij}|$ and we obtain an upper bound of the relative error on $[\tilde{\mathbf{T}}_{\Phi}]_{ij}$:

$$ \frac{4\epsilon|q(R^{\pm\pm})|}{|[\tilde{\mathbf{T}}_{\Phi}]_{ij}|} = \frac{4\epsilon}{h^2} \left| k^2 + \frac{2ik}{R} - \frac{2}{R^2} \right|^{-1}. \quad (3) $$

Finally, since the term $2/R^2$ dominates for $R < 1/k$ and the term $k^2$ dominates for $R > 2/k$, the bound can be approximated as $4\epsilon/h^2 \max\{2/R^2, k^2\}^{-1}$. $\square$

This proposition suggests that the loss of significant digits on $\tilde{\mathbf{T}}_{\Phi}$ is an issue for far basis functions at low frequency and for refined meshes.

## IV. NUMERICAL CANCELLATION: SOLUTION

The solution we propose relies on the compression of entries with the same scaling as the entries of $\mathbf{T}_{\Phi}$. Using an octree in $\mathcal{O}(N_{pts} \log N_{pts})$ complexity, we construct the minimum spanning tree linking the quadrature points and identify the edges by the index of their end points, after arbitrarily choosing a root. An edge-edge matrix $\mathbf{K}_e$ is considered, with entries $[\mathbf{K}_e]_{ij} = G(c_i^+, c_j^+) - G(c_i^-, c_j^-) - G(c_i^-, c_j^+) + G(c_i^-, c_j^-)$ where $c_i^\pm$ and $c_j^\pm$ are the origin and end points of the edges $i$ and $j$. With its definition $\mathbf{K}_e = \mathbf{Y}\mathbf{K}\mathbf{Y}^\top$, where $\mathbf{Y}$ is the sparse matrix of entries $[\mathbf{Y}]_{ij} = 1$ if $r_j$ is the origin of edge $i$, $-1$ if $r_j$ is the end of edge $i$, and $0$ otherwise. $\mathbf{Y}^{-1}$ and its transpose are not sparse, but their matrix-vector products

can be performed in linear complexity thanks to the tree structure, starting from the root. After compressing $\mathbf{K}_e$, we have $\tilde{\mathbf{T}}_{\Phi} = \mathbf{X}_{\Phi}\mathbf{Y}^{-1}\tilde{\mathbf{K}}_e(\mathbf{Y}^{\top})^{-1}\mathbf{X}_{\Phi}^{\top}-\tilde{\mathbf{T}}_{\Phi}^{\text{near}}+\mathbf{T}_{\Phi}^{\text{near}}$, where the sparse pre-correction accounts for the incompressible singularity. The row and column requiring the root vertex are filled as points before the application of the $\mathbf{X}_{\Phi}$ matrices.

## V. NUMERICAL RESULTS

In Fig. 1b we plotted the relative error on the far interaction of two equal, equally oriented, and coplanar RWG basis functions with $\Delta r = 1$ m at a frequency of 1.5 MHz. Both the vector and scalar potential entries are computed from the four cells interactions with one Gaussian quadrature point per cell. The compression at the desired tolerance of $10^{-3}$ is simulated by rounding the cell-cell entries, giving the sawtooth profile. This corroborates Proposition 1, as $\mathbf{T}_{\Phi}$ suffers from cancellation, unlike $\mathbf{T}_A$, and the resulting error scales as predicted and changes behaviour between $R = 1/k$ and $R = 2/k$.

The tree-based solution has been tested at a frequency of 1 MHz on a sphere of radius 1 m discretized at $h = 10$ cm, resulting in 4320 basis functions. In Fig. 1c we show the relative errors of a $1080 \times 1076$ contiguous block of $\mathbf{T}_{\Phi}$ corresponding to the interactions between two far regions of the sphere after compressing the whole matrix. While using an MLACA on $\mathbf{K}$ implies the loss of significant digits, our strategy remains unaffected by numerical cancellation, resulting in an error that follows the tolerance.

## ACKNOWLEDGEMENT

This work has been funded in by the European Innovation Council (EIC) through the European Union's Horizon Europe research Programme under Grant 101046748 (Project CEREBRO).

## REFERENCES

[1] K. Zhao, M. Vouvakis, and J-F Lee, "The adaptive cross approximation algorithm for accelerated method of moments computations of EMC problems," *IEEE transactions on electromagnetic compatibility*, vol. 47(4), pp. 763-773, 2005.

[2] S. Rao, D. Wilton and A. Glisson, "Electromagnetic scattering by surfaces of arbitrary shape," *IEEE Transactions on Antennas and Propagation*, vol. 30, no. 3, pp. 409-418, May 1982.