An Efficient Method for Solving Problems of Acoustic Scattering on Three-Dimensional Transparent Structures

Abstract

The article contains a study of methods for solving integral equations in the context of acoustic problems. The methodology considered is applied to describe acoustic wave propagation and scattering. Efficient discretization methods are used together with iterative methods to solve the operator equations, including an apparatus for fast multiplication of the resulting post-discretization Toeplitz matrices by a vector using the fast Fourier transform. The theoretical analysis of the proposed numerical algorithm demonstrates its efficiency in terms of the required number of arithmetic operations and the memory footprint of the computing system. The presented numerical simulation demonstrates the possibility of solving the problem of acoustic wave propagation in transparent media using the proposed methods. A visualization of the obtained solutions for a practical problem with a high level of discretization of the solution volume domain is also presented.

1. Introduction

The initial formulation for solving problems of mathematical physics using integral equations refers to the works of Fredholm, who first proposed the study and solution of classical problems of mathematical physics using the transformation of the original differential equations to the formulation of the problem for solving related integral equations, later called Fredholm equations. The solution of problems of mathematical physics, both in the form of differential and integral equations. at that time was impossible using numerical methods and methods of discretization of the initial problem via a grid of high dimensionality in the number of elementary partitions, as such problems were associated with impossible manual calculations. With the development of computer technology appeared the formulation and the first methods and algorithms for distributed computations over differential equations, allowing the modeling of real physical processes, and integral equations were used less frequently, due to problems with the complexity of calculations, in the solution of systems of linear algebraic equations (SLEs) over fully filled operator matrices.

Mikhlin, in his works, investigated different variations of problem statements for solving integral equations in the theory of elasticity [1], and also investigated the possibilities of obtaining analytical and numerical solutions of singular integral equations [2,3] related to the problems of mathematical physics. Colton and Kress considered surface and volume problems of acoustics and electrodynamics based on integral equations, as well as their numerical formulations [4]. In these papers, however, due to the reasons described above, there were no numerical results accompanying test problems. Miller, in a series of papers [5], also investigated modern formulations for the solution of electrodynamics problems on the basis of computational methods for solving bulk singular integral equations, using features of kernels of integral operators based on Green’s functions [4], satisfying solutions for the volumetric Helmholtz equation.

Mathematical formulation and development of effective computational methods for solving problems of mathematical physics are the most important areas of science and technology development today. Modern problems are associated with large requirements for accuracy and speed of computation, in view of which the requirements of the dimensionality of discretization of problems increase. The classical formulation of the solution of differential and integral equations, as well as general computational methods for solution, are becoming insufficient to meet the needs of experimenters, researchers and manufacturers. This aspect pushes us to the necessity of developing effective methods for solving specific problem formulations based on a priori factual knowledge, assumptions about the structure of the solution, and constraints on the conditions at the boundary of the domain [6].

Overcoming these limitations, recent works on computational methods, applicable to the solution of differential equations that reduce to solving systems of linear equations (SLEs) with high-dimensional sparse operator matrices, have proposed a series of approximate methods capable of addressing the problem of increasing the discretization of the solution domain, based on homotopies. Among these are the homotopy method described in detail in [7], the multigrid method [8], the multigrid–homotopy method [9], and the wavelet multiscale method [10]. These methods are traditionally used in tasks related to solving problems of mathematical physics, including, for example, [11,12,13], and conducting numerical experiments with models based on partial differential equations.

Homotopy methods [7] for approximate solutions of large-dimensional SLEs, resulting from the discretization of the solution domain of the corresponding differential equation, allow the obtaining of solutions with specified accuracy on large grids. The main idea of the method is to use the concept of homotopy for a gradual transition from a simplified version of the problem to its original, more complex form, while following a path that maintains connectivity between solutions. The homotopy grid method is especially useful in cases where direct solving of the original problem is difficult due to its complexity or the presence of multiple local minima/maxima.

The algebraic multigrid method [8] represents a class of algorithms for solving systems of linear equations, particularly effective for large sparse matrices arising from the discretization of differential equations. This method does not require explicit knowledge of the geometry of the problem and is built directly on the basis of the equation coefficients, making it applicable to a wider class of problems. The key idea of multigrid methods is to use iterations at different levels of discretization (grids) to accelerate convergence to the solution. Simplified versions of the original problem are solved on coarse grids, effectively reducing error components, which decrease slowly on finer grids. This is achieved by analyzing the structure of the system matrix and grouping variables in such a way that the problem dimension decreases at each subsequent level. Various coarsening strategies, for example, based on the strength of connections between matrix elements, are used for this purpose.

Combining these two approaches, the multigrid–homotopy method [9] leverages the advantages of both in order to solve complex nonlinear problems. The homotopy method allows effective finding of an initial approximation for a nonlinear problem or a series of such approximations, while the multigrid approach ensures rapid convergence to an accurate solution at each step of homotopy. This method is particularly useful for problems where direct application of the multigrid method to a nonlinear system proves inefficient due to the complexity of the problem or lack of a good initial approximation. Combining the homotopy approach, to find an initial approximation, with the multigrid method for its refinement can significantly improve the efficiency of the iterative process.

All the above methods can adequately solve high-dimensional problems related to setting tasks in terms of differential equations; however, they may be difficult to apply to tasks based on integral equations. This is due to the complexity of calculations with fully populated operator matrices that arise when applying various discretization methods to integral formulations of mathematical physics problems.

The formulation of problems in terms of integral equations in mathematical physics is also important, as it allows the finding of the desired solution with fewer restrictions both on the form of the solution function itself and on the number of initial and boundary conditions. Solving integral equations on a grid of large dimensionality is associated with the use of fast iterative methods capable of solving full systems of linear equations in fewer iterations with a given accuracy. For this, it is necessary to introduce special discretization grids for a more accurate approximation [14] of the original solution domain, both in volume and at the boundaries between different media. This work presents one such technique for introducing a special discretization grid for volumetric integral equations, allowing for the solution of such equations using efficient matrix–vector multiplications in the process of any iterative method of solving systems of linear equations based on fast Fourier transform (FFT).

