# On a Fast Direct Solver for High Order Integral Equations: Results for Low-to-Middle Frequencies

Damiano Franzò, Adrien Merlini, Clément Henry, Francesco Andriulli

► To cite this version:

Damiano Franzò, Adrien Merlini, Clément Henry, Francesco Andriulli. On a Fast Direct Solver for High Order Integral Equations: Results for Low-to-Middle Frequencies. 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.499-500, 10.1109/AP-S/INC-USNC-URSI52054.2024.10687162. hal-04801333

HAL Id: hal-04801333

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

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
# On a Fast Direct Solver for High Order Integral Equations: Results for Low-to-Middle Frequencies

Damiano Franzò$^{(1)}$, Adrien Merlini$^{(2)}$, Clément Henry$^{(2)}$, and Francesco P. Andriulli$^{(1)}$

$^{(1)}$ Department of Electronics and Telecommunications, Politecnico di Torino, Turin, Italy

$^{(2)}$ Microwave Department, IMT Atlantique, Brest, France

*Abstract*—In this work, we present a fast direct solver for integral equations-based electromagnetic scattering that is compatible with spatial discretization of arbitrary order. The proposed approach leverages a suitably preconditioned combined field formulation and involves the regularization of the elliptic spectrum through the extraction of a specifically chosen equivalent circulant problem. The procedure yields a compressed representation of the inverse system matrix in the form of a matrix amenable to fast inversion and a single-skeleton low-rank perturbation. In addition to proposing the fast direct solution framework, this contribution presents the algorithms and explores the intricacies that enable the high-order discretization of the underlying integral equations and is supported by numerical results illustrating its effectiveness.

## I. INTRODUCTION

Integral equation (IE) methods are widely used to solve radiation and scattering problems across a wide range of frequencies and common equations are the electric, magnetic, and combined field integral equations (EFIE, MFIE, and CFIE respectively). One of their attractive features is that they only require the discretization of the boundary of the scatterer. Nonetheless, the matrices that result from the discretization procedure are dense and ill-conditioned in several regimes. Solving classical IE schemes typically has a computational complexity greater than quadratic for high frequency regimes. One way to lower this computational complexity are fast matrix-vector products (MVP) used with iterative solvers for efficient solutions. Nevertheless, this approach becomes less attractive when the matrix systems need to be solved for numerous right-hand sides (RHS), as is required, for instance, in computing backscattering cross-sections. To overcome this, different fast-direct solver (FDS) methods have been implemented, such as [3] or references therein. Recently, a new FDS based on a preconditioned combined field integral equation (CFIE) and that enhances compressibility of the system matrix by extracting a circulant matrix, has been presented [1]. The overall scheme yields a single skeleton form to be stored for the whole extracted part, such that accurate solutions for multiple RHSs can be obtained efficiently. This method however, did not allow for high order discretization, thereby limiting its ability to harness their faster convergence rates.

Our work extends the previously introduced fast direct solver scheme [1], expanding its capabilities to include high-order functions. The new scheme yields a representation of the inverse system matrix in the form of a matrix amenable to fast inversion and a single-skeleton low-rank perturbation. The new algorithms and corresponding theoretical framework

are presented in this contribution and are accompanied by numerical results, showing the effectiveness of the high-order scheme.

## II. BACKGROUND AND NOTATION

Consider a perfectly electrically conducting (PEC) object delimited by a 2-dimensional convex contour $\gamma$, excited by a TE time-harmonic incident electromagnetic field ($E_{\text{inc}}$, $H_{\text{inc}}$) of angular frequency $\omega$. This object is situated within a medium characterized by a permittivity $\varepsilon$, permeability $\mu$, impedance $\eta = \sqrt{\mu/\varepsilon}$, and wavenumber $k = \omega\sqrt{\mu\varepsilon}$. The CFIE that establishes the relationship between the induced tangential current density $j$ on the surface $\gamma$ and the incident electric and magnetic fields is

$$ \frac{1}{ik}\mathcal{N}^k(j_t) + \left(\frac{\mathbf{I}}{2} - \mathcal{D}^k\right) (j_t) = -\frac{1}{\eta}E_t^{\text{inc}} - H_z^{\text{inc}}, \quad (1) $$

