Fast SVD-Based Linear Elastic Eigenvalue Problem Solver for Band Structures of 3D Phononic Crystals

In this article, a Fast Linear Elastic Eigenvalue Problem Solver (FLEEPS) is developed to calculate the band structures of three-dimensional (3D) isotropic phononic crystals (PnCs). In brief, FLEEPS solves in linear time complexity the smallest few eigenvalues and associated eigenvectors of the linear elastic eigenvalue problem originating from the finite difference discretization of the frequency-domain linear elastic wave equation. Notably, FLEEPS employs the weighted singular value decomposition based preconditioner to greatly improve the convergence rate of the conjugate gradient iteration, and uses the fast Fourier transform algorithm to accelerate this preconditioner times a vector, based on the structured decomposition of the dense unitary factor T of this preconditioner. Band structure calculations of several 3D isotropic PnCs are presented to showcase the capabilities of FLEEPS. The preliminary MATLAB implementation of FLEEPS is available at https://github.com/FAME-GPU/FLEEPS-MATLAB.

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

The MATLAB codes produced for the current study are available in the GitHub repository https://github.com/FAME-GPU/FLEEPS-MATLAB.

1. It is much better for the layperson to get acquainted with the component-wise velocity-stress formulation of (1), which is detailed in Sect. A, before jumping to (2).

2. The Bravais lattice is a fundamental concept to describe the repeating pattern of the crystal structure, which has been abstracted into a space-filling array of points generated by three lattice translation vectors, \(\varvec{a}_1\), \(\varvec{a}_2\) and \(\varvec{a}_3\). There are totally fourteen Bravais lattices, based on the symmetry of these periodically arranged points.

3. The reciprocal lattice can be regarded as the dual to the Bravais lattice, with lattice vectors \(\varvec{a}_1\), \(\varvec{a}_2\) and \(\varvec{a}_3\) of the Bravais lattice replaced by the reciprocal lattice translation vectors \(\textbf{b}_1\), \(\textbf{b}_2\) and \(\textbf{b}_3\) defined by \([\textbf{b}_1,\textbf{b}_2,\textbf{b}_3]=2\pi [\varvec{a}_1,\varvec{a}_2,\varvec{a}_3]^{-\top }\). The Brillouin zone is a special unit cell of the reciprocal lattice that is centered around a lattice point of the reciprocal lattice and is invariant under all symmetry operations of the reciprocal lattice and has the same volume as the primitive unit cell of the reciprocal lattice.

This work was partially supported by National Centre of Theoretical Sciences (NCTS) in Taiwan, TianHe-2 National Supercomputing Center in Guangzhou and the Big Data Computing Center in Southeast University. X.-L. Lyu was partially supported by Jiangsu Province Excellent Post-Doctoral Program 2023ZB142 in China. H. Tian sincerely thanks the hospitality of Nanjing Center for Applied Mathematics (NCAM) during the initialization of this work. T. Li was supported in parts by the National Natural Science Foundation of China (NSFC) 12371377. W.-W. Lin was partially supported by the Ministry of Science and Technology in Taiwan (MoST) 110-2115-M-A49-004-.

This work was partially supported by the National Natural Science Foundation of China (Grant No. 12371377), and the Ministry of Science and Technology in Taiwan (Grant No. 110-2115-M-A49-004-).

1. Xing-Long Lyu and Heng Tian have contributed equally to this work.

Authors and Affiliations

1. School of Mathematics and Shing-Tung Yau Center, Southeast University, Nanjing, 211189, People’s Republic of China

   Xing-Long Lyu & Tiexiang Li

2. Nanjing Center for Applied Mathematics, Nanjing, 211135, People’s Republic of China

   Xing-Long Lyu, Tiexiang Li & Wen-Wei Lin

3. College of Chemistry and Environmental Engineering, Sichuan University of Science and Engineering, Zigong, People’s Republic of China

   Heng Tian

4. Department of Applied Mathematics, National Yang Ming Chiao Tung University, Hsinchu, 300, Taiwan

   Wen-Wei Lin

Corresponding author

Correspondence to Heng Tian.

Conflict of interest

The authors have no relevant financial or non-financial interests to disclose.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The Component-Wise Formulation of (1)

In order to facilitate the reader to derive (4), we write down the component-wise velocity-stress formulation of (1) in detail as follows:

