PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

Relation between convergence tolerance for iterative solvers andcalculation results in rooms using BEM
Yosuke YASUDA1; Kota SAITO2; Hidehisa SEKINE3
1, 2, 3Kanagawa University, Japan

ABSTRACTThe boundary element method (BEM) has been widely used in the field of room acoustics. In BEM, theiterative methods are often used to solve the linear system, especially when dealing with large-scale problems.
In general, an iterative solver requires a convergence tolerance, i.e., a value for termination of the iterative
process, to control the calculation accuracy. Since the calculation accuracy and time are in a trade-off relation,
it is practically important to set this value appropriately. It is often set to a sufficiently small value for safety,
but if one focuses on important results in room acoustics, specifically reverberation decay curves, frequency
responses, sound pressure level distributions, etc., it is possible to set a necessary and sufficient convergence
tolerance for each type of results. As a result, reduction of unnecessary calculation can be expected. In thisstudy, we investigate the influence of the convergence tolerance for iterative methods in BEM on calculation
results in rooms through numerical experiments. In particular, we focus on the effect of the sound absorption
condition of the surrounding walls and the diffuseness of the sound field. Finally, a recommendation value ofthe convergence tolerance for each type of results in room acoustics is proposed.
Keywords: Boundary element method, Iterative method, Convergence

1. INTRODUCTION
In the field of architectural acoustics, the boundary element method (BEM) has been widely used asa wave-based numerical method along with the finite element method (FEM) and the finite difference
method (FDM). Compared to the time-domain methods such as the finite-difference time-domain(FDTD) method, the frequency-domain methods including BEM have a great advantage that detailed
frequency characteristics of sound absorbing surfaces can be easily given.
In BEM, the iterative methods such as Krylov subspace methods are often used to solve the linear
system, especially when dealing with large-scale problems. In general, an iterative method requires aconvergence tolerance, , which is a value for termination of the iterative process, to control the
calculation accuracy. Since the calculation accuracy and calculation time are in a trade-off relation, it
is practically important to set this value appropriately. Regarding BEM for acoustic problems, theeffects of the choice of the iterative method and preconditioner on the convergence have been
investigated through numerical experiments (1-4), whereas there are few studies on the convergencecriteria or the convergence tolerance to determine when to stop the iterative process (1, 5). The
convergence tolerance is often set to a sufficiently small value for safety (  = 10 -6 is often seen (2-4, 6,7) in the literature), but if one focuses on important results in room acoustics, specificallyreverberation decay curves, frequency responses, sound pressure level distributions, etc., it is possible
to set a necessary and sufficient convergence tolerance for each type of calculation results. As a result,
reduction of unnecessary calculation can be expected.
In this paper, we investigate the effects of the convergence tolerance of the iterative methods in
BEM on the accuracy of each type of calculation results of sound fields in rooms through numericalexperiments. Especially, we focus on the relation of the sound absorption and diffuseness in rooms to
the convergence of each type of the calculation results. Finally, we propose a recommendation value ofthe convergence tolerance for each type of the results.

1

yyasuda@kanagawa-u.ac.jpr201770171pj@jindai.jp
3sekinh01@kanagawa-u.ac.jp
27454

2. NUMERICAL SETUP
2.1 Numerical Method
We use the frequency-domain BEM. At all analysis frequencies, the boundary elements areconstant and rectangular, and their width is < 1/6 of the wavelength of the upper limit frequency.
2.2 Iterative Method and Preconditioner
We use GPBiCG, one of the product-type Krylov subspace iterative methods in which thecomputation time and required memory do not increase with iteration. Its good performance for the
acoustic frequency-domain BEM has been confirmed through numerical experiments (4). GMRes iswell known as a robust iterative method for acoustic problems and has been often used in the literature.
However, it has a disadvantage that the computation time and required memory increase with iteration.
For this reason, we did not mainly deal with GMRes in this paper. Our numerical experiment revealed
that there was no significant difference in the recommendation value of the convergence tolerance
between these two iterative methods while there were differences in the convergence process (We omit
the results because of the page limitation).
We use ILU(0) as a preconditioner. In this paper, we construct the ILU(0) preconditioner from asparse matrix, which is an approximation of the system matrix, having only 200 larger non-zero entries
in each row.
2.3 Stopping Criterion for Convergence
We use the following stopping criterion:|b – A xi||2/ ||b||2 ≤ ,(1)

where A is the system matrix, b is the right-hand-side vector, x i is the i-th approximation of thesolution, || ||2 is the l 2 norm, and  is the convergence tolerance. In this paper, we set = 10 -1 , 10 -2, 10 -4,
or 10-6, and regard the results with  = 10 -6 as sufficiently converging ones.
2.4

Calculation of Results

2.4.1 Frequency responses
The interval of the analysis frequencies is f = 0.5 Hz in all numerical experiments.
2.4.2 Reverberation decay curves
We process the frequency responses obtained by BEM through 1/3-octave band filters, and obtaintransient responses by the inverse fast Fourier transform. We time-delay the transient responses tosatisfy the law of causality, and obtain reverberation decay curves for the 1/3-octave bands by the
integrated impulse response method. We calculate reverberation times T20 and T30 and early decay
times EDT from the reverberation decay curves.
2.4.3 Sound pressure level distributions
From the single-frequency results, we obtain the sound pressure levels for the 1/3-octave bands by
the energy summation.

3. RELATION BETWEEN SOUND ABSORPTION AND CONVERGENCE
In this section, we investigate the effects of sound absorption in a room, concretely the mean soundabsorption coefficient and arrangement of sound absorption surfaces, on the convergence of thecalculation results.
3.1

Investigation Setup

3.1.1 Analysis model
The analysis model is a d 3 -m3 cubic room with a point source located at (0.9d, 0.8d,0.7d), as shownin Figure 1.

7455

z

point source (0.9d, 0.8d, 0.7d)
d

d

0

d

y

unit: [m]x

Figure 1 – Analysis model (Sec. 3).
3.1.2 Boundary conditions
Table 1 shows the boundary conditions of analysis cases, where n is the normal incidence soundabsorption coefficient of the absorbing surface, and n,ave is the mean sound absorption coefficient at
normal incidence. Cases a, c, and o have absorptive surfaces on all surfaces, only on the ceiling (z =d), and on the ceiling and a wall perpendicular to the ceiling, respectively. Other surfaces arereflective (n = 0.01). Real-number surface impedance corresponding to the normal incidence soundabsorption coefficient is given to each surface.
Table 1 – Boundary conditions of analysis cases (Sec. 3).

3.2

Case

Absorption arrangement

nn,ave

a01

all

0.1

0.100

a05

all

0.5

0.500

c05

ceiling (z = d)

0.5

0.092

o05

ceiling (z = d), a wall (y = 0)

0.5

0.173

Results and Discussion

3.2.1 Reverberation decay curves, EDT, T20, and T30
Figure 2 shows the reverberation decay curves obtained from the BEM results and the Eyring
reverberation theory (Eyring theoretical lines) at (0.3d, 0.4d, 0.2d). The Eyring theoretical lines were
calculated with the random incidence sound absorption coefficients.
When the convergence tolerance is large ( = 10 -1 , 10 -2 ), remarkable bending occurs especially inhigh frequency bands. The results with  = 10 -4 sufficiently converge over the range of 60 dB.
Between two cases where all surfaces are absorptive (a01 and a05), the bending occurs earlier andthe convergence is more insufficient in Case a01, where the sound absorption coefficient is smaller.
Between two cases where the mean sound absorption coefficients are nearly the same (a01 and c05),
the convergence is slower in Case c05, where the arrangement of the sound absorption surfaces isuneven. Between two cases where the arrangement of the sound absorption surfaces is similar (c05
and o05), the convergence is slower in Case c05, where the absorption is more uneven.
As seen in Figure 2, the all values of EDT, T20 , and T30 sufficiently converge when  = 10 -4 . The
EDT values almost converge when  = 10 -2 .
3.2.2 Frequency responses
Figure 3 shows the frequency responses of (a) the sound pressure levels and (b) their absolutedifferences from the results with  = 10-6 at (0.3d, 0.4d, 0.2d). In each case, the result with  = 10 -1 isdifferent from that with  = 10 -6, whereas the result with  = 10 -4 converges well with differences ofless than 1 dB in the whole frequency range. The effects of sound absorption on the convergence aresimilar to those in the reverberation decay curves.

7456

Relative level [dB] Relative level [dB] Relative level [dB] Relative level [dB] Relative level [dB] Relative level [dB]Eyring
0
–10
–20
–30
–40
–50
–60
0
–10
–20
–30
–40
–50
–60
0
–10
–20
–30
–40
–50
–60
0
–10
–20
–30
–40
–50
–60
0
–10
–20
–30
–40
–50
–60
0
–10
–20
–30
–40
–50
–60
0