The development of integral equations theory to date is described by a large number of works in the field of proving existence and uniqueness theorems for solutions of the problems of mathematical physics, as well as the development and study of effective numerical methods for their solution. In electrodynamics, Smirnov Yu.G. proved a large number of existence and uniqueness theorems for problems of wave propagation, diffraction and polarization using the Fredholm theory for integral operators and the theory of pseudo-differential operators [15,16,17]. In hydrodynamics and aerodynamics, integral equations can be used in modeling stationary flows, shown in the works of Setukha [18].

In this paper we consider the problems of wave propagation in three-dimensional bounded transparent structures with inhomogeneous refraction based on the Fredholm integral equation of the second kind. Modeling problems of acoustic wave propagation in an inhomogeneous medium has a wide range of applications in various fields. For example, in oceanography, such models are used to investigate underwater sound channels and to study the influence of various factors on sound propagation in water [19]. In the process of developing printed circuit boards, defects, such as delamination, can occur during the lamination process, which can be detected by the acoustic emission method discussed in [20]. In electrodynamics, models of propagation and scattering of electromagnetic waves, built on integral equations, help to make calculations for complex processes of radiation effects and load on technical devices, necessary to their design [21,22]. Such models are also used in medical diagnostics and industrial acoustics.

Numerical methods for solving integral equations, along with the growth of computing power and memory capacity, allow the solution of more and more complex structural problems with a high degree of discretization, due to optimization and the possibility of efficient parallelization of the computations of problems with the full operator matrix. Mathematical modeling of real physical problems on the basis of the apparatus for effective methods for their solution, shown in this paper, will be useful in real-time calculations for digital twins in many physical processes, natural phenomena, technical devices and their manufacturing processes. The main part of the research, problem statements, and computational methods used in the article grew out of the theoretical basis of the research on problem statements, methods, and the algorithms for their numerical solution given in the corresponding book [23].

2. Materials and Methods

In this study, we consider the integral equation within a bounded region of $Q$ of Euclidean space $E_3$ [24]:

$$\left(1+\alpha \eta \left(x\right)\right) u\left(x\right)+\underset{Q}{\int }\frac{K\left(x-y\right)}{R^m}\eta \left(y\right) u\left(y\right) dy=u^0\left(x\right),x\in Q,m\le 3.$$

(1)

In Equation (1), $x=\left(x_1,x_2,x_3\right),$ $y=\left(y_1,y_2,y_3\right)$ are the points belonging to the bounded region Q; $R=\left|x-y\right|$ is the distance between the points of the region; and $\alpha , \eta , K, u^0$ are known functions, with $K(x-y$) a differentiable coordinate function; u(x) is an unknown function.

Further, suppose that Equation (1) has a unique solution in the corresponding function space. Numerical methods turn out to be uniquely possible for solving Equation (1). In this case, Equation (1), by applying the Galerkin method or the collocation method, is approximated to an SLE with a fully filled matrix.

The equation in the form of Equation (1) describes a number of important applied problems, including propagation and scattering of acoustic waves on transparent inhomogeneous obstacles [4]. In this case, $m<3$, $\alpha =0$ and the other functions included in Equation (1) are scalar. Then, the equation is a classical Fredholm integral equation of the 2nd kind, which can be written in the following form:

$$\begin{gathered} u\left(x\right)+k^2{\int }_{Q}G\left(R\right)\widehat{n}\left(x\right)u\left(y\right)dy=u^0\left(x\right); \\ \widehat{n}\left(x\right)=\eta \left(x\right)-1; G\left(R\right)=-\frac{\exp \left(ikR\right)}{4\pi R}; u^0\left(x\right)=-{\int }_{Q}G\left(R\right)f\left(y\right)dy; x\in Q, \end{gathered}$$

(2)

where $\eta \left(x\right)$ is the function characterizing the acoustic refraction of the definition area, and $\eta \left(x\right)=1, x\notin Q$; $k$ is the wave number characterizing the properties of the simulated external radiation; $G\left(R\right)=-K(x-y)/R$ is the kernel of the integral equation, which is the solution of the corresponding differential formulation of the Helmholtz equation in three-dimensional space; $f\left(x\right)$ is the modeled external radiation function; and $u\left(x\right)$ is the unknown scalar potential field characterizing the stress for the $x\in Q$.

Previously, in [25], for the integral formulation of Equation (2), the theorem of the existence and uniqueness of the solution of the corresponding differential equation was proved under the condition of radiation at infinity, as well as under some conditions for refraction $\eta \left(x\right)$.

Equation (2) describes both the problems of acoustic wave propagation in bulk transparent media with a homogeneous scalar refraction index $\eta \left(x\right)$ and the problems of wave scattering at the boundary of media with different $\eta \left(x\right)$, where the function itself is piecewise continuous [4]. The method for solving Equation (2) does not depend on the form of $\eta \left(x\right)$.

In this paper, to demonstrate the proposed numerical methods, we propose to consider a class of acoustic problems represented by Equation (2). It is also worth noting that other classes of problems of mathematical physics can be described using integral equations [4,22,26,27].

Collocation method on a uniform grid. In practice, in order to solve integral equations, it is necessary to resort to discretization of the solution domain under consideration. To approximate the integral Equation (2), we will use the collocation method [23,28]. For three-dimensional problems, some difficulties arise in the discretization of integral equations defined in regions of complex shape.

Let us represent the area $Q$ as a union of $N_Q$ cells $\mathrm{\Omega }\left(i\right), i=1,\dots , N_Q$. Nodal points in these cells will be chosen in their centers, which are defined by the formulas [29]

$$x_l^c=\frac{{\int }_{\mathrm{\Omega }}{x}_{l}dx}{mes\mathrm{\Omega }},\hspace{1em}\hspace{1em}l=1,\dots ,3,$$

(3)

where $dx=d{x}_{1}d{x}_{2}d{x}_3$ represents the integration of the volumetric $l$-th partition of the solution domain $Q$, $x^c=(x_1^c,x_2^c,x_3^c)$ is the cell center $\mathrm{\Omega }$, and $mes\mathrm{\Omega }$ its volume. If in region $\mathrm{\Omega }$ a differentiable function of its arguments is defined $f\left(x\right)$, then approximate equality is true:

$${\int }_{\mathrm{\Omega }}f\left(x\right)dx\approx f\left(x^c\right)mes\mathrm{\Omega }.$$

(4)

Expression (4) will be an exact equality if $f\left(x\right)$ is a linear function of the arguments.

We will approximate the integral Equation (2) by SLE of dimension $\sim N_Q$ with respect to the values of the unknown function at the nodal points of the domain $Q$, located in the centers of the $x^{ci}$ cells $\mathrm{\Omega }\left(i\right), i=1,\dots , N_Q$. We will choose the cell sizes so that the desired function varies weakly within the cell. Then, redefining the corresponding SLE can be represented in the following form [23,30]:

$$\begin{gathered} u\left(i\right)+\sum _{j=1}^{N_Q}A\left(i,j\right)\eta \left(j\right)u\left(j\right)=u^0\left(i\right),\hspace{1em}\hspace{1em}i=1,\dots ,N_Q, \\ A\left(i,j\right)=-{\mathrm{k}}^2{\int }_{\mathrm{\Omega }\left(j\right)}\frac{K\left(x^{ci}-y\right)}{|x^{ci}-y|}dy; u^0\left(i\right)=-{\int }_Q\frac{K\left(x^{ci}-y\right)}{\left|x^{ci}-y\right|}f\left(x^{ci}\right)dy; \\ i,=1,\dots ,N_Q,\hspace{1em}\hspace{1em}u\left(i\right)=u\left(x^{ci}\right),\hspace{1em}\hspace{1em}{u}^0\left(i\right)=u^0\left(x^{ci}\right),\hspace{1em}\hspace{1em}\eta \left(i\right)=\eta \left(x^{ci}\right). \end{gathered}$$

(5)

To calculate the integrals in Equation (5), we can use the approximate Formula (4) or more accurate numerical integration algorithms. Note that, since the nodal points are located in the center of the cells, the accuracy of approximation of the integral operators $\sim h^2$ where $h$ is the maximum cell diameter (we define the cell diameter as the maximum distance between the boundary points). For relatively small values of $N_Q\le 10000,$ we can solve the system of Equation (5) by direct or iterative methods. Below, we will outline efficient algorithms which, using iterative methods, allow us to solve system Equation (5) with much higher dimensionality.

In the kernel of the integral Equations (1) and (2), there is a term depending on the difference in Cartesian coordinates of points x and y. However, this circumstance was not used in any way in the construction of SLE Equation (5). Below, using a uniform grid and discrete Fourier transform algorithms, we construct an efficient numerical method for solving Equation (2).

Consider a complex function $f\left(n\right)$ of discrete argument $n=0, \pm 1, \pm 2,\dots$., and let us assume that $f\left(n\right)$ is a periodic function with period $N$, i.e., $f\left(n\pm N\right)=f\left(n\right)$ for any n. The discrete Fourier transform of the function $f\left(n\right)$ is defined by the well-known formula:

$$F\left[f\right]=f^F\left(k\right)=\sum _{n=0}^{N-1}\exp \left(i\frac{2\pi }{N}kn\right)f\left(n\right); k=\overline{0,N-1},$$

(6)

where, obviously, the Fourier transformant $f^F\left(k\right)$m is also a periodic function with period N.

If we know the Fourier transformant $f^F\left(k\right)$, then we can recover the original function $f\left(n\right)$ using the inverse discrete Fourier transform:

$$F\left[f\right]=f^F\left(k\right)=\sum _{n=0}^{N-1}\exp \left(i\frac{2\pi }{N}kn\right)f\left(n\right); k=\overline{0,N-1}.$$

(7)

In general, the number of arithmetic operations $T_F\left(N\right),$ which is required to compute the discrete Fourier transform without the cost of computing functions of the form $\mathrm{e}\mathrm{x}\mathrm{p}(i2\pi kn/N),$ is estimated by the formula:

$$T_F\left(N\right)\sim N^2.$$

(8)

When using fast discrete Fourier transform algorithms, the number of arithmetic operations required is estimated by the formula [23]:

$$T_{FF}\left(n\right)\sim N\cdot LOG\left(N\right),$$

(9)

where $LOG\left(N\right)$—is the integer logarithm, i.e., the sum of all prime divisors of N. If N is a power of two, then $T_{FF}\left(n\right)\sim N\cdot lo{g}_2\left(N\right)$.

Let $A\left(l\right)$ bea periodic function of a discrete argument with period N. Consider the sums of the following form:

$$v\left(n\right)=\sum _{m=0}^{N-1}A\left(n-m\right)u\left(m\right),n=\overline{0,N-1}.$$

(10)

The sums in Equation (10) arise from multiplication of circulant matrices by a vector. Let us apply the discrete Fourier transform with period N to both parts of (10). It is not difficult to show that:

$$v^F\left(k\right)=A^F\left(k\right)u^F\left(k\right),k=\overline{0,N-1}.$$

(11)

Using Equation (11) and fast discrete Fourier transform (FFT) algorithms, one can efficiently multiply circulant matrices by a vector. However, circulant matrices rarely appear in real-world problems. However, in many problems, in particular those discussed below, one needs to compute sums of the form Equation (10), in which the function $A\left(l\right),-\left(N-1\right)\le l\le (N-1)$ is arbitrary in the specified range. Such sums arise when multiplying Toeplitz matrices by the vector [31,32]. This function $A\left(l\right)$ is defined at the $(2N-1)$ integer point. Let us further define the function $A\left(l\right)$ zero at the point $l=N$ and extend it to all integer values with period 2N. Further, for the function of the discrete argument $u\left(m\right), m=0,\dots ,\left(N-1\right),$ let us define this as zero at the points $m=N,\dots , 2N-1$. Now, consider the sums of the following form:

$$v\left(n\right)=\sum _{m=0}^{2N-1}A\left(n-m\right)u\left(m\right),n=\overline{\mathrm{0,2}N-1}.$$

(12)

It follows from the above that, at $n=0,\dots ,N-1$ function $v\left(n\right)$ from Equation (12) coincides with the values $v\left(n\right)$ from Equation (10). Further, for quick calculation of the sums in Equation (12), we will use the formula:

$$v^F\left(k\right)=A^F\left(k\right)u^F\left(k\right),k=\overline{\mathrm{0,2}N-1}.$$

(13)

In the inverse Fourier transform, only the components of $\nu \left(n\right), n=\overline{0,N-1}$. Thus, it follows from Equation (9) that the number of arithmetic operations to calculate Equation (10) is estimated by the formula:

$$T_A\sim 2NLOG\left(2N\right).$$

(14)

Moreover, it is necessary to store in the computer memory an array with the number of elements equal to:

$$M_A\sim 2N.$$

(15)

Let us proceed to the discretization of the integral Equation (2). In a rectangular Cartesian coordinate system, define a parallelepiped $\mathrm{\Pi }$, within which the region $Q$ is located. The edges of the parallelepiped are parallel to the coordinate axes, and the lengths of the edges are equal to $N_l{h}_l, l=\mathrm{1,2},3,$ where $h_l$ are the grid steps on the Cartesian coordinates. Then, the parallelepiped $\mathrm{\Pi }$ can be represented as a union of cells (elementary parallelepipeds) $\mathrm{\Pi }\left(p\right), p=\left(p_1,p_2,p_3\right), { p}_l=0, \dots ,N_l-1$. Let us define the area $\tilde{Q}$ as a union $N_Q$ of cells whose centers lie inside the area $Q$. The nodal points, in which the values of functions are defined, will be defined in the centers of cells and denoted as $x\left(p\right)$ and the values of functions in these points as $f\left(p\right)$.

The integral Equation (2) will be approximated, similarly to Equation (5), by an SLE of the following form [24]:

$$\begin{gathered} u\left(p\right)+\sum _{y\left(q\right)\in Q}A\left(p-q\right)\eta \left(q\right)u\left(q\right)=u^0\left(p\right),\hspace{1em}\hspace{1em}x\left(p\right)\in Q, \\ A\left(p-q\right)=-{\mathrm{k}}^2{\int }_{\mathrm{\Pi }\left(q\right)}\frac{K\left(x\left(p\right)-y\right)}{\left|x\left(p\right)-y\right|}dy,\hspace{1em}\hspace{1em}p\ne q,\hspace{1em}\hspace{1em}A\left(0\right)=0. \end{gathered}$$

(16)

Since the nodal points are at the center of the cells, the accuracy of the approximation of the integral operator is $\sim h^2, h=\sqrt{h_1^2,h_2^2,h_3^2}$.

It follows from Equation (16) that the main computational cost of multiplying the SLE matrix by a vector (performing one iteration) is associated with the computation of sums of the form:

$$W\left(p\right)=\sum _{y\left(q\right)\in Q}A\left(p-q\right)V\left(q\right),x\left(p\right)\in Q.$$

(17)

To calculate $W\left(p\right)$ at the nodal points $x\left(p\right)\in Q$ requires performing $\sim N_Q^2$ arithmetic operations, where $N_Q$—is the number of nodal points in the domain $Q$. To reduce the number of arithmetic operations, we will apply the technique of fast multiplication of Toeplitz matrices by vector, as described above.

Let us define the function $V\left(q\right)$ as zero at the points $x\left(q\right)$ of the parallelepiped P, not belonging to the area $Q$. Consider the following sums:

$$W\left(p_1,p_2,p_3\right)=\sum _{q_1=0}^{N_1-1}\sum _{q_2=0}^{N_2-1}\sum _{q_3=0}^{N_3-1}A\left(p_1-q_1,p_2-q_2,p_3-q_3\right)V\left(q_1,q_2,q_3\right).$$

(18)

It is obvious that at $x\left(p\right)\in Q$ values $W\left(p\right)$ from Equations (17) and (18) coincide. In Equation (18) the matrix function of the discrete argument $A\left(p\right)$ is defined for the values $-\left(N_1-1\right)\le p_1\le \left(N_1-1\right);-\left(N_2-1\right)\le p_2\le \left(N_2-1\right);-\left(N_3-1\right)\le p_3\le \left(N_3-1\right).$

Let us denote by ${\mathrm{\Pi }}_2$ a parallelepiped with sides $2N_1{h}_1, 2N_2{h}_2$ ᴎ $2N_3{h}_3.$ Let us continue the matrix function of the discrete argument $A(p_1,p_2,p_3)$ to all integer values $p_1,p_2,p_3$ assuming it to be periodic for each variable, with periods, respectively, $2N_1,2N_2,2N_3.$ At the same time, let us define the function $A(p_1,p_2,p_3)$ as zero at the points where it is not defined. Next, let us define the function of the discrete argument $V(p_1,p_2,p_3)$ as zero at all nodal points ${\mathrm{\Pi }}_2$, not belonging to $\mathrm{\Pi },$ and extend this to all integer values $p_1,p_2,p_3$, assuming it is periodic for each variable, with periods, respectively, $2N_1,2N_2,2N_3$.

Consider the expression:

$$W\left(p_1,p_2,p_3\right)=\sum _{q_1=0}^{2N_1-1}\sum _{q_2=0}^{2N_2-1}\sum _{q_3=0}^{2N_3-1}A\left(p_1-q_1,p_2-q_2,p_3-q_3\right)V\left(q_1,q_2,q_3\right).$$

(19)

Considering the above, it is clear that at $x\left(p\right)\in Q$ function $W(p_1,p_2,p_3)$ from Equation (19) coincides with the values from Equation (17). Below, by $\mathrm{\Pi }$ ᴎ ${\mathrm{\Pi }}_2$ we denote integer parallelepipeds with the number of discrete arguments on each axis $N_1,N_2,N_3$ ᴎ $2N_1,2N_2,2N_3$, respectively. Now, performing a discrete Fourier transform on each variable from both parts of Equation (19), we obtain the following equality:

$$W^F\left(k_1,k_2,k_3\right)=A^F\left(k_1,k_2,k_3\right)V^F\left(k_1,k_2,k_3\right),\hspace{1em}\hspace{1em}k\in {\mathrm{\Pi }}_2.$$

(20)

Thus, to perform one iteration when solving the SLE Equation (16), it is necessary to perform the direct Fourier transform of the function $V(p_1,p_2,p_3)$ for each variable and the inverse transform of the function $W^F(k_1,k_2,k_3)$ (the transformation of the function $A(p_1,p_2,p_3)$ is performed once before starting the iteration procedure). The number of arithmetic operations and the amount of memory required to perform one iteration are estimated by the formulas:

$$T_A\sim 10N LOG\left(N\right),\hspace{1em}\hspace{1em}{M}_A\sim 10N,\hspace{1em}\hspace{1em}N=N_1{N}_2{N}_3.$$

(21)

When choosing grid steps and values $N_1,N_2,N_3$, it is necessary to be guided by the following criteria: first, the desired function does not change much within the cells; second, the region $\tilde{Q},$ consisting of cells whose centers are inside $Q,$ describe $Q$ well enough.

Performance indicators of numerical algorithms. Having mentioned earlier that the solution of integral equations by numerical methods reduces to the solution of SLE with fully filled matrices, let us explain the main criteria for the efficiency of the algorithms, and demonstrate the visual advantage of iterative methods over direct methods in these classes of problem.

The main efficiency parameters of numerical algorithms are the number of operations T required to solve the initial problem, the amount of memory M required to implement the algorithm, the storage of $N^2$ elements of the SLE matrix with $N$ unknowns, etc. It is obvious that a large number of computational resources is required to solve the problem under consideration. In the case of the iterative method, these properties of the algorithm can be estimated using the relations:

$$T\sim L{T}_A,\hspace{1em}\hspace{1em}M\sim M_A,$$

(22)

where $T_A$ is the number of arithmetic operations required to multiply the SLE matrix by a vector, which is estimated in terms of complexity as $T_A\sim N^2$; L is the number of iterations required to obtain a solution with a given accuracy; and $M_A$ is the number of unique elements of the matrix.

Let us denote by $T_G, M_G$ the characteristics of the direct methods for solving SLE. In comparison to this, when using the direct Gauss method to solve SLE, it is necessary to perform $T_G\sim N^3$ arithmetic operations and store in the computer memory at the order of $M_G\sim N^2$ numbers, which demonstrates the preference for using iterative methods in solving these problems.

Note that only iterative methods can be used to solve SLE Equation (16) using the considered algorithm. This is due to the fact that iterative algorithms are based on multiplications of the SLE matrix by a vector. The number of arithmetic operations and the amount of memory required to solve SLE Equation (16) are estimated by Formulas (21), (22). At the same time, the number of iterations required to obtain the solution is usually much smaller than the dimensionality of the SLE. Thus, it is possible to numerically solve the integral Equation (1) and, in the special case, Equation (2), which is reduced to the SLE of dimensionality $N_Q\ge {10}^6$.

Iterative method for solving the problem. In the numerical formulation of Equation (5), and in particular Equation (16) for the case with approximation of the volumetric domain of the solution $Q$ by a system of parallelepipeds $\mathrm{\Pi }\left(q\right), q\in Q,$ we will solve the initial problem Equations (1) and (2) by the iterative method of gradient descent [33]. Let us describe the iterations and peculiarities of the method for solving a particular problem.

The task of finding a solution to an SLE is to solve the operator equation:

$$\left(H\right)u=b,$$

(23)

where $\left(H\right)$ is a known, in the general case a complex, matrix operator; $b$ is a known, in the general case complex, vector of the right-hand side; and $u$ is the vector whose values we want to find. When solving the problem Equations (1) and (2) by numerical methods, the values of the elements of the matrix operator $H$ are found as a result of discretization of Equations (5) and (16) of the initial integral equation and calculation of integrals over the obtained partitions $\mathrm{\Pi }\left(q\right)$.

In addition, the action of the integral operator Equation (1) or Equation (2) with respect to the unknown vector $u$ is not the usual multiplication of a high dimensional matrix by a vector. Strictly speaking, it follows from Equation (2) that:

$$\begin{gathered} \left(H\right)u=u+A\left(\widehat{\eta }u\right), \\ \left(H\right)u=u\left(p\right)+\sum _{y\left(q\right)\in Q}A\left(p-q\right)\widehat{\eta }\left(q\right)u\left(q\right),p\in Q, \end{gathered}$$

(24)

where $A\in {\mathbb{C}}^{N_Q\times N_Q}$ is the matrix of coefficients of the integral equation kernel obtained as a result of calculation of Equations (4) and (5); $\widehat{\eta }\in {\mathbb{C}}^{N_Q}$ is the vector of refraction values at each point of the $p$ of the given area partition $p\in Q$; $u\in {\mathbb{C}}^{N_Q}$, also in the general case, is the complex vector of unknown values of the scalar field strength at each discrete point of the given region. Each matrix-to-vector multiplication of any iterative algorithm that will solve the given problem of Equation (2) must be further defined according to Equation (24), resulting in iterated methods.

In the iterations of the modified iterative gradient descent method, the action of the conjugate operator is also presented $\left(H^{*}\right)$ on the vector $u$. Considering that, by definition for the matrix operator $H^{*}={\overline{H}}^T$ where $\overline{H}$ is a matrix with complex-conjugate elements to the elements of matrix $H$, then the action of the operator Equation (24) can be rewritten in the form:

$$\begin{gathered} \left(H^{\ast }\right)u=u-\overline{\eta }\left({\overline{A}}^{T}u\right), \\ \left(H^{\ast }\right)u=u\left(p\right)+\overline{\eta }\left(p\right)\sum _{y\left(q\right)\in Q}\overline{A}\left(q-p\right)u\left(q\right),p\in Q \end{gathered}$$

(25)

where $\overline{\eta }$ is also a refraction vector with complex-conjugate elements to vector $\widehat{\eta }.$ The iterative method aims to find an approximation $u_m\approx u$ of the unknown desired function. Iterations of the modified iterative method of gradient descent are defined as follows:

$$r_0=\left(H\right)u_0-b;u_1=u_0-\frac{{‖\left(H^{\ast }\right)r_0‖}^2}{{‖\left(H\right)\left(H^{\ast }\right)r_0‖}^2}(H^{\ast })r_0;$$

(26)

$$\begin{gathered} r_m=\left(H\right)u_m-b;\hspace{1em}\mathsf{\Delta }{r}_m=r_m-r_{m-1}; \\ {t}_m{‖\mathsf{\Delta }{r}_m‖}^2+h_{m}Re\left(\mathsf{\Delta }{r}_m,\left(H\right)\left(H^{\ast }\right)r_m\right)=Re\left(r_m,\mathsf{\Delta }{r}_m\right), \\ {t}_{m}Re\left(\mathsf{\Delta }{r}_m,\left(H\right)\left(H^{\ast }\right)r_m\right)+h_m{‖{\left(H\right)(H}^{\ast })r_m‖}^2=Re\left(r_m,\left(H\right)\left(H^{\ast }\right)r_m\right), \\ {u}_{m+1}=u_m-t_m\left(u_m-u_{m-1}\right)-h_m\left(H^{\ast }\right)r_m,\hspace{1em}\hspace{1em}m=\mathrm{1,2},\dots . \end{gathered}$$

(27)

The only restriction on the iterative method for Equations (24) and (27) is the existence of a bounded inverse operator to $\left(H\right)$. Proofs of convergence of iterations for Equation (26) and convergence analysis of iterations and increasing dimension of the matrix are presented in [33].

As a criterion for stopping iterations for Equations (24) and (27), we choose the metric ${\delta }_m$ as the relative error of the approximated vector at step $m$ iterations:

$${\delta }_m=\frac{‖u_m-u_{m-1}‖}{‖b‖}<\epsilon ,$$

(28)

where $u_m \mathsf{ᴎ} u_{m-1}$ are the obtained approximations of the unknown scalar field $u$ at the partition points; and $\epsilon$ is the given accuracy of iterations, which is most often set as equal to ${10}^{-d+1},$ where $d$—is the number of significant digits of the computer representation of floating-point numbers.

As a result, we have a method that effectively copes with the problem of solving operator equations. The study and comparison of the iterative method for the real problem of solving SLE with a fully filled matrix was carried out in our previous work [25].

3. Results

Let us present the formulation of the conducted numerical experiment devoted to the solution of the acoustics problem Equation (2).

The rectangular solution domain $Q$ is characterized by linear dimensions $l=1$ on each of the Cartesian axes $x=(x_1,x_2,x_3)$ and the point of the center of the cubic region $c=(0, 0, 0)$ at the origin of coordinates. The propagation of a plane wave in a medium is characterized by the value of the wave number $k=15$ as well as the vector of the wave propagation direction $\overrightarrow{v}=(1/\sqrt{3} , 1/\sqrt{3} , 1/\sqrt{3})$. The external function $f\left(x\right)$ from Equation (2) is modeled as complex harmonic oscillations:

$$f\left(x\right)=\exp \left(-ik\left(x,\overrightarrow{v}\right)\right),\hspace{1em}\hspace{1em}x\in Q,$$

(29)

with a given value of the wave number $k$ and vector of propagation direction $\overrightarrow{v}$.

Let us set for the time of the experiment the following function $\eta \left(x\right)$ of refraction of the transparent bulk medium:

$$\eta \left(x\right)=\left\{\begin{array}{c}2.5+1.5i={\eta }_1\left(x\right),\hspace{1em}\hspace{1em}where x_1\in (-l/4,l/4),\\ 1.0+0.0i={\eta }_2\left(x\right),\hspace{1em}\hspace{1em}where x_1\in \left[-\frac{l}{2},-\frac{l}{4}\right]\cup \left[\frac{l}{4},\frac{l}{2}\right].\end{array}\right.$$

(30)

In fact, this formulation of the refraction function does not correspond in any way to a concrete real object, such as a transparent screen with appropriate physical properties. This value of refraction was chosen on the basis of considerations of conducting a numerical experiment with two media having different physical properties, in order to obtain the effect of refraction of a wave and attenuation of its energy when passing through a point $l/4$. A flat screen is defined by a stepwise refraction transition along the Cartesian coordinate $x_1$. As a result of Equation (30), we have a piecewise constant given function $\eta \left(x\right)$ of refraction of the region $Q,$ which defines a transition from the medium with constant real refraction ${\eta }_2\left(x\right)$ to a constant complex refraction ${\eta }_1\left(x\right)$.

Discretization of the area $Q$ is according to Equation $\left(16\right),$ into parallelepipeds of equal volume with equal sides. Modeling the problem in the cubic space, we will vary $N$, the number of such cubes, which are located along each of the axes of the Cartesian coordinate system, and the number of which along each of the coordinates is the same. Then the total number of partitions is strictly equal to $N_Q=N^3$. Within the results of this paper, we will compare the influence of the degree of discretization of the problem on the final quality of the solution of the problem for Equation (2).

In a series of numerical experiments, we computed (16) based on the modified iterative gradient descent method for Equations (24)–(27), in the private formulation of Equations (29)–(30), with specifically specified parameters. In each of the selected discretizations, for convenience, we set the stopping criterion of the iterative method $\epsilon ={10}^{-4}$ for the metric Equation (28).

The obtained solution characterizes the unknown value of the scalar field of intensity $u$ at each point of the center of partitioning of the initial region $Q$ into parallelepipeds of equal volume. Let us visualize the obtained results with sufficiently small $N=25$ in Figure 1. The total number of obtainable elementary partitions in such a case is $N_Q=15626$. The figure shows the obtained solution $u$ as a result of convergence of iterations for Equations (24)–(27). The first sub-drawing in Figure 1a shows the values of the $\left|u\right|$ of the potential field modulus in the center of each cube, which vary in the range of $\left|u\right|\in (0.0025, 0.0325)$. The gradation of these values is also highlighted in the graph. The sub-drawing Figure 1b shows the real part of the complex scalar field $u$ whose values vary in the range $Re u\in \left[-0.025, 0.025\right].$

According to the character of the obtained picture, we obtain the expected type of the solution. In Figure 1b, one can notice the refraction of the plane wave Equation (29) incident from the point of $x=(0.5,0.5,0.5)$ against the boundary located at the point along the axis $x_1=l/4=0.125.$ In the limit $-l/4<x_1<l/4,$ how the refracting wave changed the direction of its propagation can be seen, and as it passed within this region it experienced attenuation, due to which both the modulus of its value and the amplitude of its real values decreased. After passing through the point $x_1=-l/4=-0.125,$ the propagation resumes due to the run-up of the plane wave front, as in Equation (29).

We also visualize the obtained results with a sufficiently large $\mathrm{N}=100$ in Figure 2. The limits of the scalar field values remained the same, and also the obtained structure for the solution remained unchanged. The total number of obtainable elementary partitions in this case amounted to ${\mathrm{N}}_{\mathrm{Q}}={10}^6$.

As a result, having increased the linear partitioning of the domain by approximately a factor of 4, we obtain a 64-times increase in the discretization of the domain $N_Q,$ in view of which the obtained solutions show more features at the refraction boundary. Let us note the smoothness of the obtained solution, as a result of which it is possible to take into account small features both on the refraction boundary and on possible inclusions, whose sizes do not exceed the linear sizes of the obtained partitions $\mathrm{\Pi }\left(q\right)$.

On the center slices on the axis $x_3\approx 0$ areas $Q$ in Figure 3, we also see the picture of the refraction and attenuation of the wave, both in the modulus value, Figure 3a, and in the values of amplitudes of the real Figure 3b and the complex part of the wave Figure 3c. On the graph of the modulus approximation of the scalar field values $\left|u_m\right|,$ Figure 3a, we can clearly see the transition in the range of $-l/4<x_1<l/4$, which is a consequence of the transition of the refraction value in this region. In Figure 3b,c, we see refraction in the pattern of propagation of the plane wave, as well as the attenuation of the amplitude of the scalar field values within $-l/4<x_1<l/4$.

Performing calculations with different degrees of discretization of the solution domain, we measure the efficiency of the presented iterative method, Equations (24)–(27), using values of different metrics. For solutions with different discretization of the domain $Q,$ we measured the total number of matrix multiplications by the vector $m$ during the algorithm operation, the relative error of approximations of iterations ${\delta }_m$ counted by (28), the norm of incoherence of the obtained approximate solution $\left|\right|\left(H\right)u_m-b\left|\right|$, and the maximum value of the modulus of the scalar field $\underset{p\in Q}{\max }\left|u_m\right(p\left)\right|$. In addition to the above metrics, we measured the program execution time in seconds ($t, sec$) for a given discretization ($N$) of the volume area of the solution on a personal computing device with an intel core i5-10400f processor and DDR4 RAM, with a frequency of 3200 MHz and a volume of 32 GB. Time was measured via the implementation of a C++ program in a sequential single-threaded execution of instructions.

Table 1. Table of performance and convergence quality of the method with different numbers of partitions.

$\mathbf{N}$	${\mathbf{N}}_{\mathbf{Q}}$	m	${\mathsf{\delta }}_{\mathbf{m}}$	$\left|\right|\left(\mathbf{H}\right){\mathbf{u}}_{\mathbf{m}}-\mathbf{b}\left|\right|$	$\underset{\mathbf{p}\in \mathbf{Q}}{\mathbf{max}}\left|{\mathbf{u}}_{\mathbf{m}}\right(\mathbf{p}\left)\right|$	$\mathbf{t},\mathbf{s}$
10	1000	48	$9.652\times {10}^{-5}$	$4.221\times {10}^{-5}$	$5.618\times {10}^{-6}$	0.833
15	3375	48	$3.587\times {10}^{-5}$	$3.854\times {10}^{-5}$	$3.198\times {10}^{-6}$	2.895
20	8000	42	$5.366\times {10}^{-5}$	$2.984\times {10}^{-4}$	$1.491\times {10}^{-5}$	6.131
25	15,625	45	$4.784\times {10}^{-5}$	$2.601\times {10}^{-4}$	$9.216\times {10}^{-6}$	13.065
30	27,000	57	$7.434\times {10}^{-5}$	$4.631\times {10}^{-4}$	$1.100\times {10}^{-5}$	29.069
35	42,875	45	$8.941\times {10}^{-5}$	$2.603\times {10}^{-4}$	$6.935\times {10}^{-6}$	36.775
40	64,000	42	$6.731\times {10}^{-5}$	$3.751\times {10}^{-4}$	$6.805\times {10}^{-6}$	51.703
45	91,125	42	$5.311\times {10}^{-5}$	$3.083\times {10}^{-4}$	$5.210\times {10}^{-6}$	74.258
50	125,000	54	$2.992\times {10}^{-5}$	$5.041\times {10}^{-4}$	$6.793\times {10}^{-6}$	132.335
100	1,000,000	42	$9.680\times {10}^{-5}$	$1.273\times {10}^{-4}$	$5.321\times {10}^{-6}$	1039.317

demonstrates the values of the above metrics for different numbers of partitions $N_{Q.}$ source region $Q$.

According to the results, we note that there is no dependence of the investigated values of the metrics on the number of partitions of the solution domain. It is shown that the number of matrix multiplications by vector $m$ necessary to obtain an approximation $u_m$ with a predetermined accuracy $\epsilon$ does not increase depending on the number of partitions of the domain $Q$. However, we observe an oscillation of this value, which may be due to a poor choice of cube sizes, which resulted in an inaccurate description of the refraction region.

Metrics ${\delta }_m$, $\left|\right|\left(H\right)u_m-b\left|\right|$ do not change depending on the number of partitions; we can only note the low value of the norm of non-convexity of the integral operator with respect to the obtained solution at small discretizations $N_Q=1000$ ᴎ $N_Q=3375$. This circumstance can also be related to the number of final calculations or to a weak description of the solution domain of the problem. The index $\underset{p\in Q}{\max }\left|u_m\right(p\left)\right|$ also varies within the error limits.

We display the graphs of the metric Equation (28) of the relative error ${\delta }_m$ of approximations during the iterative method for Equations (24)–(27) depending on the number of applications of the matrix-to-vector multiplication operation $m$. Figure 4 shows the value of ${\log }_{10}{\delta }_m$ depending on $m,$ due to the exponential nature of the reduction of this metric. By prologarithmizing ${\delta }_m,$ we see the convergence pattern of the method in Equations (24)–(27) for the problem in Equations (29)–(30), in the values of powers of 10 of the base of the number system.

In the linearized dependence shown, we see by the convergence pattern in Figure 4 an oscillating graph of metric Equation (28), demonstratively converging (in natural values) to the required value of approximation accuracy $\epsilon$. It can be seen from the graph that, with increasing discretization (n on the graph is equal to $N$ or $N_1$, $N_2$ or $N_3$) of the initial domain, the number of iterations (matrix-to-vector multiplications) does not increase, which cardinally contributes to the speed in solving the initial formulation of the integral equation at large discretizations.

4. Discussion

The solution of problems using the most general approaches to approximation of complex regions is widely used in various works on mathematical modeling of wave effects on physical bulk multilayer objects of different shapes with different boundary conditions. To date, there are many ways to solve wave propagation and scattering problems in volumetric media, in particular for the problem of acoustics, in which, often, complex discretizations based on tetrahedral meshes are applied. The approach based on tetrahedra as units of elementary partitioning of the three-dimensional solution domain is due to the need for the most accurate description of complex domain boundaries to the detriment of the possibility of fast computations of problems on uniform mesh approximations.

The computational approach demonstrated in this paper will be useful in the case of a large number of different refractive inclusions that can be described by rectangular parallelepipeds. This method will also be useful in modeling subtle wave processes due to the possibility of performing efficient computations over fully filled matrices of large size, $N_Q\ge {10}^6$. The condition of applicability of this computational method is the requirement of the presence of the kernel of the integral operator depending on the coordinate difference; therefore, the main computations are associated with the computation of the coefficients of the operator in the form of the Teuplice matrix. The numerical experiment presented demonstrates an effective approach to solving such a problem formulation, as well as qualitative pictures for solving the problems of acoustic wave propagation in a medium with attenuation and wave scattering at the refraction boundary of transparent bulk media.

It is noteworthy that there are no requirements for the structure of the refraction function, except for the physical limitations of the modeled materials, which makes it possible to perform mathematical modeling of the problems presented above for areas of high complexity. The main limitation will be the size of the area on which the function will be set.

The modified iterative method of gradient descent proposed for solving the problem has shown itself to be an effective method for solving the operator equation with a fully filled operator matrix. According to the results of numerical experiments, it is shown that it converges steadily when the dimensionality of the problem changes and also converges steadily for a reasonable number of matrix-to-vector multiplications (about 100 operations).

Moreover, when the dimensionality of the problem increases, the number of such operations does not change, which is associated only with the physical properties of the problem formulation. The obtained numerical scheme allows us to seriously accelerate computations in general formulations of the solution of SLE or operator equations and, with the use of special discretization, will also allow us to solve a large number of unknowns in applied modeling problems.

It is also worth noting that, due to the currently prevailing most accurate approach to describing a three-dimensional domain based on complex elementary partitions, comparative results with other studies cannot be given.

5. Conclusions

Today, the problem of economical fast computations for mathematical models of physical processes to provide high quality and fast computations over digital doubles of real in-situ experiments is acute. The development of efficient numerical methods for solving a narrow class of problems in mathematical physics makes it possible to solve such problems with increasing requirements on the degree of discretization due to the reduced complexity of calculations and the required memory for the implementation of the algorithm. This can be achieved by using the peculiarities of the kernels of the initial formulation in the form of solving Fredholm integral equations describing mathematical models of some real processes.

This paper shows how, taking into account the features of the integral equation kernel function and introducing a special spatial grid, it is possible to obtain a multiple acceleration of calculations and a reduction in the required memory for the realization of calculations.

The paper proposed a formulation for solving problems in mathematical physics using integral equations and iterative solution methods. The theoretical efficiency of application of fast Fourier transform in the optimization of multiplication of the Toeplitz matrix by a vector is shown. Numerical results are demonstrated on the model problem of acoustic wave propagation in a transparent bulk medium with an inhomogeneous refraction index.

In future works, it is assumed that the use of this numerical method for discretizing the volume domain of the solution, the technique of multiplying a Toeplitz matrix by a vector for optimization of iterative methods according to a given discretization form, as well as they setting up of the solution of the integral equation problem, will be applied in related tasks in the propagation and scattering of waves in electromagnetism, which will also find its application in applied areas such as radio-physics, remote sensing, radio-spectroscopy, and crystal growth.

Continuing this work will also be considered towards optimizing the numerical method of solving linear systems with Toeplitz-type matrices, both in terms of specifically optimizing multiplication with the use of multilevel Toeplitz matrices and in terms of optimizing the iterative method proposed. This direction of work will be associated with attempts to simplify the use of the proposed methodology for solving volumetric integral equations in various related tasks, as well as to speed up computations and reduce estimates of the required memory for iterations.