=Paper=
{{Paper
|id=Vol-1839/MIT2016-p39
|storemode=property
|title= Using of the mosaic-skeleton method for numerical solution of three-dimensional scalar diffraction problems
|pdfUrl=https://ceur-ws.org/Vol-1839/MIT2016-p39.pdf
|volume=Vol-1839
|authors=Sergei Smagin,Aleksei Kashirin,Mariia Taltykina
}}
== Using of the mosaic-skeleton method for numerical solution of three-dimensional scalar diffraction problems==
https://ceur-ws.org/Vol-1839/MIT2016-p39.pdf
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
Using of the Mosaic-skeleton Method for
Numerical Solution of Three-Dimensional Scalar
Diffraction Problems
Sergei Smagin, Aleksei Kashirin, and Mariia Taltykina
Computing Center of FEB RAS, Khabarovsk, Russia
smagin@as.khb.ru,elomer@mail.ru,taltykina@yandex.ru
http://www.ccfebras.ru/en/general
Abstract. In the paper three-dimensional stationary scalar diffraction
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 equa-
tions 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.
Keywords: diffraction problem, Helmholtz equation, integral equation,
mosaic-skeleton method, GMRES
1 Introduction
Three-dimensional stationary diffraction problems of acoustic waves are widely
used in science and technology. Most often such problems are solved numeri-
cally, since their analytical solutions can be constructed only in simplest cases.
The method of numerical solution of the initial problems should take into con-
sideration that the solutions are found in unbounded domains, must satisfy the
radiation condition at infinity, and are rapidly oscillating functions at large wave
numbers. Keeping these limitations in mind, it is preferable to use methods of
the potential theory for creating the algorithms for numerical solution. In this
case, the three-dimensional problem in an unbounded domain can be reduced to
a two-dimensional problem formulated on a closed surface of inclusion.
Using the methods of potential theory, two equivalent weakly singular bound-
ary Fredholm integral equations of the first kind with one unknown function
(density) [1β3] are constructed for the problem of diffraction. A special method of
averaging kernels of weakly singular integral operators [4, 5] is used to build their
discrete analogs. In this case, the unknown density is sought as a linear combi-
nation of continuously differentiable finite functions forming a partition of unity
on the surface of inclusion. In the process of discretization of the integral equa-
tion, the surface integrals are approximated by expressions containing integrals
over space R3 , which are then calculated analytically. This allows calculation of
458
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
coefficients of systems of linear algebraic equations (SLAE), approximating cor-
responding 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 time-
consuming 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 mosaic-
skeleton 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 numeri-
cal 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 calcu-
lation 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 nu-
merical solution of such problems. It also presents the results of numerical ex-
periments, which allow us to judge the effectiveness of the approach.
2 The initial problem and its equivalent integral
equations
Let us consider a three-dimensional Euclidean space R3 with orthogonal coordi-
nate system ππ₯1 π₯2 π₯3 , filled with a homogeneous isotropic medium with density
ππ , with speed of propagation of acoustic oscillations ππ and absorption coeffi-
cient πΎπ , 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,
by πΊπ and πΊπ (πΊπ = R3 βπΊ Β―π ).
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
459
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
waves in domain πΊπ , and transmitted waves in domain πΊπ . Therefore, the com-
plex amplitude of the complete field of pressures π’ can be presented as:
{οΈ
π’π , π₯ β πΊπ ,
π’=
π’π + π’0 , π₯ β πΊπ ,
where π’π , π’π are complex amplitudes of the pressure field of transmitted and
reflected wave fields.
The inclusion, initial, reflected and transmitted pressure wave fields are
shown in the Figure 1.
Fig. 1. Wave propagation in the inclusion and in the containing medium
Letβs formulate the initial problem.
Problem 1. In bounded domain πΊπ of three-dimensional Euclidean space R3 and
unbounded domain πΊπ = R3 βπΊ Β―π , separated by closed surface π€ β πΆ π+π½ , π + π½ >
1, find complex-valued functions π’π(π) β π» 1 (πΊπ(π) , π₯), which satisfy integral
identities
β«οΈ β«οΈ
* 2 *
(οΈ )οΈ
βπ’π(π) βπ£π(π) ππ₯ β ππ(π) π’π(π) π£π(π) ππ₯ = 0 βπ£π(π) β π»01 πΊπ(π) , (1)
πΊπ(π) πΊπ(π)
conjugation conditions on the material interface between πΊπ and πΊπ
β¨οΈ β β©οΈ
π’π β π’+
π , π π€ = β¨π0 , πβ©π€ βπ β π» β1/2 (π€ ), (2)
β¨οΈ β©οΈ
π, ππ π β π’π β ππ π + π’π π€ = β¨π, ππ π1 β©π€ βπ β π» 1/2 (π€ ),
as well as the radiation condition at infinity
(οΈ )οΈ
β1
ππ’π /π |π₯| β πππ π’π = π |π₯| , |π₯| β β, (3)
if functions π0 β π» 1/2 (π€ ) and π1 β π» β1/2 (π€ ) are set on the boundary π€ .
460
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
Here π£ * is a complex conjugate function to π£, β¨Β·, Β·β©π€ is duality ratio on
π» (π€ ) Γ π» β1/2 (π€ ), which generalizes the inner product in π» 0 (π€ ), π’Β± β‘ πΎ Β± π’,
1/2
πΎ β : π» 1 (πΊπ ) β π» 1/2 (π€ ), πΎ + : π» 1 (πΊπ ) β π» 1/2 (π€ ) are trace operators, π β :
π» 1 (πΊπ , π₯) β π» β1/2 (π€ ), π + : π» 1 (πΊπ , π₯) β π» β1/2 (π€ ) are operators of normal
derivatives [12], π0 = π’+ +
0 , π1 = π π’ 0 ,
2
(οΈ )οΈβ§ΈοΈ β2 β1
ππ(π) = π π + ππΎπ(π) π2π(π) , Im(ππ(π) ) β₯ 0, ππ(π) = π2π(π) ππ(π) ππ(π) ,
π is a wave circular frequency, ππ(π) > 0, ππ(π) > 0, πΎπ(π) β₯ 0. The definitions of
functional spaces used hereinafter are available in [12].
1
Remark 1. If Im(ππ ) = 0, then π’π β π»loc (πΊπ , π₯).
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
(οΈ )οΈ β¨οΈ β©οΈ (οΈ )οΈ β¨οΈ β©οΈ
π΄π(π) π (π₯) β‘ πΊπ(π) (π₯, Β·), π π€ , π΅π(π) π (π₯) β‘ ππ₯ πΊπ(π) (π₯, Β·), π π€ , (4)
(οΈ )οΈ β¨οΈ β©οΈ (οΈ )οΈβ§ΈοΈ
*
π΅π(π) π (π₯) β‘ π(Β·) πΊπ(π) (π₯, Β·), π π€ , πΊπ(π) (π₯, π¦) = exp πππ(π) |π₯ β π¦| (4π |π₯ β π¦|).
The solution to problem 1 will be sought in the form of potentials
π’π (π₯) = (π΄π π) (π₯), π₯ β πΊπ , (5)
(οΈ (οΈ + )οΈ (οΈ )οΈ)οΈ
π’π (π₯) = πππ π΄π π π’π + π1 β π΅π* π’+ π + π0 (π₯), π₯ β πΊπ ,
where π β π» β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 (1) and radiation condition at infinity (3). In addition, the
fulfillment of the first of the matching conditions for them (2) automatically en-
tails the fulfillment of the second matching condition. Substituting potentials (5)
in the first matching condition, we obtain a weakly singular Fredholm integral
equation of the first kind for defining unknown density π:
β¨πΆπ, πβ©π€ = β¨π2 , πβ©π€ βπ β π» β1/2 (π€ ), (6)
where
πΆ = (0.5 + π΅π* ) π΄π + πππ π΄π (0.5 β π΅π ) , π2 = β (0.5 + π΅π* ) π0 + πππ π΄π π1 .
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
π’π (π₯) = (π΄π π) (π₯), π₯ β πΊπ , (7)
461
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
(οΈ (οΈ )οΈ (οΈ )οΈ)οΈ
π’π (π₯) = π΄π π1 β πππ π β π’π β π΅π* π0 β π’β
π (π₯), π₯ β πΊπ ,
where π β π» β1/2 (π€ ) is an unknown density, π0 β π» 1/2 (π€ ), π1 β π» β1/2 (π€ ),
πππ = ππ /ππ .
In this case problem 1 is reduced to equation
β¨π·π, πβ©π€ = β¨π0 , πβ©π€ βπ β π» β1/2 (π€ ), (8)
π· = (0.5 β π΅π* ) π΄π + πππ π΄π (0.5 + π΅π ) .
The next theorem is proved [1].
Theorem 2. Let π0 β π» 1/2 (π€ ), π1 β π» β1/2 (π€ ), πΎπ > 0 or π is not the eigen
frequency of the problem
π₯π’ + ππ2 π’ = 0, π₯ β πΊπ , π’β = 0.
Then equations (6) and (8) are correctly solvable in the class of densities π β
π» β1/2 (π€ ) and formulae (5) and (7) provide the solution to problem 1.
Remark 2. In cases where we are more interested in the wave field in domain πΊπ ,
it is preferable to use equation (6), which allows the calculation of the reflected
field by a simpler formulae. For a similar reason, if we are interested in the
transmitted wave field in domain πΊπ , it is preferable to use equation (8).
3 Numerical method
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 {οΈ (οΈ β§ΈοΈ 2 )οΈ3
βοΈ 2
1 β ππ βπ , ππ < βπ ,
ππ (π₯) = πβ²π (π₯) πβ²π (π₯) , πβ²π (π₯) =
0, ππ β₯ βπ ,
π=1
where π₯ β π€ , ππ = |π₯ β π₯β²π |, ππ β πΆ 1 (π€ ) when π€ β πΆ π+π½ , π + π½ > 1.
Approximate solutions of equations (6) and (8) will be sought on grid {π₯π },
β«οΈ β«οΈ
1
π₯π = π₯ππ ππ€ , πΒ―π = ππ ππ€ ,
πΒ―π
π€ π€
the nodes are the centers of gravity of functions ππ . We assume that for all
π = 1, 2, . . . , π inequalities are satisfied
0 < ββ² β€ |π₯π β π₯π | , π ΜΈ= π, π = 1, 2, . . . , π,
462
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
ββ² β€ βπ β€ β, β/ββ² β€ π0 < β,
where β, ββ² are positive numbers depending on π , π0 does not depend on π .
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,
(οΈ β§ΈοΈ )οΈ
ππ (π₯) = (πππ2 )β3/2 exp β(π₯ β π₯π )2 ππ2 , ππ2 = 0.5πΒ―π .
In paper [3] it is shown that for any functions π β π» 1 (R3 ) and π β π» 1 (π€ )
equality
π
βοΈ
β¨π, πβ©π€ = ππ πΒ―π (ππ , π)R3 + π(β2 ).
π=1
Approximation of the density of a single layer potential by a volume density
allows us to obtain simple formulae for approximating integral operator π΄π(π)
from (4). The theoretical justification of the presented approach is given in [4].
The integral operators from (4) on π€ are approximated by expressions [1, 5].
π
βοΈ
β¨οΈ β©οΈ
π΄π(π) π, ππ π€ β π΄ππ
π(π) ππ , π = 1, 2, ..., π, (9)
π=1
(οΈ )οΈ
π΄ππ
π(π) β‘ π΄ππ ππ(π) ,
πΒ―π πΒ―π (οΈ )οΈ (οΈ (οΈ β )οΈ (οΈ + )οΈ)οΈ
π΄ππ (π) = exp π2ππ β πΎππ
2
π€ π§ππ β π€ π§ππ , π ΜΈ= π,
8ππππ
(οΈ β (οΈ )οΈ)οΈ
πΒ―2π (οΈ 2 )οΈ 2π πΒ―π π 2 ππ
3
π΄ππ (π) = exp πππ πππ€ (πππ ) + + 2ππ β ,
4π πΒ―π πππ 3
2 2 Β±
πππ = ππ + ππ2 , πππ = 0.5ππππ , π§ππ = πππ Β± ππΎππ ,
πΎππ = πππ /πππ , π2 = β1,
β«οΈβ
2π (οΈ )οΈ (οΈ )οΈ
π€ (π§) = β β exp βπ§ 2 exp π‘2 ππ‘,
π
π§
π
βοΈ
β¨οΈ β©οΈ ππ
ππ + π΅π(π) π, ππ π€ β π΅π(π) ππ , π = 1, 2, ..., π, π = Β±0.5, (10)
π=1
463
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
β¨ β© π
βοΈ
* ππ
ππ + π΅π(π) π, ππ β π΅π(π) ππ , π = 1, 2, ..., π, (11)
π€
π=1
ππ πππ (οΈ )οΈ (οΈ )οΈ
π΅π(π) = 2
exp πππ(π) πππ ππ π(π) πππ β 1 πΒ―π πΒ―π , π ΜΈ= π,
4ππππ
3
βοΈ βοΈ πππ πΒ―ππ
ππ π₯ππ β π₯ππ
π΅π(π) = (β |π| + π + Gsπ ) πΒ―π , πππ = πππ , Gsπ = β 2
,
πππ 4ππππ
π=1 πΜΈ=π
πππ are components of the unit vector of the external normal to the surface at
point π₯π .
The operators on the left sides of equations (6) and (8) are approximated by
the composition of operators (9)β(11):
π
βοΈ π
βοΈ
ππ ππ
β¨πΆπ, ππ β©π€ β πΆππ ππ , β¨π·π, ππ β©π€ β β πΆππ ππ , π = 1, 2, ..., π, (12)
π=1 π=1
ππ
πΆππ = π΅πππ π΄ππ
π β πππ π΄ππ ππ
π π΅π ,
and the right sides of equations (6) and (8) by formulae
π
βοΈ
β¨π2 , ππ β©π€ β (πππ π΄ππ ππ
π π1π β π΅π π0π ), β¨π0 , ππ β©π€ = πΒ―π π0π , (13)
π=1
πππ = β¨ππ , ππ /πΒ―π β©π€ , π = 0, 1, π = 1, 2, ..., π.
Solving the corresponding SLAE, we find the approximate values of the den-
sity of integral equations at discretization points. After that, an approximate
solution of the diffraction problem using integral representations can be calcu-
lated at any point in space.
4 Mosaic-skeleton method
Let us introduce the definitions necessary for describing the mosaic-skeleton
method [7]. Let π΄π be block matrix π΄πΓπ , and π±(π΄π ) be matrix of size π Γ π,
resulting from π΄π by adding zeros up to π΄.
Definition 1. The systems of blocks {π΄π } will be called covering π΄, if
βοΈ
π΄= π±(π΄π ),
π
βοΈ
and mosaic partitioning π΄, if, in addition, π΄π = β
.
π
Definition 2. The mosaic rank of matrix π΄ β CπΓπ , corresponding to some
covering, is number βοΈ
mr π΄ = mem π΄π /(π + π), (14)
π
where the sum is taken over all blocks of covering π΄π β Cππ Γππ , and mem π΄π =
min{ππ ππ , (ππ + ππ ) rank π΄π }.
464
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
The mosaic rank determines the memory requirements for storing the mosaic
partitioning of matrix π΄, and also the computational complexity of matrix-vector
multiplication.
An important indicator of the effectiveness of this method is the compression
factor πΌ [7]. It is defined by the following formulae
mem π΄π
πΌ= 100, (15)
mem π΄
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 (12).
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 approx-
imation 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 ππ .
465
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
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
will be denoted by ππ+1 .
4. Along the cross with centre (ππ , ππ ) we build a skeleton so that the coefficients
of matrix
βοΈπ
(π) π’πΌ π£πΌ*
π΄ = (16)
π΄
πΌ=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 (16) is considered to be sufficiently accurate if inequality
(stopping criterion)
||π΄ β π΄(π) ||πΉ β€ π||π΄(π) ||πΉ ,
where π is a relative approximation error, || Β· ||πΉ is a Frobenius norm [15], is
satisfied.
To verify this criterion, it is necessary to calculate ππ of the elements of ma-
trix π΄, which is very time-consuming. The number of operations can be reduced
by using the upper estimate of ||π΄ β π΄(π) ||πΉ , obtained in [8]
||π΄ β π΄(π) ||πΉ β€ (π β π)||π’π π£π* ||πΉ /|π΄ππ ππ | = (π β π)||π’π ||πΉ ||π£π ||πΉ /|π΄ππ ππ |,
where π = min(π, π).
Here the stopping criterion takes the form of
(π β π)||π’π ||πΉ ||π£π ||πΉ /|π΄ππ ππ | β€ π||π΄(π) ||πΉ . (17)
Using recurrence formulae we can obtain next result
(οΈ βοΈ
πβ1
(π’* π’πΌ )(π£ * π£π ) )οΈ ||π’π ||2πΉ ||π£π ||2πΉ
π πΌ
||π΄(π) ||2πΉ = ||π΄(πβ1) ||2πΉ + 2Re + ,
πΌ=1
π΄ππΌ ππΌ π΄*ππ ππ |π΄ππ ππ |2
||π’1 ||πΉ ||π£1 ||πΉ
||π΄(1) ||πΉ = , π β₯ 2.
|π΄π1 π1 |
5 Numerical results
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 mul-
tiprocessor computing systems. Intel Fortran Compiler performs as a compiler,
466
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
Intel Math Kernel Library (Intel MKL) and implementation of the OpenMP
standard for the compiler are also used. The above mentioned software is used
in the computing cluster of CC FEB RAS (http://hpc.febras.net/). The cluster
consists of one subcontrol and thirteen compute nodes. The compute node with
four Six Core AMD OpteronTM 8431 processors, with clock frequency of 2.40
GHz and 96 GB of RAM has been involved in testing. To read more information
about the program, see [16].
Example 1. The diffraction problem of a plane acoustic wave on a triaxial el-
lipsoid 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 = π π’0 , parameters of the media: πΌ) ππ = 8, ππ =
3, ππ = 5.5, ππ = 1; πΌπΌ) ππ = 15.5, ππ = 5, ππ = 9, ππ = 4; πΌπΌπΌ) ππ = 21, ππ =
7, ππ = 30.5, ππ = 9.5.
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 (6) and (8) was carried out by formulas (12) and (13), the order
of matrix π 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 solu-
tions 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 (6) and (8) are shown in Figures 2. The
solid lines indicate errors π’π , and the dashed lines represent errors π’π . According
to the numerical experiments, the method has the second order of accuracy, as
the errors fall twofold when the order of matrix is doubled (π(π β1 ) = π(β2 )).
The errors were calculated in the norm of grid functions π»β0 (πΊπ(π) ).
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 (6). The results for equation (8) are similar. Hereinafter, the time π‘ is
presented in seconds. The dashed lines show the time for solving SLAE using the
mosaic-skeleton method, and the solid lines indicate the time for solving SLAE
without it. When using the mosaic-skeleton method in GMRES, the solution
time is 115 times less (at π = 64139 for the third set of parameters for equation
(6)) than without the fast method. When the number of nodes is doubled, the
time for solving SLAE using the mosaic-skeleton method increases 2.5 times on
average.
Tables 1 and 2 show the values of the mosaic rank and compression factor
for SLAE, approximating the integral equation (6), depending on the number
of grid points. For equation (8) almost the same results were obtained. This is
because the kernels of equations (6) and (8) have similar properties as they are
467
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
Fig. 2. Errors of solutions π’π (solid lines) and π’π (dashed lines) for the problem in
Example 1, obtained with the help of equations (6) and (8), respectively
compositions of similar integral operators. The mosaic rank was calculated using
formula (14), and the compression factor by formula (15). The numerical experi-
ments show that the mosaic rank increases like π(ln3 (π )), and the compression
factor, from a certain π , decreases like π(ln3 (π )/π ).
Table 1. Mosaic ranks for the problem in Example 1.
kβM 1032 2096 4010 8022 16033 32120 64139
I 493.1 751.9 991.0 1244.1 1517.1 1845.1 2241.0
II 509.9 818.5 1103.0 1376.8 1666.3 2004.4 2410.3
III 518.0 948.5 1327.9 1677.7 2011.6 2380.7 2810.7
Table 2. Compression factors for the problem in Example 1.
kβM 1032 2096 4010 8022 16033 32120 64139
I 95.6 71.7 49.4 31.0 18.9 11.5 7.0
II 98.8 78.1 55.0 34.3 20.8 12.5 7.5
III 100.0 90.5 66.2 41.8 25.1 14.8 8.8
Example 2. The problem of diffraction on a triaxial ellipsoid with semi-axes
(0.75, 1, 0.5) centered at the origin is considered. The incident field is generated
by a point source
π’0 (π₯) = exp (ππ |π₯ β π₯0 |)/|π₯ β π₯0 |,
468
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
Fig. 3. SLAE solution time for the problem with the mosaic-skeleton method (dashed
lines) and without it (solid lines) for equation (6) in Example 1
where π₯0 = (1.125, 1, 0.75). The sets of parameters of the containing medium
and inclusion are the same as in Example 1.
The relative errors of solutions π’π(π) for the problem in Example 2 formulated
in the form of equation (6) and (8) are shown in Figures 4. The solid lines indicate
errors π’π , and the dashed lines represent errors π’π . All the errors are calculated
in the norm of grid functions π»β0 (πΊπ(π) ). When we increase the order of the
matrices twofold, they also halve, as in the above example.
The values of mosaic rank and compression factor for the problem in Example
2 formulated in the form of equation (6) are shown in Tables 3 and 4, respectively.
The mosaic rank, just as in the previous example, grows at the rate of π(ln3 (π )),
and the compression factor falls like π(ln3 (π )/π ). Almost the same results were
obtained for equation (8).
Table 3. Mosaic ranks for the problem in Example 2
kβM 1032 2096 4010 8022 16033 32120 64139
I 493.1 751.9 991.0 1244.1 1517.1 1845.1 2241.0
II 509.9 818.5 1103.0 1376.8 1666.3 2004.4 2410.3
III 518.0 948.5 1327.9 1677.7 2011.6 2380.7 2810.7
Figures 5 show the real and imaginary parts of solution π’π for the first set
of media parameters in Example 2 on the interval β5 6 π₯1 6 5, π₯2 = 0, π₯3 = 1
and on the interval π₯1 = 0, β5 6 π₯2 6 5, π₯3 = 1, respectively. It is evident
that solution π’π oscillates rapidly and slowly decreases while distancing from the
inclusion.
469
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
Fig. 4. Errors of solutions π’π (solid lines) and π’π (dashed lines) for the problem in
Example 2, obtained with the help of equations (6) and (8), respectively
Fig. 5. Real and imaginary parts of solution π’π in Example 2
470
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
Table 4. Compression factors for the problem in Example 2
kβM 1032 2096 4010 8022 16033 32120 64139
I 95.6 71.7 49.4 31.0 18.9 11.5 7.0
II 98.8 78.1 55.0 34.3 20.8 12.5 7.5
III 100.0 90.5 66.2 41.8 25.1 14.8 8.8
6 Conclusion
We have studied the possibilities of using the mosaic-skeleton method for nu-
merical solution of three-dimensional diffraction problems of acoustic waves, for
which two new formulations have been offered in the form of boundary Fredholm
integral equations of the first kind with one unknown function. The method has
been implemented on the stage of solving linear systems. All of the most labor-
intensive computing processes have been parallelized using OpenMP. The obvi-
ous advantage of this approach compared to the others, for example, to a fast
multipole method, is that the mosaic-skeleton method can be used in already
existing software for the numerical solution of integral equations of diffraction
theory. Thereat only the algorithms for calculating the coefficients of linear sys-
tems and matrix-vector multiplication are changed. Using the multipole method
and other similar methods involves creating virtually new computing algorithms
and programs.
Numerical experiments have shown that the modification of numerical algo-
rithms for solving diffraction problems with the mosaic-skeleton method leads
to a significant acceleration of their work while maintaining the same accuracy
of calculations.
References
1. Kashirin, A.A.: Research and Numerical Solution of Integral Equations of Three-
dimensional Stationary Problems of Diffraction of Acoustic Waves. PhD thesis,
Khabarovsk (2006) (in Russian)
2. Kashirin, A.A., Smagin, S.I.: Generalized Solutions of the Integral Equations of a
Scalar Diffraction Problem. Diff. Equat. 42(1), 79β90 (2006)
3. Kashirin, A.A., Smagin, S.I.: Numerical Solution of Integral Equations of a Scalar
Diffraction Problem. Dokl. Akad. Nauk, 458(2), 141β144 (2014)
4. Smagin, S.I.: Numerical Solution of an Integral Equation of the First Kind with a
Weak Singularity for the Density of a Simple Layer Potential. Comp. Math. Math.
Phys., 28(11), 1663β1673 (1988)
5. Kashirin, A.A., Smagin, S.I.: Potential-based Numerical Solution of Dirichlet Prob-
lems for the Helmholtz Equation. Comp. Math. Math. Phys. 52(8), 1492β1505 (2012)
6. Saad, Y., Schultz, M.: GMRES: a Generalized Minimal Residual Algorithm for
Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical
Computing. 7(3), 856β869 (1986)
7. Tyrtyshnikov, E.E.: Methods for Fast Multiplication and Solution of Equations.
Matrix Methods and Computations, INM RAS, Moscow. 4β41 (1999) (in Russian)
471
Mathematical and Information Technologies, MIT-2016 β Mathematical modeling
8. Savostyanov, D.V.: Polilinear Approximation of Matrices and Integral Equations.
PhD thesis, Moscow (2006) (in Russian)
9. Goreinov, S.A.: Mosaic-skeleton Approximations of Matrices Generated by Asymp-
totically Smooth and Oscillatory Kernels. Matrix Methods and Computations, INM
RAS, Moscow. 41β76 (1999) (in Russian)
10. Tyrtyshnikov, E.E.: Mosaic-skeleton Approximations. Calcolo, 33(1), 47β57 (1996)
11. Tyrtyshnikov, E.E.: Incomplete Cross Approximations in the Mosaic-skeleton
Method. Computing, 64(4), 367β380 (2000)
12. McLean, W.: Strongly Elliptic Systems and Boundary Integral Equations. Cam-
bridge University Press, Cambridge (2000)
13. Bebendorf, M.: Approximation of Boundary Element Matrices. Numer. Math.
86(4), 565β589 (2000)
14. Kashirin, A.A., Smagin, S.I., Taltykina, M.Y.: Mosaic-Skeleton Method as Ap-
plied to the Numerical Solution of Three-Dimensional Dirichlet Problems for the
Helmholtz Equation in Integral Form. Computational Mathematics and Mathemat-
ical Physics, 4(56), 612β625 (2016)
15. Meyer, C.D.: Matrix analysis and applied linear algebra. SIAM, Philadelphia
(2000)
16. Kashirin, A.A., Smagin, S.I., Taltykina, M.Y.: Realization of Mosaic-skeleton
Method in Dirichlet problems. Informatics and control systems, 4(46), 32β42 (2015)
472