ε = 10–2

ε = 10–1
a05

a01

ε = 10–4

ε = 10–6
o05

c05

63 / d Hz

63 / d Hz

63 / d Hz

63 / d Hz

125 / d Hz

125 / d Hz

125 / d Hz

125 / d Hz

250 / d Hz

250 / d Hz

250 / d Hz

250 / d Hz

500 / d Hz

500 / d Hz

500 / d Hz

500 / d Hz

1000 / d Hz

1000 / d Hz

1000 / d Hz

1000 / d Hz

2000 / d Hz

2000 / d Hz

2000 / d Hz

2000 / d Hz

0.5dTime [s]

d

0

0.5d
Time [s]

d

0

0.5d
Time [s]

d

0

0.5d
Time [s]dFigure 2 – Reverberation decay curves for 1/3-octave bands (Sec. 3).

(b) Absolute differences from results with ε = 10-6
Difference [dB]

Relative level [dB]10(a) sound pressure levels

45
ー
d

a05

Diference [dB]c05

Difference [dB]o05

Difference [dB]

Relative level [dB]10

Relative level [dB] Relative level [dB]
10
10

a01

100
ー
d

Frequency [Hz]

1000
ー
d

2828
ー
d

15
10
5
0
–5
–10
–15

a01

15
10
5
0
–5
–10
–15
15
10
5
0
–5
–10
–15
15
10
5
0
–5
–10
–15

a05

c05

o05
45
ー
d

100
ー
d

Frequency [Hz]

1000
ー
d

2828
ー
d

Figure 3 – Frequency responses (Sec. 3): (a) sound pressure levels, and (b) the absolute differences fromresults with  = 10-6.
3.2.3 Sound pressure level distributions
Figure 4 shows distributions on the receiving plane x = 0.9d in Case c05, where the convergenceof the reverberation decay curves and frequency responses was the slowest: (a) the sound pressure7457

levels for 1/3-octave bands, and (b) the absolute level differences of results with  = 10 -1 , 10 -2, 10 -4
from those with  = 10-6 . The level differences are sufficiently small, less than 1 dB, when  ≤ 10 -2 .
This indicates that a larger convergence tolerance can be chosen for the sound pressure level
distributions for bands than for frequency responses and reverberation decay curves. The same wastrue in other cases.
(a) SPLs
125/d Hzε = 10–6dz [m]5 [dB]

0

0.2

0.4

0.6

0.8

1 [dB](b) Absolute SPL differencesfrom results with ε = 10–6
ε = 10–1
ε = 10–2
ε = 10–4

0.5d
0500/d Hz
z [m]

d
0.5d
0

2000/d Hz
z [m]

d
0.5d
0
0

0.5d
y [m]

d 0

0.5d
y [m]

d 0

c05

0.5d
y [m]

d0

0.5d
y [m]

d

Figure 4 – Distributions on receiving plane x = 0.9d (Case c05, Sec. 3): (a) sound pressure levels for
1/3-octave bands, and (b) the absolute differences from results with  = 10-6.

4. RELATION BETWEEN SOUND DIFFUSENESS AND CONVERGENCE
In this section, we install sound diffusive surfaces in a cubic room and investigate the relation
between the diffuseness of the sound field in the room and the convergence of the calculation results.
4.1

Investigation Setup

4.1.1 Analysis cases
Figure 5 shows analysis cases. The volume of each room is about d 3 m3 . A point source is located at(0.9d, 0.8d, 0.7d).
4.1.2 Diffuser shapes
Cases H and V have horizontally and vertically diffusive surfaces on one wall, respectively. Cases
HH and VV have horizontally and vertically diffusive surfaces on two orthogonal walls, respectively.
The symbol ‘r’ attached to these symbols of the cases is related to the size of the diffuser period; theperiod is designed to diffuse the sounds at more than 1000/d Hz in the cases with ‘r’ and at more than
2000/d Hz in the cases without ‘r.’
d

N-a
N-c

H-a
H-c

Hr-a
Hr-c

HH-c

HrHr-c

V-c

Vr-c

VV-c

VrVr-cd0dpoint source(0.9d, 0.8d, 0.7d)zxy
unit: [m]
0.25d

0.05d
detail

0.1d
0.5d

Figure 5 – Analysis cases (Sec. 4).

7458