where $\mathcal{N}^k(j_t) := -\frac{\partial}{\partial n} \int_{\gamma} \frac{\partial}{\partial n'} g_k(\mathbf{r}, \mathbf{r}') j_t(r') dr'$, $\mathcal{D}^k(j_t) := \int_{\gamma} \frac{\partial}{\partial n'} g_k(\mathbf{r}, \mathbf{r}') j_t(r') dr'$, and $g_k(\mathbf{r}, \mathbf{r}') := i/4H_0^{(1)}(k ||\mathbf{r} - \mathbf{r}'||)$, where $H_0^{(1)}$ indicates the Hankel function of the first kind. The boundary elements used to expand the current density $j_t$ are scalar piecewise polynomials $\lambda_i^p$, forming a complete space of order $p$, in which cell-to-cell continuity is appropriately enforced. After Galerkin testing, the discrete system reads

$$ \frac{1}{ik}\mathbf{N}^k\mathbf{j} + \left(\frac{\mathbf{G}}{2} - \mathbf{D}^k\right)\mathbf{j} = -\frac{1}{\eta}\mathbf{e} - \mathbf{h}, \quad (2) $$

with $[\mathbf{e}]_i = \int_{\gamma} \lambda_i^p(\mathbf{r}) E_t^{\text{inc}}(\mathbf{r}) d\mathbf{r}$, $[\mathbf{h}]_i = \int_{\gamma} \lambda_i^p(\mathbf{r}) H_t^{\text{inc}}(\mathbf{r}) d\mathbf{r}$, $[\mathbf{N}^k]_{ij} = \langle \lambda_i^p, \mathcal{N}^k \lambda_j^p \rangle$, $[\mathbf{D}^k]_{ij} = \langle \lambda_i^p, \mathcal{D}^k \lambda_j^p \rangle$, $[\mathbf{G}^k]_{ij} = \langle \lambda_i^p, \lambda_j^p \rangle$, with $\langle f, g \rangle = \int_{\gamma} f(\mathbf{r}) g(\mathbf{r}) d\mathbf{r}$.

## III. METHODS

The matrix system (2) can be preconditioned with a modified version of the preconditioner proposed in [2], which results in the following discrete system

$$ \mathbf{C}\mathbf{j} = -\frac{ik}{\eta}\mathbf{S}^{\tilde{k}}\mathbf{G}^{-1}\mathbf{e} - \left(\frac{\mathbf{G}}{2} + \mathbf{D}^{\tilde{k}}\right)\mathbf{G}^{-1}\mathbf{h}, \quad (3) $$

with

$$ \mathbf{C} = \mathbf{S}^{\tilde{k}}\mathbf{G}^{-1}\mathbf{N}^k + \left(\frac{\mathbf{G}}{2} + \mathbf{D}^{\tilde{k}}\right)\mathbf{G}^{-1}\left(\frac{\mathbf{G}}{2} - \mathbf{D}^k\right), \quad (4) $$

where $\tilde{k} := k + 0.4ik^{1/3}a^{-2/3}$, and $a$ is a suitable average of the radii of curvature along $\gamma$.
The FDS scheme is based on the separation of the system matrix **C** into a main part for which some fast algebra arithmetic is available, and a compressible rank-deficient reminder. The main part arises from the discretization of an equivalent circular problem generating a matrix **C**<sub>c</sub> that entails the computation of the same operators appearing in (4) on a circular curve γ<sub>c</sub> of the same perimeter as γ. Finally, **C**<sub>c</sub> := **C** − **C**<sub>c</sub> represents the remaining extracted part that exhibits a rank deficiency.

Because $\mathbf{C}_c$ is built out of products and linear combination of matrices with structured circularity, generically designated as $\mathbf{O}_c$ in the following, it can be stored an multiplied to a vector in quasi-linear complexity. In particular, when discretized with order one functions, $\mathbf{O}_c$ is circulant, which enables fast matrix operations through the fast Fourier transform (FFT). However, for higher basis order $p$, these matrices are no longer circulant. To recover fast manipulations of $\mathbf{O}_c$, we organize the matrix elements into a matrix composed of circulant blocks

$$ \mathbf{O}_c = \begin{bmatrix} \mathbf{O}_c^{11} & \mathbf{O}_c^{12} & \cdots \\ \vdots & \vdots & \\ \mathbf{O}_c^{p1} & \mathbf{O}_c^{pp} \end{bmatrix}, \quad (5) $$

where each block $\mathbf{O}_c^{ij}$ is internally circulant. Despite not being exactly circulant, with suitable generalization, this structure can be exploited to enable fast matrix operations via FFT operations performed at a block level. The main idea stems from the possibility to diagonalize $\mathbf{O}_c$ in quasi linear time with the following expression

$$ \mathbf{K}_{(p,n)}^T \tilde{\mathbf{F}}_p^* \mathbf{O}_c \tilde{\mathbf{F}}_p \mathbf{K}_{(p,n)} = \begin{bmatrix} \mathbf{B}^{11} & 0 & \cdots \\ \vdots & \vdots & \\ 0 & \mathbf{B}^{nn} \end{bmatrix}, \quad (6) $$

where $\mathbf{B}^{ij} \in \mathbb{C}^{p \times p}$, $\mathbf{K}_{(p,n)}$ is the commutation matrix, $n$ is the number of vertices, and $\tilde{\mathbf{F}}_p = \mathbf{I}_p \otimes \mathbf{F}$, with $\mathbf{F}$ being the discrete Fourier transform matrix, $\otimes$ the Kronecker product and $\mathbf{I}_p$ an identity matrix with $p$ diagonal elements. The matrix $\tilde{\mathbf{F}}_p$ can be implicitly manipulated using the FFT. Consequently, by noting also that $\mathbf{K}_{(p,n)}$ benefits from permutation matrix properties (such as $\mathbf{K}_{(p,n)}^{-1} = \mathbf{K}_{(p,n)}^T$).

To finalize the scheme, the extracted part is compressed into a non-hierarchical low-rank form $\mathbf{C}_e \approx \mathbf{U}\mathbf{V}^T$, and the inversion of **C** is performed through the Woodbury identity

$$ \mathbf{C}^{-1} = \mathbf{C}_c^{-1} - \mathbf{C}_c^{-1} \mathbf{U} (\mathbf{I} + \mathbf{V}^T \mathbf{C}_c^{-1} \mathbf{U})^{-1} \mathbf{V}^T \mathbf{C}_c^{-1}, \quad (7) $$

which can efficiently computed by exploiting the low-rank form and the fast algebra operations described above.

## IV. NUMERICAL RESULTS

To validate the new direct solution scheme, it is applied on an ellipse of with a semi-major axis $a = \cosh(1)$ and a semi-minor axis $b = \sinh(1)$. The numerical results reported in Fig. 1 and Fig. 2 show (i) for an appropriate compression threshold $\varepsilon$ the proposed scheme delivers high accuracy solution and (ii) that increasing the order $p$ does not

<figure><img src="image_2.png" /><figcaption>Fig. 1: Current accuracy for $k = 40 \text{ m}^{-1}$, where $\varepsilon$ is the compression threshold set on $\mathbf{C}_e \approx \mathbf{U}\mathbf{V}^T$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2: $\varepsilon$-rank of $\mathbf{C}_e \approx \mathbf{U}\mathbf{V}^T$ for $\varepsilon = 10^{-2}$ and $k = 40 \text{ m}^{-1}$.</figcaption></figure>

significantly impact the $\varepsilon$-rank (and compressibility) of the extracted part This shows that our FDS scheme is capable of preserving the better accuracies of high order schemes, while maintaining a reduced complexity, and simplicity of the compressed representation.

## 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), and by the Italian PNRR action funded by the European Union through the NextGenerationEU Initiative, Program RESTART, Structural Project S12 "Smart propagation environments" (PNRR M4C2, Investimento 1.3 - Avviso n. 341 del 15/03/2022 - PE0000001 RESearch and innovation on future Telecommunications systems and networks, to make Italy more smart (RESTART) - CUP E13C22001870001)

## REFERENCES

[1] D. Consoli, C. Henry, A. Dély, L. Rahmouni, John Erik Ortiz Guzman, T. L. Chhim, S. B. Adrian, A. Merlini, and F. P. Andriulli "On the Fast Direct Solution of a Preconditioned Electromagnetic Integral Equation", International Conference on Electromagnetics in Advanced Applications (ICEAA), Cape Town, South Africa, 2022, pp. 193-195

[2] F. P. Andriulli, I. Bogaert and K. Cools, "On the high frequency behavior and stabilization of a preconditioned and resonance-free formulation," 2015 International Conference on Electromagnetics in Advanced Applications (ICEAA), Turin, Italy, 2015, pp. 1321-1324

[3] A. Sharshевsky, Y. Brick, and A. Boag, "Direct Solution of Scattering Problems Using Generalized Source Integral Equations," IEEE Transactions on Antennas and Propagation, vol. 68, no. 7, pp. 5512-5523, Jul. 2020.