In the paper three-dimensional stationary scalar difraction problems are considered. They are formulated in the form of boundary weakly singular Fredholm integral equations of the first kind with a single function, each of which is equivalent to the original problem. These equations are approximated by systems of linear algebraic equations, which are then solved numerically by an iterative method. In order to reduce the computational complexity, the mosaic-skeleton method is used at the stage of numerical solution of the systems.
coefficients of systems of linear algebraic equations (SLAE), approximating corresponding integral equations, by very simple formulas. Once the approximate solution of integral equation is found, the required approximate solution of the diffraction problem by using integral representations is restored in any point of space.
Matrices of the obtained SLAE are dense. The complexity of solving such SLAE by direct methods is ( 3), where is the order of the system. It was established earlier [1] that the use of the generalized minimal residual method (GMRES) [6] allows to lower the complexity to ( 2), where is the number of iterations, which is much less than and barely depends on it. The most timeconsuming part of the algorithm of GMRES is the use of multiple matrix-vector multiplication. The time expenses on multiplication can be lowered using the ”quick” method, the complexity of which is ( 2), when → ∞. The mosaicskeleton method [7] – [11] can be used for already developed software for solving problems of mathematical physics. Its basic idea is to approximate the large blocks of dense matrices by low-rank matrices. As a result, the SLAE matrix is not fully computed and stored, but its approximation is used in the procedure of matrix-vector multiplication. In this case, the complexity of multiplication of this approximation by a vector is almost linear.
To implement the method in an already established program for numerical solution of diffraction problems, it is only necessary to add the procedure of creating and storing the approximate matrix and to change the module of matrix-vector multiplication, whereas the method of sampling and the calculation procedure of the elements of the original matrix remain the same. An additional advantage of this method is that, due to the cost effective storage of the approximate matrix, the requirements for computer resources are reduced.
The present work outlines the numerical method for solving three-dimensional problems of diffraction formulated as boundary Fredholm integral equations of the first kind and describes the use of the mosaic-skeleton method for the numerical solution of such problems. It also presents the results of numerical experiments, which allow us to judge the effectiveness of the approach. 2
Let us consider a three-dimensional Euclidean space R3 with orthogonal coordinate system 1 2 3, filled with a homogeneous isotropic medium with density , with speed of propagation of acoustic oscillations and absorption coefficient , which has a homogeneous isotropic inclusion, limited by the arbitrary closed surface , with density , speed of sound and absorption coefficient . Let’s denote domains R3, occupied by the inclusion and containing medium, 3 ¯ ). by and ( = R ∖
Suppose there are harmonic sound sources in space , stimulating initial pressure wave field 0 in the containing medium. Sound waves come through space and scatter, reaching the inclusion. As a result, there appear reflected waves in domain , and transmitted waves in domain . Therefore, the complex amplitude of the complete field of pressures can be presented as: = ︂{ , + 0, ∈ , ∈ , reflected wave fields. shown in the Figure 1. where , are complex amplitudes of the pressure field of transmitted and The inclusion, initial, reflected and transmitted pressure wave fields are conjugation conditions on the material interface between and as well as the radiation condition at infinity ︀⟨ − − +, ⟩︀ = ⟨ 0, ⟩
∀ ∈ −1/2( ), ︀⟨ , − − + ⟩︀ = ⟨, 1⟩
∀ ∈ 1/2( ), / | | − = | |−1)︁ , ︁( | | → ∞, if functions 0 ∈ 1/2( ) and 1 ∈ −1/2( ) are set on the boundary . identities
︁∫ ( )
︁∫ ( )
Let’s formulate the initial problem.
Problem 1. In bounded domain of three-dimensional Euclidean space R
3 and
unbounded domain = R ∖
1, find complex-valued functions ( ) ∈ 1( ( ), ), which satisfy integral
3 ¯ , separated by closed surface
∈ + , + >
∇ ( )∇ *( )
− 2( )
( )
*
( )
= 0
∀ ( ) ∈ 01 (︀ ( )︀) ,
(
In works [1], [2] the next theorem was proven.
Theorem 1. Problem 1 has no more than one solution.
Let’s introduce the following notation
︀( ( ) )︀ ( ) ≡ ︀⟨ ( )(, ·), ⟩︀ , (︀ ( ) )︀ ( ) ≡ ︀⟨ ( )(, ·), ⟩︀ , ︁( ( ) * ︁)
( ) ≡ ︀⟨ (·) ( )(, ·), ⟩︀ , ( ) (, ) = exp (︀ ( ) | − |)︀⧸ (4 | − |).
The solution to problem 1 will be sought in the form of potentials
(
+ 1)︀ − * (︀ + + 0)︀ ( ), ∈ ,
∈ −1/2( ) is an unknown density, 0 ∈ 1/2( ), 1 ∈ −1/2( ),
The kernels of these integral operators are fundamental solutions of the
Helmholtz equations and their normal derivatives. Therefore, as shown in [2],
they satisfy identities (
⟩ = ⟨ 2, ⟩ ∀ ∈ −1/2( ),
Problem 1 allows another equivalent formulation in the form of Fredholm integral equation of the first kind with a weak singularity in the kernel [1]. We shall seek its solution in the form
( ) = ( ) ( ), ∈ ,
Theorem 2. Let 0 ∈ 1/2( ), 1 ∈ −1/2( ), > 0 or is not the eigen
+ 2 = 0, ∈ ,
− = 0.
−1/2( ) and formulae (
The applied method of numerical solution is presented in work [3] and appears to be the development of the technics proposed and tested for the first time in [4]. Let us briefly describe the general scheme of its implementation. Let us construct a surface coating by system { } =1 of neighborhoods of nodes ′ ∈ , lying within the spheres of radii ℎ centered at ′ , and denote its subordinate partition of unity by { }. In place of let us take functions ( ) = ′ ( ) ′ ( )
′ ( ) = ︃( ︁∑ =1 ︃) −1 , ︂{ (︀ 1 − 2 ⧸︀ ℎ 2 )︀ 3 , < ℎ , ≥ ℎ , where ∈ , = | − ′ |, ∈ 1( ) when ∈ + , + > 1.
Approximate solutions of equations (
Instead of unknown function set on we shall seek generalized function , acting according to the rule ( , )R3 = ⟨, ⟩ ∀ ∈ 1(R3).
We shall approximate this function by expression ( ) ( ) ≈
¯ ( ), ︁∑ =1 where are unknown coefficients, equality
In paper [3] it is shown that for any functions ∈ 1(R3) and ∈ 1( ) ⟨, ⟩ =
¯ ( , )R3 + (ℎ 2).
Approximation of the density of a single layer potential by a volume density
allows us to obtain simple formulae for approximating integral operator ( )
from (
The integral operators from (
√2 (︂ ¯ , ± = + 2 − ± , 2 3 )︂ )︃ 3 ,
( ) = ( ) = 8 ¯ ¯ exp (︀ 2 4 ¯2 exp (︀ 2 ︀) ︃( ( ) ≡
︀( ( )︀) , (
) + 2 = 2 + 2 ,
= 0.5 ( ) = − √
/ , 2 = −1, exp (︀ − 2)︀ ∫︁ ∞
exp (︀ 2)︀ , =1 ︁∑ =1 ︁∑ =
2
︁∑
=1
︀⟨ ( ), ≈
︀⟩
( ) , = 1, 2, ..., ,
(
4
2 ( ) = (− | | + + Gs ) ¯ , point .
⟨,
⟩ ≈
⟨
︁∑
=1
⟩
≈
︁∑
=1
+ *( ),
( ) ,
= 1, 2, ..., ,
(
− 1︀) ¯ ¯ , ̸= ,
=
︁∑
3
are components of the unit vector of the external normal to the surface at
The operators on the left sides of equations (
,
⟩ ≈ −
,
= 1, 2, ..., ,
(
Solving the corresponding SLAE, we find the approximate values of the density of integral equations at discretization points. After that, an approximate solution of the diffraction problem using integral representations can be calculated at any point in space. 4
resulting from by adding zeros up to .
Let us introduce the definitions necessary for describing the mosaic-skeleton method [7]. Let be block matrix × , and ( ) be matrix of size × ,
The systems of blocks { } will be called covering , if and mosaic partitioning , if, in addition, ⋂︀ = ∅.
The mosaic rank of matrix covering, is number ︁∑ = ( ),
︁∑ mr =
mem /( + ),
min{ , ( + ) rank }.
where the sum is taken over all blocks of covering ∈ C × , and mem =
∈ C × , corresponding to some
(
An important indicator of the effectiveness of this method is the compression factor [7]. It is defined by the following formulae where mem is the amount of memory required to store a matrix in a low-rank format, mem is the amount of memory for storing the original matrix. Definition 3. We shall call a matrix of * kind a skeleton, where and are column vectors, * denotes the row-vector Hermite-conjugate to vector .
From definition 3 it follows that rank ( *) = 1.
The mosaic-skeleton method consists of three phases: constructing a cluster tree, creating a list of blocks and finding a low-rank approximation. For the first phase, the input data is a grid, on which discretization is carried out. The set of all points of the grid is called a zero-level cluster (or a root of the tree). At each step, the original cluster is split into several disjoint subclusters. This continues until the level of the tree reaches the prescribed one. As a zero-level cluster, we can take a cube, and plunge the domain on the border of which the grid is built there. When moving to the next level of the tree, each edge of this cluster bisects, resulting in 8 subcubes. Then, all the points belonging to the original cube are distributed across 8 subcubes, forming 8 subclusters (some of them might be empty). On reaching the maximum level of , the number of clusters equals 8 .
The next step is creating the list of blocks. Any two clusters determine a
block in a matrix. If the points of the clusters are geometrically separated from
each other, the block enters a far zone and an approximation will be constructed
for it, or, otherwise, it enters a near zone, and the elements of the blocks will be
calculated by formulae (
At the final stage, the blocks of the far zone are approximated by low-rank matrices. Such matrices can be sought in different ways [10, 11]. In paper [13] the cross of the matrix row and column is chosen as a template. A current approximation is built on the cross at each step. This algorithm is called incomplete cross approximation [8, 14].
Algorithm 1. Incomplete cross approximation.
Approximation of matrix of size × by matrix ( ), being the sum of skeletons. 1. Let be the number of the calculated skeleton. At this step, we assume that = 1 and select an arbitrary number from the column of matrix . 2. In the column with number we calculate all the elements of matrix , then, subtract the elements of all previously obtained skeletons in these positions from them. In the obtained column we find the maximum element in modulus. Let it be located in the row with number . will be denoted by +1.
of matrix 3. We calculate the row of matrix
with number and subtract from it the elements of all already found skeletons in their respective positions. In the obtained row * we find the maximum element in modulus. It should be noted that an element from the column with number cannot be selected again. The number of the column in which the maximum element is found 4. Along the cross with centre ( , ) we build a skeleton so that the coefficients ( ) = ︁∑
=1 * coincide with the coefficients of the original matrix in positions of computed columns 1, ..., and of calculated rows 1, ..., . 5. We check the stopping criterion. If it is satisfied, the calculation stops. If it is not satisfied, we set := + 1 and repeat the algorithm from step 2.
An approximation (
(17)
Using recurrence formulae we can obtain next result
||
( )||2 = || ( −1)||2 + 2Re︁( ∑︁
−1 ( * )( * ) ︁)
=1
*
+ || ||2 || ||2 ,
| |
2
||
(
| 1 1 | 5
The program for the numerical solution of diffraction problems is written in
Fortran 90 and has the form of a console application, designed to run on
multiprocessor computing systems. Intel Fortran Compiler performs as a compiler,
(
Example 1. The diffraction problem of a plane acoustic wave on a triaxial ellipsoid with semi-axes (0.75, 1, 0.5) centered at the origin of coordinates and having three different sets of parameters of the host medium and inclusion. The complex amplitude of the initial wave field of pressures has the form 0( ) = exp( 3), 0 = 0+, 1 =
+ 3, = 5.5, = 1; 7, = 30.5, = 9.5.
0, parameters of the media: ) = 8, = ) = 15.5, = 5, = 9, = 4; ) = 21, =
In all the examples, the diffraction problem was solved twice: using the
mosaic-skeleton method in GMRES and without it. In the process discretization
of equations (
ranged from 1032 to 64139. Hereinafter the accuracy of GMRES amounted to 10−7, and the accuracy of the low-rank approximations was 10−5. The level of the cluster tree varied from 4 to 5. Further, squares denote the first set of media parameters, circles stand for the second set and triangles mark the third set in all the graphs.
In this case, the solutions found approximately were compared with the
solutions found on a denser grid, i.e. with 64139 sampling points. This is because the
analytical solution of Example 1 is unknown. The relative errors of solution of
the problem in Example 1 using equation (
The time spent on the solution of SLAE depending on the order of matrix
with the mosaic-skeleton method and without it is shown in Figures 3 for
equation (
= 64139 for the third set of parameters for equation
(
The relative errors of solutions ( ) for the problem in Example 2 formulated
in the form of equation (
The values of mosaic rank and compression factor for the problem in Example
2 formulated in the form of equation (
Numerical experiments have shown that the modification of numerical algorithms for solving diffraction problems with the mosaic-skeleton method leads to a significant acceleration of their work while maintaining the same accuracy of calculations.