4.1.3 Boundary conditions
Boundary conditions ‘-a’ and ‘-c’ denote that all surfaces are absorptive with n = 0.1 and that only
the ceiling (z = d) is absorptive with n = 0.5, respectively. Other surfaces are reflective (n = 0.01).
The mean sound absorption coefficients at normal incidence are nearly the same (0.087 ≤ n,ave ≤0.100) among all cases. Real-number surface impedance corresponding to the normal incidence soundabsorption coefficient is given to each surface.
4.2 Evaluation of Diffuseness
Here we evaluate the diffuseness of the sound field in each case by using the sufficiently
converging results with  = 10 -6. We use the diffuseness indexes based on the coefficients of variation(8, 9) and a newly introduced index based on the finite difference series of the frequency response.
4.2.1 Evaluation indexes
The coefficient of variation in the time domain (CVTD) is defined as:

CVTD 

t2
1[ g 2 (t )  g 2 (t )]2 dt
tt 2  t1 1

g 2 (t ) ,

(2)

where

g 2 (t ) t2
1
g 2 (t )dt ,tt 2  t1 1

(3)

g(t) is the decay-cancelled impulse response (9), defined asg (t )  p(t ) p ( )d ,
2

(4)

t

and p(t) is the impulse response. The lower and upper limits t 1 and t 2 for the evaluation time are set
corresponding to the limits of the whole evaluation range 0 to − 60 dB because significant bending
was hardly seen in the reverberation decay curves regardless of the case.
The coefficient of variation in the frequency domain (CVFD) is defined as:

CVFD 

f2
1[ P( f )  P( f )]2 df

f 2  f1 f1P( f ) ,

(5)

where

P( f ) f2
1
P ( f )df ,
f 2  f 1 f1

(6)

P(f) = |H(f)|2 is the power spectrum, and H(f) is the frequency response. The upper and lower limitsof the frequency band, f1 and f 2, are set to the cut-off frequencies of the 1/3-octave band.
In order to evaluate the degree of peaks and dips of a frequency response, an index calculated
from a finite difference series of the frequency response of sound pressure level is introduced: LH  1 N LH ( f n  f )  H ( f n )  
,f
 f  N n 1(7)where LH (f n) is the sound pressure level at the analysis frequency f n in the frequency band, f is theinterval of the analysis frequencies, and N is the number of the analysis frequencies in the frequency
band.
4.2.2 Comparison of indexes
Figure 6 and 7 show calculation results of these indexes. The results show different aspects ofdiffuseness from each other.

7459

63/d Hz500/d Hz125/d Hz
1000/d Hz

250/d Hz
2000/d Hz

32

2

125/d Hz
1000/d Hz

250/d Hz
2000/d Hz

VrVr-c

VV-c

Vr-c

V-c

HrHr-c

Hr-c

HH-c

H-c

N-c

Hr-a

N-a

VV-c

VrVr-cVr-cV-c

HrHr-c

Hr-c

HH-c

H-c

N-c

0

Hr-a

0

H-a

1

N-a

1

H-a

CV FD

CV TD

3

63/d Hz500/d Hz

Figure 6 – Coefficients of variation: (left) in time domain, and (right) in frequency domain (Sec. 4).
63/d Hz500/d Hz125/d Hz
1000/d Hz

250/d Hz
2000/d Hz

1

H

( L
f )

0.8
0.6
0.4

VrVr-c

Vr-c

VV-c

V-c

HrHr-c

Hr-c

HH-c

H-c

N-c

Hr-a

H-a

0

N-a

0.2

Figure 7 – Mean values of finite differences in frequency domain (Sec. 4).
4.3

Relation between Convergence Tolerance and Room Diffuseness4.3.1 Evaluation of convergence in reverberation decay curve
When the convergence is insufficient, the reverberation decay curve greatly bends due to theintegration of small errors in the late part of the reverberation decay by the integrated impulse
response method. Here we propose a new index Dm,n to evaluate the convergence of a reverberation
decay curve.
First, we evaluate the curvature of the calculated reverberation decay curve by

C m , n  Tm ,n / EDTn ,(8)

where Tm,n and EDTn are Tm and EDT calculated from the reverberation decay curve obtained under
the condition of  = 10n, respectively. m represents the evaluation range (– 5 to – m – 5 dB) of thereverberation decay curve to obtain the reverberation time.
One can evaluate the curvature of the reverberation decay curve using this index, whereas one
cannot determine whether it is caused by the convergence of the calculation or the sound fieldcharacteristics itself. In order to evaluate only the convergence of the reverberation decay curve, we
define the degree of divergence (Dm,n) as follows:

Dm , n C m,nC m , 6

1 .

(9)

When the reverberation decay curve bends attributed to the sound field itself, and when thereverberation decay curve sufficiently converges, C m,n has a large value but Dm,n is nearly 0 because
C m,n ≈ C m,-6.
4.3.2 Correlation between degree of divergence and diffuseness indexes
Figure 8 shows the calculation results of the degrees of divergence Dm,-2. The correlation of D30,-2
and D40,-2 with CVTD and CVFD is partly observed, but not clear. On the other hand, the correlation of
D30,-2 and D40,-2 with the mean value of LH/f is relatively high. This difference is probably causedby that the coefficients of variation represent the deviation of the data from the mean value while the7460mean value of the difference series represents the fluctuation of the waveform, i.e., the degree ofpeaks and dips. In other words, the convergence of the calculated reverberation decay curve by BEM
has a comparatively high correlation with the degree of peaks and dips of the frequency response.
63/d Hz125/d Hz250/d Hz1

500/d Hz

1000/d Hz

2000/d Hz5
m = 40, n = – 2

0.6

3

VV-c

VrVr-cVr-cV-c

HrHr-c

HH-c

H-c

Hr-c

N-c

VV-c

VrVr-cVr-cV-c

HrHr-c

HH-c

H-c

Hr-c

0
N-c

0
Hr-a

1
H-a

0.2

Hr-a

2

H-a

0.4

N-a

D m,n

4

N-a

Dm,n

m = 30, n = – 2
0.8

Figure 8 – Degrees of divergence in reverberation decay curves (Sec. 4).
4.3.3 Reverberation curves, frequency responses and sound pressure level distributions
As in Sec. 3, the reverberation decay curves and frequency responses converged well when  =10 -4 and the sound pressure level distributions converged well when  = 10-2 regardless of the case.
We omit the results because of the page limitation.

5. CONCLUSIONSWe recommend  = 10-4 for reverberation decay curves, T20, T30, EDT, and frequency responses ofthe sound pressure level, along with  = 10 -2 for the sound pressure level distribution for bands, as
practical convergence tolerances for room sound field analyses using BEM. In the cases of highly
diffuse sound fields, larger values than these are acceptable.

REFERENCES
1. Chen S, Liu Y, Dou X. A unified boundary element method for the analysis of sound and shell-likestructure interactions. II. Efficient solution techniques. J Acoust Soc Am. 2000;108(6):2738-2745.
2. Marburg S, Schneider S. Performance of iterative solvers for acoustic problems. Part I. Solvers andeffect of diagonal preconditioning. Eng Anal Boundary Elem. 2003;27(7):727-750.
3. Schneider S, Marburg S. Performance of iterative solvers for acoustic problems. Part II. Acceleration by
ILU-type preconditioner. Eng Anal Boundary Elem. 2003;27(7):751-757.
4. Yasuda Y, Sakamoto S, Kosaka Y, Sakuma T, Okamoto N, Oshima T. Numerical analysis of large-scale
sound fields using iterative methods part I: application of Krylov subspace methods to boundary
element analysis. J Comput Acoust. 2007;15(4):449-471.
5. Gumerov NA, Duraiswami R. A broadband fast multipole accelerated boundary element method for the
3D Helmholtz equation. J Acoust Soc Am. 2009;125:191-205.
6. Yasuda Y, Ueno S, Kadota M, Sekine H. Applicability of locally reacting boundary conditions to porous
material layer backed by rigid wall: wave-based numerical study in non-diffuse sound field with
unevenly distributed sound absorbing surfaces. Appl Acoust. 2016;113:45-57.
7. Karimi M, Croaker P, Kessissoglou N. Acoustic scattering for 3D multi-directional periodic structuresusing the boundary element method. J Acoust Soc Am. 2017;141:313-323.
8. Eda K, Sakuma T. Estimation of sound field diffuseness from temporal structure of room impulse
responses. Proc Annual Meet, Environ Eng 1 (Archit Inst Jpn); Japan 2015. p. 233-234. (in Japanese)
9. Hanyu T, Hoshi K, Nakakita T. Assessment of sound diffusion in rooms for both time and frequency
domain by using a decay-cancelled impulse response. Proc Euronoise 2018; Crete, Greece 2018, p.
2011-2018.

7461