Therefore, in some sense, only six components of \(\tau (\textbf{r})\), i.e.,  \(\tau _{xx}(\textbf{r}),\tau _{yy}(\textbf{r}),\tau _{zz}(\textbf{r})\), \(\tau _{yz}(\textbf{r})\), \(\tau _{xz}(\textbf{r})\) and \(\tau _{xy}(\textbf{r})\), are actually needed to solve (1).

Definitions of \(J_2\) and \(J_3\) for 14 Bravais Lattices

For all 14 Bravais lattices, we can compute the QR factorization with column pivoting of \([\varvec{a}_1,\varvec{a}_2,\varvec{a}_3]\) such that

where \(\Pi \) is a suitable permutation, Q is orthonormal and \(R_3\) is upper triangular. Define \(S:=\text{ diag }\big \{\text{ sign }(R_3(1,1)),\text{ sign }(R_3(2,2)),\text{ sign }(R_3(3,3))\big \}\) and

whereby three nonnegative scalars \(\rho _i\) (\(i = 1, 2, 3\)) are defined as follows:

Furthermore, we define

where round denotes rounding to the nearest integer. Given a wave vector \(\textbf{k}\), using these notations, matrices \(J_2\) and \(J_3\), with those in (10b) and (11b) for the FCC lattice as the special case, for all 14 Bravais lattices are defined as follows:

1. (i)

   for \(\widetilde{\varvec{t}}_3(2) \ge 0\) and \(m_{4} \equiv m_2 - m_1 \ge 0\),

   $$\begin{aligned} \hspace{-10mm} J_3 = { \begin{bmatrix} &{} e^{-\imath 2\pi \textbf{k}\cdot ({\widetilde{\varvec{t}}}_{2} + (\rho _{1} - \rho _{2}){\widetilde{\varvec{t}}}_{1})}I_{m_{3}}\otimes J_{2,m_4} \\ e^{\imath 2 \pi \rho _2 \textbf{k}\cdot \widetilde{\varvec{t}}_1} I_{n_{2} - m_{3}}\otimes J_{2,m_2} &{} \end{bmatrix}}; \end{aligned}$$
2. (ii)

   for \(\widetilde{\varvec{t}}_3(2) \ge 0\) and \(m_{4} < 0\),

   $$\begin{aligned} J_3&={ \begin{bmatrix} &{} e^{-\imath 2\pi \textbf{k}\cdot ({\widetilde{\varvec{t}}}_{2} + (\rho _{1}-\rho _{2} - 1){\widetilde{\varvec{t}}}_{1})}I_{m_{3}}\otimes J_{2,n_1+m_4} \\ e^{\imath 2 \pi \rho _2 \textbf{k}\cdot \widetilde{\varvec{t}}_1} I_{n_{2} - m_{3}}\otimes J_{2,m_2} &{} \end{bmatrix}}; \end{aligned}$$
3. (iii)

   for \(\widetilde{\varvec{t}}_3(2) < 0\) and \(m_{5} \equiv m_1 + m_2 \le n_{1}\),

   $$\begin{aligned} \hspace{-10mm} J_3 = { \begin{bmatrix} &{} e^{\imath 2 \pi \rho _2 \textbf{k}\cdot \widetilde{\varvec{t}}_1} I_{m_{3}}\otimes J_{2,m_2}\\ e^{\imath 2\pi \textbf{k}\cdot (\widetilde{\varvec{t}}_2+(\rho _2+\rho _1)\widetilde{\varvec{t}}_1)}I_{n_{2} - m_{3}}\otimes J_{2,m_5} &{} \end{bmatrix}}; \end{aligned}$$
4. (iv)

   for \(\widetilde{\varvec{t}}_3(2) < 0\) and \(m_{5} > n_{1}\),

   $$\begin{aligned} J_3&={ \begin{bmatrix} &{} e^{\imath 2 \pi \rho _2 \textbf{k}\cdot \widetilde{\varvec{t}}_1} I_{m_{3}}\otimes J_{2,m_2}\\ e^{\imath 2\pi \textbf{k}\cdot (\widetilde{\varvec{t}}_2+(\rho _2+\rho _1-1)\widetilde{\varvec{t}}_1)}I_{n_{2} - m_{3}}\otimes J_{2,m_5-n_1} &{} \end{bmatrix}}, \end{aligned}$$

