Appendix
Printone: Interactive Resonance Simulation for
Free-form Print-wind Instruments Design
Nobuyuki Umetani
September 18, 2016ABoundary Formulation of Acoustics

To make the paper self-contained, we briefly explain the boundary formulationof the Helmholtz equation. We refer readers to the book [2] for more details ofthe BEM implementation. The Helmholtz equation (1) has a kernel
G(x, y) =exp(+ikr),4πr

where r = ||x − y||,(A.1)

which is the fundamental solution to the Dirac delta function δ(x − y). Using
this kernel function and the Neumann boundary condition, the second Stoke’s
theorem leads to the equation which the sound pressure on the surface p(x)needs to satisfyZΩ(x)p(x) +4π∂G(x, y)p(y)ds(y)S ∂n(y)=G(x, xsrc ),x∈S,(A.2)where the ∂G(x, y)/∂n(y) derivative the kernel with respect to change of y ∈ S
in the normal direction of the surface is
exp(+ikr)(x − y) · n∂G(x, y)=(1 − ikr).
2∂n(y)
4πr||x − y||(A.3)

The Ω(x) is a solid angle which takes 2π on a smooth surface, and is computedfor triangle mesh using a formula presented in [3]. In our implementation,
the sound pressure is stored at the vertices of a triangle mesh and linearly interpolated over the triangle faces. We discretize equation (A.2) using a typicalcollocation method, which formulates a linear system (3) by satisfying the equation at every vertex. We use a fifth-order Gaussian quadrature to compute
this surface integration.
Once the reflection pressure at the vertices p in (3) is solved, the pressure
value at the observation point xobs inside medium Ω is computed with thesurface integration
Zp(xobs ) = −∂G(xobs , y)p(y)ds(y) + G(xobs , xsrc ),∂n(y)S1xobs ∈ Ω.

(A.4)

Our implementation is specifically categorized as the conventional boundary
integration method (CBIM), in contrast to a more sophisticated model suchas the Burton-Miller method [1]. The CBIM often suffers from errors in thefrequency where the complementary region of the media Ω̄ = {x ∈ R3 |x ∈/ Ω}
has a fictitious resonance mode. In our simulation the complementary region Ω̄
is the solid region of the musical instrument. Since our complementary regionΩ̄ is small compared to the cavity, the fictitious resonance mode is much higher
compared to the fundamental cavity resonance frequency, and thus CBIM isadequate.
The off-diagonal (i, j)-entry of the resulting coefficient matrix Aij is approximately written as:#"
rij · ni(A.5)
Aij '3 ∆j exp(+ikrij ) (1 − ikrij ),
4πrij
|
{z}
g(γ)

where rij is a vector between i- and j-vertices, rij = ||rij ||, the ni is the unit
normal vector, Ωi is the solid angle at the i-vertex, and ∆j is one third of
the area of triangles around j-vertex. Notice the nonlinearity of the coefficientmatrix with respect to wavenumber k (see Sec. 5.1). Furthermore, the nonlinear
dependent part g(γ) is a function of γ = krij and if it is small, the linear
approximation over the wavenumber is reasonable (see Sec. 6.2). Finally, theentry is invariant under the scaling geometry with s and scaling the wave numberwith 1/s i.e., rij → srij and k → k/s (see Sec. 8).

B

The Minimum Eigenvalue Bounds the Magnitude of System’s Output from Below

Here, we show that if the minimum eigenvalue of the system’s coefficient matrix
A ∈ CN ×N is very small, the magnitude of the output p = A−1 f becomes verylarge for almost arbitrary inputs f ∈ CN . We denote the eigenpair of smallest
and second-smallest magnitude eigenvalues of A as (λ0 , p0 ) and (λ1 , p1 ). The
following relationships holds according to matrix norm theory:
1|A−1 p0 |=,|λ0 ||p0 |1|A−1 p|≥max.
|λ1 ||p|p∈CN(B.1)

0hp,p i=0Note that we assume A is non-Hermitian and invertible, which is typicallytrue for exterior acoustic problems. For arbitrary input f , the vector fk =hf , p0 ip0 /|p0 |2 is a projection of f in the direction of p0 and the vector f⊥ = f −fk
is the remaining component. The output magnitude is bounded from below as|p| = |A−1 f | = A−1 (fk + f⊥ )≥ A≥−1

fk − A|fk ||f⊥ |−.
|λ0 | |λ1 |−1(B.2)
f⊥

(B.3)
(B.4)

This relationship shows that the magnitude of output |p| will become larger than
the input as long as the magnitude of λ0 is much smaller than the magnitudeof λ1 , except for the very rare case where f is perpendicular to p0 .
2CSensitivity Derivation

To derive the sensitivity of the eigenvalue and resonance wavenumber (Equation (12)), we compute one iteration of the inverse power iteration in Algorithm 1 which is w = DA−1 v. Let the matrix D and A be perturbed by ageometric change as D + ∆D and A + ∆A, where  is a small number. Thenw changes asw + ∆w = (D + ∆D)(A + ∆A)−1 v,= (D + ∆D){A(I + A−1∆A)}

v,

(C.2)

−1v,

(C.3)

= (D + ∆D)(I + A−1∆A)' (D + ∆D)(I − A−1∆A)x,−1= DA| {z v} + |∆Dx − DA
{z=w−1=∆w(C.1)−1A−1(C.4)
∆Ax .
}(C.5)

We use the Neuman expansion for the transformation from (C.3) to (C.4). Weignore the term ∆Dx since we observed that its contribution is very small.

References[1] Burton, A., and Miller, G. The application of integral equation methods to the numerical solution of some exterior boundary-value problems. InProceedings of the Royal Society of London A: Mathematical, Physical andEngineering Sciences (1971), vol. 323, The Royal Society, pp. 201–210.
[2] Kirkup, S. M. The Boundary Element Method in Acoustics: A Development in Fortran (Integral Equation Methods in Engineering). Integrated
Sound Software, 10 1998.
[3] Van Oosterom, A., and Strackee, J. The solid angle of a plane triangle.
IEEE Transactions on Biomedical Engineering 2, BME-30 (1983), 125–126.

3