where \(J_{2,m} \equiv { \begin{bmatrix} &{} e^{-\imath 2\pi \textbf{k}\cdot \widetilde{\varvec{t}}_1} I_{m} \\ I_{n_1-m} &{} \end{bmatrix}}\).

Some Details on Derivations of Lemmas 3, 4 and Theorem 5

To derive eigendecompositions in Lemmas 3, 4 and Theorem 5 in Sect. 4.1, here we provide some necessary proofs and matrix properties.

Lemma 18

\({\begin{bmatrix} &{}\hspace{-0.5cm} e^{-\imath 2\pi \textbf{k}\cdot \varvec{a}_1} I_{q}\\ I_{n_1-q} &{} \end{bmatrix}}\!\!=\! (K_1^*)^q,q=1,\ldots ,n_1\), with \(K_1\) defined in (9). Hence, \(J_2\) defined in (10b) satisfies

Proof

By applying mathematical induction with respect to q. \(\square \)

Corollary 19

\((I_{n_2}\!\otimes \! K_1)K_2\!=\!K_2(I_{n_2}\!\otimes \!K_1),C_1C_2=C_2C_1\) and \(C_1C_2^*=C_2^*C_1\).

Proof

By direct verification and using (C1), we have \((I_{n_2}\!\otimes \! K_1)K_2\!=\!K_2(I_{n_2}\!\otimes \!K_1)\) and \((I_{n_2}\!\otimes \! K_1)K_2^*\!=\!K_2^*(I_{n_2}\!\otimes \!K_1)\), therefore, \(C_1C_2=C_2C_1\) and \(C_1C_2^*=C_2^*C_1\), using the definitions of \(C_1\) and \(C_2\). \(\square \)

Lemma 20

It holds that \({\begin{bmatrix} &{}\hspace{-1cm} e^{-\imath 2\pi \textbf{k}\cdot \varvec{a}_2} I_{q}\otimes J_2^*\\ I_{n_2-q}\otimes I_{n_1} &{} \end{bmatrix}}=(K_2^*)^q\), \(q=1,\ldots ,n_2\), with \(K_2\) and \(J_2\) in (10b).

Proof

By applying mathematical induction with respect to q. \(\square \)

Corollary 21

With \(K_2\) and \(J_3\) in (10b) and (11b), respectively, we have

Proof

\(J_3 (I_{n_2}\!\otimes \!K_1^{n_1/2})\!=\! J_3(I_{n_2}\!\otimes \! J_2^*)\!=\!(I_{n_2}\!\otimes \! K_1^{n_1/2})\! J_3\! =\! (I_{n_2}\!\otimes \! J_2^*)\! J_3\), which further equals \({\begin{bmatrix} &{} \hspace{-1.2cm}e^{-\imath 2\pi \textbf{k}\cdot \varvec{a}_2}\!I_{n_2/3}\otimes J_2^*\\ I_{2n_2/3}\otimes I_{n_1} &{} \end{bmatrix}} = (K_2^*)^{n_2/3}\), by Lemma 20. \(\square \)

Corollary 22

\(C_3\) commutes with \(C_1, C_1^*, C_2\) and \(C_2^*\).

Proof

By (C2), \(J_3\) commutes with \(I_{n_2}\otimes K_1, I_{n_2}\otimes K_1^*, K_2\) and \(K_2^*\), hence, by direct verification, we have \(C_3 C_\ell =C_\ell C_3\) and \(C_3C_\ell ^*=C_\ell ^*C_3\), \(\ell =1,2\). \(\square \)

Theorem 23

\(\{C_\ell , C_\ell ^*\}_{\ell =1}^3\) is a set of commutative matrices.

Proof

By Corollary 19, Corollary 22 and the normality of \(C_1,C_2\) and \(C_3\). \(\square \)

Therefore, \(C_\ell \) and \(C_\ell ^*\), \(\ell =1,2,3\), can be simultaneously diagonalized by a common unitary matrix. which is determined instantly.

Lemma 24

([30]) Given \(p,q\in \mathbb {N}\), let \(M(x)=\sum _{s=0}^{q-1}x^s M_s + x^q I_p\) with \(M_s\in \mathbb {C}^{p\times p},s=0,1,\ldots ,q-1\), then \(\det M(x) = \det (xI_{mq} - C_{BF})\), with

Moreover, if \(\textbf{x}\in \mathbb {C}^p\) and \(\beta _0\in \mathbb {C}\) satisfy \(M(\beta _0)\textbf{x}=0\), then the eigenvector of \(C_{BF}\) associated with the eigenvalue \(\beta _0\) is \([\beta _0^0,\ldots ,\beta _0^{q-1}]^\top \otimes \textbf{x}\).

Proof of Lemma 3

By direct verification, and noting that \(K_1\) is the companion matrix of the monic polynomial \((x^{n_1} - e^{\imath 2\pi \textbf{k}\cdot \varvec{a}_1})\). Moreover, for \(i,i'\in \mathbb {N}_1\) but \(i\ne i'\), we have \(e^{\imath \theta _i}\ne e^{\imath \theta _{i'}}\), then \(\textbf{x}_i^*\textbf{x}_{i'}=0\), since \(K_1\) is normal. \(\square \)

Proof of Lemma 4

Note that \(K_2\) defined in (10b) can be seen as the block companion matrix of the monic matrix polynomial \((x^{n_2}I_{n_1}-e^{\imath 2\pi \textbf{k}\cdot \varvec{a}_2}J_2)\). Hence, by (C1), Lemma 24 and Lemma 3, \(\det (x^{n_2}I_{n_1}-e^{\imath 2\pi \textbf{k}\cdot \varvec{a}_2}J_2)=0\) if and only if \(x^{n_2}= e^{\imath n_1\theta _i/2}\) for \(i\in \mathbb {N}_1\), i.e.,  \(x=e^{\imath \theta _{ji}}\) with \(\theta _{ji}\) defined in (19a), and the associated eigenvector of \(K_2\) is just \(\textbf{y}_{ji}\) for \(i\in \mathbb {N}_1, j\in \mathbb {N}_2\). Due to the orthogonality of \(\textbf{x}_i\), we have \((\textbf{y}_{j'i'}\otimes \textbf{x}_{i'})^*(\textbf{y}_{ji}\otimes \textbf{x}_i) =0\) if \(i'\ne i\). Moreover, for \(j,j'\in \mathbb {N}_2\) but \(j\ne j'\), we have \(e^{\imath \theta _{ji}}\ne e^{\imath \theta _{j'i}}\), then \((\textbf{y}_{ji}\otimes \textbf{x}_{i})^*(\textbf{y}_{j'i}\otimes \textbf{x}_{i})=0\), since \(K_2\) is normal. The eigen-decomposition of \(K_3\) can be similarly found. \(\square \)

Proof of Theorem 5

By Lemma 4, T consists of n orthonormal columns, hence is unitary. Moreover, (21) can be directly verified using (9), (10a), (11a), Lemmas 3, 4 and the properties of the Kronecker product. \(\square \)

Computating the WSVD of \(\Phi _m\)

Here, we briefly describe the numerical procedures to compute the WSVD of \(\Phi _m\) in (24b). Specifically, given a HPD matrix \(\Psi \!\in \!\mathbb {C}^{6\times 6}\), after multiplying \(\Phi _m\) with \(\Psi ^{1/2}\), we compute the SVD of \(\Phi _m\Psi ^{1/2}\), i.e., 

then we compute \(Q_m:=\Psi ^{-1/2}{\widetilde{Q}}_m\). Now, it is easy to see that \(Q_m\) satisfies \(Q_m^{*}\Psi Q_m=I_3\) and \(\Phi _m=P_m\texttt{diag}(\sigma _{1,m},\sigma _{2,m},\sigma _{3,m})Q_m^{*}\). In other words, the WSVD (25) of \(\Phi _m\) is obtained. In particular, if \(\Psi \) is set to \((\lambda _0 {\textbf{1}}_{3\times 3}+2\mu _0I_3)\oplus \mu _0 I_3\) following (38b), then its HPD square root and inverse square root are

respectively.

Proof of Corollary 11

Recall that it is claimed in Corollary 11 that \(\Psi \) in (38b) along with (36b) changes the RHS of (31) into (39), i.e., 

which is proved as follows.

Proof

Denote \(\Psi _{3\times 3}:=\lambda _0 {\textbf{1}}_{3\times 3}+2\mu _0I_3\). Following (33), eigenvalues of \(\mathcal W^{-1/2}(\mathcal {V}_v\oplus \mathcal {V}_f){\mathcal {W}}^{-1/2}\) are \(\mu _1/\mu _0, \mu _2/\mu _0\) and union of eigenvalues \(\beta _m\) of \(\Psi _{3\times 3}^{-1/2}(L_{v,mm}{\textbf{1}}_{3\times 3}+2M_{v,mm} I_3)\Psi _{3\times 3}^{-1/2}\) for \(m=1,\ldots ,n\). It is easy to see

thus \(\beta _m= \tfrac{3L_{v,mm}+2M_{v,mm}}{3\lambda _0+2\mu _0},\; \tfrac{M_{v,mm}}{\mu _0}\), where \(L_{v,mm},M_{v,mm}\) can be either \(\lambda _1,\mu _1\) or \(\lambda _2,\mu _2\). Hence, the equality in (39) holds by definition.

Without loss of generality, let \(\mu _1<\mu _2\), which, in view of (35), implies

If \(\mu _0=\mu _1\) and \(\lambda _0=\lambda _1\), we have \(\mu _2 / \mu _0 = \mu _2/\mu _1>1\), \((3\lambda _2+2\mu _2)/(3\lambda _0+2\mu _0)=(3\lambda _2+2\mu _2)/(3\lambda _1+2\mu _1)\) and \(\mu _1/\mu _0=1=(3\lambda _1+2\mu _1)/(3\lambda _0+2\mu _0)\). Then \(\mu _\ell /\mu _0 \) and \((3\lambda _\ell +2\mu _\ell )/(3\lambda _0+2\mu _0)\), \(\ell =1,2\), can be arranged in three possible orders as follows:

1. 1.

   \(1\le (3\lambda _2+2\mu _2)/(3\lambda _1+2\mu _1)\le \mu _2/\mu _1\), which, with (E5), implies

   $$\begin{aligned} \kappa ({\mathcal {W}}^{-1/2}(\mathcal {V}_v\oplus \mathcal {V}_f)\mathcal W^{-1/2}) = \tfrac{\mu _2}{\mu _1}< \tfrac{1}{2}\left( \tfrac{2\mu _2}{\mu _1}+\tfrac{3\lambda _2}{\mu _1}\right) \le \tfrac{\kappa (\mathcal {V}_v\oplus \mathcal {V}_f)}{2}; \end{aligned}$$
2. 2.

   \(1<\mu _2/\mu _1\le (3\lambda _2+2\mu _2)/(3\lambda _1+2\mu _1)\), which, with (E5), implies

   $$\begin{aligned} \kappa ({\mathcal {W}}^{-1/2}(\mathcal {V}_v\oplus \mathcal {V}_f)\mathcal W^{-1/2})= \tfrac{3\lambda _2+2\mu _2}{3\lambda _1+2\mu _1}< \tfrac{1}{2}\tfrac{3\lambda _2+2\mu _2}{\mu _1}\le \tfrac{\kappa (\mathcal {V}_v\oplus \mathcal {V}_f)}{2}; \end{aligned}$$
3. 3.

   \((3\lambda _2+2\mu _2)/(3\lambda _1+2\mu _1)\le 1<\mu _2/\mu _1\), which, with (E5), implies

   $$\begin{aligned} \kappa ({\mathcal {W}}^{-1/2}(\mathcal {V}_v\oplus \mathcal {V}_f){\mathcal {W}}^{-1/2})&= \tfrac{\mu _2}{\mu _1}\tfrac{3\lambda _1+2\mu _1}{3\lambda _2+2\mu _2}=\left( \tfrac{3\lambda _2+2\mu _2}{\mu _2}\right) ^{-1}\tfrac{3\lambda _1+2\mu _1}{\mu _1} \\&<\tfrac{1}{2}\tfrac{3\lambda _1+2\mu _1}{\mu _1}\le \tfrac{\kappa (\mathcal {V}_v\oplus \mathcal {V}_f)}{2}. \end{aligned}$$

On the other hand, if \(\mu _0=\mu _2\) and \(\lambda _0=\lambda _2\), we can similarly show that \(\kappa (\mathcal A)<\kappa (\mathcal {V}_v\oplus \mathcal {V}_f)/2\). Therefore, the inequality in (39) holds. \(\square \)

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions