Structure-preserving Arnoldi-type algorithm for solving eigenvalue problems in leaky surface wave propagation

Structure-preserving Arnoldi-type algorithm for solving

eigenvalue problems in leaky surface wave propagation

Tsung-Ming Huang

a

, Wen-Wei Lin

b,⇑

, Chin-Tien Wu

b

a

Department of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan

b

Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan

a r t i c l e

i n f o

Keywords: Leaky SAW Structure-preserving

Palindromic quadratic eigenvalue problem GTSHIRA

Mesh refinement

a b s t r a c t

We study the generalized eigenvalue problems (GEPs) that arise from modeling leaky sur-face wave propagation in an acoustic resonator with an infinite amount of periodically arranged interdigital transducers. The constitutive equations are discretized by finite ele-ment methods with mesh refineele-ments along the electrode interfaces and corners. The non-zero eigenvalues of the resulting GEP appear in reciprocal pairs ðk; 1=kÞ. We transform the GEP into a T-palindromic quadratic eigenvalue problem (TPQEP) to reveal the important reciprocal relationships of the eigenvalues. The TPQEP is then solved by a structure-pre-serving algorithm incorporating a generalized T-skew-Hamiltonian implicitly restarted Arnoldi method so that the reciprocal relationship of the eigenvalues may be automatically preserved. Compared with applying the Arnoldi method to solve the GEPs, our numerical results show that the eigenpairs produced by the proposed structure-preserving method not only preserve the reciprocal property but also possess high efficiency and accuracy.

Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction

Waveguide devices have been widely used in controlling and interconnecting guided electromagnetic waves. Advances in thin film technology and efficient transducers further encourage investigations of more sophisticated waveguide concepts in acoustic systems. Acoustic waveguide devices are widely employed in applications including telecommunication filters

[8,22]and sensor technologies[2]. One of the basic elements in most of the acoustic wave filters is a resonator that generally consists of reflectors that are externally coupled through two interdigital transducers (IDT) where the first IDT generates the surface waves and the second IDT detects the wave and filters the waves in desired frequency range out. The IDT is primarily

made by depositing periodic metallic grating electrodes on a piezoelectric film substrate, as shown inFig. 1(a). Extensive

theoretical and experimental works have been performed, especially on the Rayleigh surface acoustic wave (SAW)

[4,9,21,22]. Finite element simulations of piezoelectric devices in two dimensions (2D) and three dimensions (3D) have been

studied by Allik and Hughes[1], Buchner et al.[6], Koshiba et al.[16], and Lerch[18]and others. In filter design, it is

impor-tant to know the stop band width and the center frequency fcfor a given layout of the second sensing IDT, where fc¼

v

s=ks

and

v

sand ksare the wave velocity and the wave length of the incident wave, respectively. The center frequency and stop

band width can be determined visually by plotting the dispersion diagram in which an eigenvalue problem associated with each frequency in the search range must be solved.

Because of the slower propagation velocity of the Rayleigh SAW, filters based on Rayleigh SAW design are usually limited to an operational frequency range that is less than 1 GHz. For a frequency higher than 1 GHz, more recent attention has been

0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved.

http://dx.doi.org/10.1016/j.amc.2013.03.120

⇑Corresponding author.

E-mail addresses:min@ntnu.edu.tw(T.-M. Huang),wwlin@math.nctu.edu.tw(W.-W. Lin),ctw@math.nctu.edu.tw(C.-T. Wu).

Contents lists available atSciVerse ScienceDirect

Applied Mathematics and Computation

paid to the so-called leaky surface acoustic wave (LSAW) because of its faster propagation speed in crystal cuts such as 64°

YX-LiNbO3and 36° YX-LiTaO3and its higher electromechanical coupling and minimal propagation loss in crystal cuts such as

40°–42° YX-LiTaO3[8]. Searching for a better crystal cut among various piezoelectric substrates (PZT) to increase the LSAW

velocity has become one of the major issues in high frequency filter design. For each crystal cut, one must solve many eigen-value problems to plot the dispersion diagram. An efficient and accurate algorithm for solving eigeneigen-value problems resulting from the mathematical model of an LSAW resonator is desirable.

The eigenvalue problem[13]obtained from the finite element modeling of the SAW or LSAW resonance can be

repre-sented as a T-palindromic quadratic eigenvalue problem (TPQEP) of the form k2ðF>_M1 1 GÞ > w‘þ kðG >_M1 1 G þ F >_M1 1 F  M2Þw‘þ ðF >_M1 1 GÞw‘¼ 0; ð1Þ

or a generalized eigenvalue problem (GEP) of the form

M1 G F> 0   wi w‘   þ k 0 F G> M2   wi w‘   ¼ 0; ð2Þ where M>

1¼ M12 Cnn, M>2¼ M22 Cmm, F and G 2 Cnmwith m  n, and the superscript ‘‘>’’ denotes the complex

trans-pose. The scalar k 2 C is called the eigenvalue of(1) and (2), and the nonzero vectors w‘and w

> i w

> ‘

 >

are the associated

eigenvectors of(1) and (2), respectively. It is easily seen from(1)that the eigenvalues appear in reciprocal pairs ðk; 1=kÞ

(including 0 and 1).

The TPQEP(1)can be linearized as a GEP. The linearized GEP, as well as the GEP(2), is then solved by traditional methods

such as the QZ algorithm[13]or Arnoldi method. However, the reciprocal property of the eigenvalues can be easily

de-stroyed, and large numerical errors can be generated during the computation[15]. To preserve the reciprocal property of

the eigenvalues, a structure-preserving doubling algorithm for solving the TPQEP was developed in[11]via the computation

of a solvent of a nonlinear matrix equation. Another structure-preserving algorithm based on the ðS þ S1_{Þ-transform[19]}

and Patel’s approach[24]was developed in[14]. The drawback of solving(1)by using these structure-preserving algorithms

is that they need to explicitly compute the coefficient matrices in(1)which require large computational costs. To remedy

these drawbacks, we transform the GEP(2)into a new TPQEP of the form

PðkÞwi ðk 2 A>1þ kA0þ A1Þwi¼ 0 ð3aÞ with A> 1 ¼ FM 1 2 G >_{; A} 0¼ FM12 F > þ GM12 G >  M1 ð3bÞ

so that the structure-preserving algorithm [14]using the ðS þ S1_{Þ-transform and the generalized T-skew-Hamiltonian}

implicitly restarted Arnoldi method (GTSHIRA) can be applied to search eigenvalues in a specified region of interest. On

the basis of the shift-and-invert technique, the desired eigenpairs can be easily extracted. For giving a shift value

s

, in each

iteration of GTSHIRA, the shift-and-invert technique leads to solve a linear system ð

s

2_A>

1þ

s

A0þ A1Þx ¼ b. From (3b),

although the coefficient matrices A0and A1are large, they are not sparse. In order to avoid the computations of A0and A1

in(3b), we apply the Sherman–Morrison–Woodbury formula to directly solve the linear system. Compared with the

tradi-tional Arnoldi method for solving GEP(2), our proposed structure-preserving method not only preserves the reciprocal

prop-erty but also possesses a high efficiency and accuracy.

This paper is organized as follows. We first introduce finite element modeling for a simple resonator in Section2. In

Sec-tion3, we introduce the efficient structure-preserving algorithm to solve the large and sparse generalized eigenvalue

prob-lems resulting from our FEM model. Our numerical experiments in Section4show that the proposed structure-preserving

algorithm for solving the GEP in(2)is efficient and accurate. Finally, we conclude this paper in Section5.

2. Finite element model for SAW

In contrast to the well-known Rayleigh waves, which consist of partial longitudinal waves and shear waves, the LSAW propagates mainly in the shear direction on the sagittal plane and is trapped at the substrate surface and satisfies the

stress-free boundary condition on the surface. These properties allow one to reduce the general mode analysis in 3D to a 2D

problem, as shown inFig. 1(b)[25]. Furthermore, the boundary conditions for displacement can naturally be set to be rigid

on the bottom boundary and stress-free on the top surface, and the boundary conditions for the electric potential can be set

to be short-circuited for the electrodes on the top boundary and open-circuited elsewhere[10]. As proved in Auld’s book[4],

these boundary conditions guarantee mode orthogonality and further ensure that the mode excitation is determined by the applied traction force and the potential on the free surface. Therefore, on the sagittal plane, the usual 2D mode analysis can be applied to analyze the LSAW on the resonators with IDTs. In the following, we consider only the LSAW resonator on a 2D plane (the sagittal plane associated with the crystal cuts).

To model the wave propagation in an infinite domain with periodically arranged electrodes because of the Floquet–Bloch Theorem, we can reduce the problem to a single cell domain with one IDT by assuming that the wave w is quasi-periodic of the form wðx1;x2Þ ¼ wpðx1;x2ÞeðaþibÞx1, where x1is the wave propagation direction, p is the length of the unit cell (i.e., the

peri-odic interval),

a

and b are the attenuation and phase shifts along the wave propagation direction, respectively, and wp

sat-isfies wpðx1þ p; x2Þ ¼ wpðx1;x2Þ. LetXdenote the piezoelectric transducer with a single IDT as shown inFig. 2, and letCland

Crdenote the left and right boundary segments ofX. For the general anisotropic PZT substrates, under the assumption of

linear piezoelectric coupling, the elastic and electric fields interact following the general material constitutions below: T ¼ cE_{S  e}>_E;

D ¼ eS þ

e

S_E; ð4Þ

where vectors T; S; D and E are the mechanical stress, strain, dielectric displacement and electric field, respectively, and the

matrices cE_;

_e

S_{and e are the elasticity constant, the dielectric constant and the piezoelectric constant matrices, respectively,}

measured at constant electric and constant strain fields at a constant temperature. For various crystal cut of the PZT, the

material constant matrices cE_;

_e

S_{and e depend on the Euler angle h of the cut. By applying the Bond strain transformation}

matrix Nh[5]and the usual coordinates transformation matrix Mhto the strain field and electric field, respectively, the

mate-rial constant matrices for the cut angle h can be obtained by cE_{:¼ ½N} h > cE 0½Nh; e :¼ ½Mh > e0½Nh; and

e

S:¼ ½Mh >

e

S 0½Mh; here cE

0;e0, and

e

S0denote the material constant matrices of the crystal cut at Euler angle h ¼ ½0;0;0.

By applying the virtual work principle to Eq.(4), the equilibrium state under the external body force f, the electrical field g

and the above mentioned boundary conditions of the LSAW resonator, we have Z X ½dS>CE½S dV þ Z X ½dS>e>_½

_r

_{/ dV þ} Z X ½

r

d/>e½S dV  Z X ½

r

d/>

e

S ½

r

/ dV þ Z X dq>

_q

_{€q dV} ¼ Z X dq>_{f dV þ}Z X ½

r

d/>g dV þ Z Cl[Cr dq>_{ðT  ~nÞ dA þ}Z Cl[Cr d/>ðD  ~nÞ dA: ð5Þ

Here,

q

is the mass density, ~n is the boundary normal, q ¼ ½u;

v

;w>_{is the displacement vector, / is the electric potential}

that satisfiesr/¼ E, S ¼ ½@u

@x;@@yv;@w@z;@@zvþ@w@y;@w@xþ@u@z;@u@yþ@@xv

>_{, and dq, d/ and dS are the corresponding virtual displacement,}

potential and strain vectors, respectively. The equation can then be discretized by finite element method[1,10]. Following

the usual free mode analysis, we consider f ¼ 0, g ¼ 0 and a time harmonic quasi periodic solution vector w_xðx; tÞ ¼ wðxÞeixt_.

The spatial function wðxÞ ¼ ½qðxÞ; /ðxÞ satisfies the boundary conditions shown inFig. 2in which the periodic boundary

con-ditions, proposed by Buchner[6],

wr¼ wðx1þ p; x2Þ ¼

c

wðx1;x2Þ ¼

c

wl; ð6Þ

Tr nr¼ TðwrÞ  nr¼ 

c

TðwlÞ  nl¼ 

c

Tl nl; ð7Þ

Dr nr¼ DðwrÞ  nr¼ 

c

DðwlÞ  nl¼ 

c

Dl nl; ð8Þ

are enforced on the left and right boundaries,ClandCr, here nland nrare the normal vector ofClandCrrespectively and

c

¼ eðaþibÞ_{. By plugging w}

winto(5), the equation can be rewritten in the following matrix form:

Fig. 2. A 2D single cell domain of a LSAW resonator and boundary conditions where the region DEFG is a aluminium electrode and the region ACHB is the PZT substrate.

Kqq

x

2_Mqq Kq/ K/q K// " # q /   ¼ Flþ Fr Qlþ Qr   ; ð9Þ where Kqq¼ Z X dq>_B> SC E BSq dV; Mqq¼ Z X dq>

_q

_{q dV; and K}//_¼Z X d/>B> /eBSq dV; Kq/ ¼ Z X dq>_B> Se>B//dV and K/q¼ Z X d/>B> /eBSq dV; Fl¼ Z Cl dq>_T l~nldA and Fr¼ Z Cr dq>_T rn~rdA; Ql¼ Z Cl d/>Dln~ldA and Qr¼ Z Cr d/>Drn~rdA;

and the matrices BS¼

@ @x 0 0 0 @ @z @ @y 0 @ @y 0 @z@ 0 @x@ 0 0 @ @z @ @y @ @x 0 2 6 4 3 7 5 > and B/¼ @x@ @y@ @z@ h i>

. Mechanical damping effects can also be considered

by using the Rayleigh damping assumption in which the matrix Kqq_

_x

2_Mqq _{in (9) are modified into}

Kqqþ i

x

ð

j

1Kqqþ

j

2MqqÞ 

x

2Mqq. Here

j

1and

j

2are coefficients associated with the viscous damping and mass damping,

respectively.

To obtain the generalized eigenvalue problem that is associated with the propagation parameter

c

, following Hofer’s

ap-proach[13], the nodal unknowns are split into the inner nodes wi¼ ½qi;/i, the left boundary nodes wl¼ ½ql;/l and the right

boundary nodes w_r¼ ½qr;/r. The matrix Eq.(9)can be recast into the following form:

Kii Kil Kir Kli Kll Klr Kri Krl Krr 2 6 4 3 7 5 wi wl wr 2 6 4 3 7 5 ¼ 0 Rl Rr 2 6 4 3 7 5; ð10Þ

where Rland Rrare vectors obtained from the discretization of the terms Flþ Qland Frþ Qr, respectively. From the periodic

boundary conditions(6), (7) and (8), (10)becomes the following:

Kii Kil Kir Kli Kll Klr Kri Krl Krr 2 6 4 3 7 5 Ii 0 0 Il 0

c

Il 2 6 4 3 7 5 w_wi l   ¼ 0 Il 

c

Il 2 6 4 3 7 5Rl: ð11Þ

Furthermore, by multiplying the matrix

Ii 0 0

0

c

Il Il

 

with(11), the GEP associated with the propagation parameter

c

is obtained:

Kii Kil K> ir 0   þ

c

0 Kir K> il Kllþ Krr     _w i wl   ¼ 0: ð12Þ

Because the viscosity is small for most crystalline solids, the attenuation factor

a

is close to zero. As a result, the

propa-gation factors k near the unit circle, denoted by U, are desired. Moreover, for the frequency

x

in the stopping band, the

fre-quency shift parameter b will be close to

p

when the periodic interval p (i.e., the domain width here) equals half of the

incident wave length k0. Therefore, for the eigenvalue problem(12), we are interested in finding the eigenvalues k that

are close to U, especially, those that are near 1 on the complex plane.

Note that it is well known that the solution of a general elliptic problem has singularities near to corners[12], and in

addi-tion, the solution can become less regular near the interface between the electrode and the PZT substrate. It is inevitable that the error from discretization could be amplified when the reciprocal pairs of eigenvalues are computed. Therefore, it is important to minimize the accuracy deterioration from singularities and low regularity in finite element solutions. One can resolve the singularity by constructing the singular elements in which the mesh points are clustered to the singular

source according to the order of the singularity[3]. In our calculation, we simply employ locally refined meshes. An

addi-tional benefit from using locally refined meshes is that we can discretize Eq.(9)using linear elements instead of using high

order finite element discretization[6]. Due to the simple discretization scheme on locally refined meshes, the GEP becomes

very large sparse, an accurate and robust eigensolver for preserving the reciprocal eigenpair structure is needed. We present our algorithm and numerical results in the following sections.

3. Structure-preserving Arnoldi-type algorithm

Observing that the imaginary part of Kllþ Krr in (12)is symmetric and positive definite, by the Bendixson Theorem,

wl¼  1

c

ðKllþ KrrÞ 1 K> irwi ðKllþ KrrÞ1K>ilwi: ð13Þ Letting M1¼ Kii; M2¼ Kllþ Krr; F ¼ Kir; G ¼ Kil; x ¼ wi; k¼

c

; ð14Þ

and substituting(13)into the first equation of(12), we obtain the TPQEP in(3). To solve(3)in a structure-preserving way, we

first transform the TPQEP in(3a)into a >-skew-Hamiltonian pencil ð bK; cN Þ through the following procedure:

1. The TPQEP is linearized into a special GEP[14],

ðM  kLÞ x y   ¼ 0; ð15Þ where ky ¼ A1x, and M ¼ A1 0 A0 I   ; L ¼ 0 I A>1 0   : ð16Þ

The reciprocal eigenpairs ðk; 1=kÞ are kept in the matrix pencil ðM; LÞ because the matrix pencil M  kL is >-symplectic,

i.e., it satisfies MJM>_{¼ LJL}>_{where J ¼} 0 In

In 0

 

.

2. Using the ðS þ S1_{Þ-transform, the matrix pencil M  kL is further transformed to a >-skew-Hamiltonian pencil K }

_l

_{N ,}

i.e., ðKJ Þ>

¼ KJ , ðNJ Þ>¼ NJ and

l

is the eigenvalue of the pencil:

K 

l

N  ½ðLJM>_{þ MJL}>_{Þ }

_l

_LJL>_J>_{: ð17Þ}

From the relationship

l

¼ k þ 1=k , one can relate the two eigenvalues k and

l

and further implies that the multiplicity of

the eigenvalue

l

is even.

3. Let

s

be a shift value and

s

R

r

ðM; LÞ where

r

ðA; BÞ denotes the set of all eigenvalues of any matrix pair ðA; BÞ.

Since

l

0

s

þ 1=

s

R

r

ðK; N Þ, one can define the shift–invert transformation K b

l

bcN for K 

l

N with

b

l

¼ ð

l



l

0Þ 1 where b K  

s

N ¼

s

A > 1 0 0 A1 " # ; ð18aÞ c N  

s

ðK 

l

0N Þ ¼ ðM 

s

LÞJ ðM >_

_s

_L>_ÞJ>_{; ð18bÞ}

and bK and cN are >-skew-Hamiltonian.

The relationship between eigenpairs of the TPQEP in(3)and the >-skew-Hamiltonian pencil ð bK; cN Þ in(18)is stated in the

following theorem.

Theorem 3.1 [14]. Let ð bK; cN Þ be defined in(3)and

s

be a shift value with

s

R

r

ðM; LÞ. If ðb

l

;½z> 1;z>2

>

Þ with z1;z22 Cnis an

eigenpair of ð bK; cN Þ and

m

satisfies

s

þ1

sþ_b1_l¼

m

þ1m, then z1þ1mz2and z1þ

m

z2are eigenvectors of the TPQEP in(3)corresponding

to eigenvalues

m

and1

m, respectively.

From the definition of cN in(18b), cN can be factorized as cN ¼ N1N2, where

N1¼ M 

s

L; N2¼ J ðM>

s

L>ÞJ> ð19Þ

are nonsingular and satisfy N>2J ¼ JN1.Let B  N11 KNb 1

2 and u1be an initial vector. Define the Krylov matrix with respect

to u1by

Kn Kn½B; u1 ¼ ½u1;Bu1; . . . ;Bn1u1:

Theorem 3.2 guarantees that the Arnoldi process can be executed in a way that the >-skew-Hamiltonian structure of the matrix pencil is preserved. As a result, a generalized >-skew-Hamiltonian implicitly restarted Arnoldi (GTSHIRA) algorithm

proposed in[14]can be employed to solve the eigenvalue problem bKz ¼b

l

cN z.

Theorem 3.2 [14]. Let B ¼ N11 KNb

1

2 with cN ¼ N1N2be >-skew-Hamiltonian and Kn Kn½B; u1 be the Krylov matrix with

rankðKnÞ ¼ n. Then there are unitary matrices U and V satisfying V ¼ J>UJ , Ue1¼ u1and Ve1¼ N1u1=kN1u1k2such that

V>_{KU ¼b} Hbn bSn 0 Hb> n " # ; V>_{cN U ¼} bRn bTn 0 bR> n " # ; ð20Þ

The unitary matrices U and V in Theorem 3.2 can be generated in a structure-preserving way (GTSHIRA) as follows. Recall

that the >-bi-isotropic orthonormal matrices bZj; bYj2 C2njare computed iteratively according to the following

structure-preserving Arnoldi process: b K bZj¼ bYjHbjþ bhjþ1;jbyjþ1e>j ð21Þ and c N bZj¼ bYjbRj ð22Þ with b YH jbyjþ1¼ 0 and bZ>jJbyjþ1¼ 0; ð23Þ

where bHj; bRj2 Cjj are unreduced upper Hessenberg and nonsingular upper triangular, respectively. By defining

Uj ½bZj;JYbj and Vj ½Ybj;J bZj whereYbjdenotes the conjugate matrices of bYj, it is easily seen that

V> j KUb j¼ b Hj bYHjKJb  b Yj 0 Hb> j 2 4 3 5; V> jcN Uj¼ bRj bYHj cN J  b Yj 0 bR> j 2 4 3 5

which implies that the >-skew-Hamiltonian property is preserved in each iteration step.

Note that Theorem 3.1 indicates that although the number of the eigenvectors associated with the eigenvalue

l

is even,

only half of the eigenvectors are needed to extract all the eigenvectors corresponding to the eigenvalues

m

and1

m. Further-more, through the above mentioned structure-preserving Arnoldi process, the >-skew-Hamiltonian structure of the matrix

pencil in Theorem 3.2 is preserved and the even multiplicity of the eigenvalue

l

is automatically obtained. Therefore, the

required halves of the eigenvectors associated with the eigenvalue

l

can be easily computed when the desired eigenpairs

are convergent. In fact, the desired eigenpairs ðb

l

i;ziÞ of ð bK; cN Þ can be computed from the matrix pair ð bHj; bRjÞ with

b

Hj^si¼

l

bibRj^siand zi¼ bZj^si. From Theorem 3.1, one can compute the desired eigenpairs of ðM; LÞ from ðb

l

i;ziÞ and preserve

the reciprocal relationship of the eigenvalues of the GEP algebraically.

We summarize the above procedures for computing the reciprocal eigenpairs of the GEP(12)in Algorithm 1. Note that, in

step 1 of Algorithm 1, the linear systems

N1

v

1¼ b1; N2

v

2¼ b2; ð24Þ

have to be solved in order to obtain bZjfrom(22). This is indeed the most time-consuming step in the proposed

structure-preserving algorithm. In the following, we discuss how to solve(24)efficiently.

Algorithm 1. Structure-preserving algorithm for solving GEP(12)

Input: matrices F; G; M2and M1, shift value

s

and the number m of desired eigenvalues.

Output: eigenpairs fð

c

j;½ðw ð1Þ i;j Þ > ;ðwð1Þ_l;jÞ>>Þ, ð

c

1 j ;½ðw ð2Þ i;jÞ > ;ðwð2Þ_l;jÞ>>Þgm

j¼1of the GEP in(12)where

c

jþ

c

1j for j ¼ 1; . . . ; m are the closest to shift value

s

þ

s

1.

1: Compute eigenpairs fð

l

bj;zj ½z>j1;z>j2 >

Þgm

j¼1of ð bK; cN Þ by using GTSHIRA.

2: Compute eigenvalues

c

jandc1jof TPQEP in(3)by solving

c

2_{ ð}

_s

_þ

_s

1_þb

_l

1 j Þ

c

þ 1 ¼ 0; Compute eigenvectors wð1Þi;j  1

c

j zj1 zj2; wð2Þi;j 

c

jzj1 zj2

corresponding to

c

j;cj1, respectively, for j ¼ 1; 2; . . . ; m.

3: Compute wð1Þ_l;j ¼ M12 ð

c

1j F > wð1Þi;j þ G > wð1Þi;jÞ; w ð2Þ l;j ¼ M 1 2 ð

c

jF > wð2Þi;j þ G > wð2Þi;jÞ for j ¼ 1; . . . ; m.

By the definitions of M and L in(16), we have I 

s

I 0 I   ðM 

s

LÞ ¼

s

2_A> 1þ

s

A0þ A1 0 A0

s

A>1 I " # ð25aÞ and I A0

s

A1 0 I   ðM>

s

L>Þ ¼

s

2_A 1þ

s

A0þ A>1 0 

s

I I " # : ð25bÞ

From(19)and(25), we see that solving(24)is equivalent to solving

ð

s

2_A> 1þ

s

A0þ A1Þ

v

11¼ b11

s

b12; ð26Þ

v

12¼ b12 ðA0þ

s

A>1Þ

v

11; and ð

s

2_A 1þ

s

A0þ A>1Þ

v

22¼ b22þ ðA0þ

s

A1Þb21; ð27Þ

v

21¼

s

v

22 b21; where

v

1¼ ½

v

>11;

v

>12 > ;

v

2¼ ½

v

>21;

v

>22 > ;b1¼ ½b>11;b > 12 > and b2¼ ½b>21;b > 22 >

. By the definitions of A0and A1, it holds that

s

2_A> 1þ

s

A0þ A1¼ ðG þ

s

FÞM12 ðF > þ

s

G> Þ 

s

M1 ð28Þ and

s

2_A 1þ

s

A0þ A>1 ¼ ðF þ

s

GÞM 1 2 ðG > þ

s

F> Þ 

s

M1: ð29Þ

Let M1¼ LU be the LU factorization of M1. Set

E1¼ L1

1

s

G þ F

 

; E2¼ U>ðF þ

s

GÞ:

By the Sherman–Morrison–Woodbury formula,(28) and (29)can be further factorized as following,

ð

s

2_A> 1þ

s

A0þ A1Þ1¼  1

s

U 1 ½I  E1M12 E > 2 1_L1 ¼ 1

s

U 1 ½I þ E1ðM2 E>2E1Þ1E>2L 1_; ð30Þ and ð

s

2_A 1þ

s

A0þ A>1Þ 1 ¼ 1

s

L > ½I þ E2ðM2 E>1E2Þ1E>1U >_: ð31Þ

Now, obviously, the solutions of(26) and (27)can be obtained by two forward substitutions (L1_{), two backward}

substitu-tions (U1_{) and solving small linear systems ðM}

2 E>2E1Þ1and ðM2 E>1E2Þ1. As a result, Algorithm 1 is very efficient.

4. Numerical results

In this section, we first conduct numerical results to validate the convergence of our finite element model. Second, we report numerical comparisons with our structure-preserving method and the traditional Arnoldi method for solving the

GEP(12)to demonstrate the accuracy and efficiency of the proposed eigenvalue solver. All of the computations are

per-formed in MATLAB 2010b on an HP workstation with an Intel Quad-Core Xeon X5570 2.93 GHz and 60 GB of main memory, using IEEE double-precision floating-point arithmetic.

The configuration of our computational domain, shown inFig. 2, is as follows. The domain width AB and the height CD are

set to be 106and 3  106, respectively. The ratio of the electrode width EF versus the domain width is set to be1

2and the

ratio of the electrode thickness DE versus the domain height is1

15. The material constants of LiTaO3and LiNbO3are taken from

the measurements obtained by Kushibiki, Takanaga and Sannomiya [17]. Additionally, the viscous damping coefficient

j

1 Oð108Þ in the 10 KHZ operation range and

j

1 Oð1010Þ in the MHZ operation range for a family of PZT materials

has been shown in[20,23]. In general,

j

1depends on the operation frequency

x

. The viscous damping coefficient is

extrap-olated to the GHz operation range according to the reciprocal rule

j

1/_x1[7]. In our numerical studies, the viscous damping

coefficient

j

1is set to be 1014, and the mass damping coefficient

j

2is taken as 1 

j

1to account for the effect of the

elec-trode weight.

4.1. Accuracy and convergence of finite element approximation

First, we show that our finite element model gives accurate results when predicting the center of the stopping band of

YX-LiNbO3. The dispersion diagrams of the attenuation constant

a

and the propagation constant b associated with the

eigen-value kð

x

Þ, which is the closest to 1 on the complex plane for the frequency

x

around the stopping bands, are shown in

Fig. 3(a) and (b) for the crystals 36° YX-LiTaO3and 64° YX-LiNbO3, respectively. A typical shear wave displacement

associ-ated with the eigenvector of the computed eigenvalue is shown inFig. 4which mesh is locally refined twice near the

inter-face over a uniform mesh.

To measure the convergence of the eigenvalues, tests over three successively refined meshes with an initial mesh length

of h ¼ p

20are performed. The dimensions of matrices M1and M2associated with uniform meshes and locally refined meshes

are listed in the second and third columns ofTable 1, respectively. We set k½i;nto be the eigenvalue obtained from the meshes

with a mesh length of p=ð10  iÞ. Here, the index n ¼ u and n ¼ ‘ denote that the mesh is uniform without and with local

refinement, respectively. Using k½16;u as an exact value, the convergence of the eigenvalues can be verified from

jk½16;u k½2i_;uj and jk½16;u k½2i_;‘j for i ¼ 1; 2; 3. The values of jk½16;u k½2i_;uj and jk½16;u k½2i_;‘j ði ¼ 1; 2; 3Þ for 36° YX-LiTaO3

and 64° YX-LiNbO3at

x

sand

x

eare shown inTable 2, where

x

s and

x

eare the frequencies for which the stopping band

starts and ends, respectively. FromTable 2, it can be seen that the accuracy of the eigenvalue is increased as the mesh length

is reduced and is improved with locally refined meshes. The average of the error reduction rate is about 3.68 which is cal-culated by averaging the ratio jk½16;uk½2i ;uj

jk½16;uk_{½2iþ1 ;u}j, for i ¼ 1; 2, over both cases. Moreover, it is known that the wave propagation

veloc-ity is approximately 4112 m/s for 36° YX-LiTaO3and approximately 4478 m/s for 64° YX-LiNbO3. It is clear that because the

domain width p ¼ 1  106_{, the center of the stopping band is approximately 2.056 GHz and 2.239 GHz for 36° YX-LiTaO}

3

and 64° YX-LiNbO3, respectively. We compute the center of the stopping band by averaging

x

sand

x

eon different mesh

lengths and show the computational results inTable 3. Obviously, one can see that the central frequency monotonically

con-verges to a constant when the mesh length is reduced. The numerical error from our finite element simulations is less than

0.2% and 1.2% for 36° YX-LiTaO3and 64° YX-LiNbO3, respectively.

4.2. Comparison of Algorithm 1 and traditional Arnoldi method

FromTables 2 and 3, we have already shown that the accuracy of the computed eigenvalues and the central frequency of the stopping band obtained from the locally refined mesh with a mesh length of p=80 is almost the same as the values ob-tained from a uniform mesh with a mesh length of p=160. Therefore, in the following numerical computations, we consider

only the coefficient matrices in the GEP(12)that are generated by the finite element discretization on the mesh that is

lo-cally refined twice over the uniform mesh with a mesh length of p=80.

Let the pair ðkk;I;kk;OÞ, k ¼ 1; . . . ; N, denote the reciprocal pairs of eigenvalues of(12), where kk;Iand kk;Olie inside and

out-side of U, respectively. Fig. 5 displays the eigenvalues fk1;I; . . . ;k9;I;k1;O; . . . ;k5;Og of the LiNbO3 at a frequency of

x

¼ 2:180 GHz in which reciprocal pairs ðkk;I;kk;OÞ for k ¼ 1; . . . ; 5 close to U may be of interest. The notation OðkÞ represents

the set of all of the eigenvalues that cluster at the origin of the complex plane. Suppose that 2N eigenvalues near U are

de-sired. The Arnoldi process in(21) and (22)for GTSHIRA is set to restart if the desired eigenpairs are not convergent when the

dimension j of the subspace spanfbYjg grows more than 5N. The number for restarting the Arnoldi process is denoted by

‘‘#Iter’’ in the following.

A standard iterative approach for solving the GEP(12)is to apply the Arnoldi method on the equation directly. However,

the reciprocal property of the eigenvalues is not guaranteed to be preserved in the computation. For Algorithm 1, based on

the ðS þ S1_{Þ-transform, if k and}

_l

_{are the eigenvalues of(12) and (17), respectively, then k and}

_l

_{satisfy the relationship}

l

¼ k þ k1_{. As a result, we can obtain the k-th reciprocal pair ðk}

k;I;kk;O 1=kk;IÞ by solving the algebraic equation

−0.05 0 0.05 2.01 2.02 2.03 2.04 2.05 2.06 2.07 2.08 2.09 x 109 α ω 0.98 1 1.02 2.01 2.02 2.03 2.04 2.05 2.06 2.07 2.08 2.09x 10 9 β / π ω −0.1 0 0.1 2.16 2.18 2.2 2.22 2.24 2.26 2.28x 10 9 α ω 0.95 1 1.05 2.16 2.18 2.2 2.22 2.24 2.26 2.28x 10 9 β / π ω

l

k¼ kk;Iþ k1k;I after the k-th eigenvalue

l

kof Kz ¼

l

N z is computed. Hence, the reciprocity is automatically preserved. Two

numerical comparisons on preserving the reciprocal property, between Algorithm 1 and the traditional Arnoldi method, are

listed in the following, where the eigenvalues of 64° YX-LiNbO3at a frequency of

x

¼ 2:180 GHz are computed.

The traditional Arnoldi method does not guarantee that half of the computed eigenvalues lie inside of the unit circle and

the others lie outside. For example, when we use the Arnoldi method to compute four eigenvalues (i.e., 2N ¼ 4) of(12)

that are near to 1, the four convergent eigenvalues are k1;I;k1;O;k2;I and k3;I. Clearly, the reciprocal property of the

eigenvalues is lost.

Fig. 4. A shear wave displacement.

Table 1

Dimension information of matrices M1and M2obtained from FEM discretization.

Mesh length Uniform Local refine

M1 M2 M1 M2 p=20 3554 183 4968 183 p=40 14,548 363 17,192 363 p=80 58,856 723 63,960 723 p=160 236,752 1443 Table 2

The values of jk½16;u k½2i;nj for different mesh lengths at frequenciesxsandxe.

x(GHz) 36° YX-LiTaO3 64° YX-LiNbO3 xs¼ 2:028 xe¼ 2:075 xs¼ 2:177 xe¼ 2:257 jk½16;u k½2;uj 0.0222 0.0161 0.0955 0.0797 jk½16;u k½2;‘j 0.0178 0.0141 0.0857 0.0751 jk½16;u k½4;uj 0.0076 0.0056 0.0299 0.0458 jk½16;u k½4;‘j 0.0036 0.0042 0.0166 0.0329 jk½16;u k½8;uj 0.0016 0.0015 0.0060 0.0087 jk½16;u k½8;‘j 0.0001 0.0006 0.0007 0.0016 Table 3

Computed center frequency fc(GHz) of stopping bands of LSAW on various meshes. Here, huand h‘denote the mesh length of meshes without and with local

refinement, respectively.

Mesh length hu¼40p h‘¼40p hu¼80p h‘¼80p hu¼160p

fc LiTaO3 2.05246 2.05222 2.05206 2.05191 2.05183

Suppose that one wants to compute the five reciprocal pairs ðkk;I;kk;OÞ for k ¼ 1; . . . ; 5 and takes a shift value

s

¼ 2:89.

The eigenvalues in OðkÞ are closer to

s

than the eigenvalues k5;Iand k5;O. As a result, the desired reciprocal pair ðk5;I;k5;OÞ

near U cannot be discovered by the Arnoldi method. In contrast, in Algorithm 1, according to the relationship

l

¼ k þ k1_,

the eigenvalue

l

of Kz ¼

l

N z is far away from the shift value

s

þ 1=

s

when k is close to the origin. Naturally, Algorithm 1

will not converge to those unwanted eigenvalues in OðkÞ. As a result, all of the desired eigenvalues can be discovered

more easily by Algorithm 1 than with the traditional Arnoldi method. Our numerical results inTable 4show that not only

all of the desired eigenvalues are found by Algorithm 1, even when the number of desired eigenvalues is set to 2N ¼ 18, but Algorithm 1 also converges much faster than the traditional Arnoldi method. In fact, only two restarting steps are

needed for Algorithm 1 to converge for all of the cases shown inTable 4. In addition, from the rightmost column ofTable 4,

one can see that all of the computed eigenvalues indeed preserve the reciprocal property. In contrast, the reciprocity of the convergent eigenvalues obtained by the Arnoldi method are diminished by approximately 3 significant digits.

From the above comparison, Algorithm 1 preserves the reciprocal property of the eigenvalues of the GEP(12)effectively.

For measuring the accuracy of Algorithm 1, let us define the relative residual of an eigenpair ðk; uÞ of(12), where u ¼ ½wi;wl>,

as the following: M1 G F> ₀   u  k 0 F G> _M 2   u         F M1 G F> ₀           F kukFþ jkj 0 F G> M2           F kukF ;

where k  kF is the Frobenius matrix norm. The maximal relative residuals of the 10 desired eigenpairs for 64° YX-LiNbO3,

with various frequencies, are shown inFig. 6. From these numerical results, one can see that the eigenpairs produced by

Algorithm 1 possess high accuracy in terms of the relative residual error.

−1.5 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

Real part of eigenvalue

Imaginary part of eigenvalue

λ3,I λ3,O λ4,I λ4,O λ1,I λ1,O λ2,I λ2,O λ5,I λ5,O λ6,I λ7,I λ8,I λ_9,I Inside Outside

Fig. 5. Distribution of eigenvalues for 64° YX-LiNbO3at frequency 2.180 GHz.

Table 4

Convergent reciprocal pairs and the associated errors of reciprocity for 64° YX-LiNbO3at frequency 2.180 GHz versus different eigensolvers with various ‘‘2N’’

which denotes the number of interested eigenvalues. Here, ~kk2 OðkÞ for k ¼ 1; . . . ; 4.

Method 2N Computed eigenvalues #Iter maxfjkk;Ikk;O 1jg

Arnoldi 8 _fðkk;I;kk;OÞg4k¼1 2 1:7  1013

10 _fðkk;I;kk;OÞg4k¼1, f~kkg2k¼1 5 1:7  1013

12 _fðkk;I;kk;OÞg4k¼1, f~kkg4k¼1 4 1:7  1013

Algorithm 1 8 _fðkk;I;kk;OÞg4k¼1 2 1:1  1016

10 _fðkk;I;kk;OÞg5k¼1 2 1:1  1016

12 _fðkk;I;kk;OÞg6k¼1 2 1:1  1016

Next, let us compare the efficiency of Algorithm 1 and the traditional Arnoldi method. The CPU times for computing 10

desired eigenpairs (i.e., 2N ¼ 10) for 64° YX-LiNbO3with various frequencies with Algorithm 1 and the traditional Arnoldi

method are shown inFig. 7. On average, Algorithm 1 takes only 476 s of CPU time to compute the desired eigenpairs for

all of the frequencies

x

in the search range. Obviously, the proposed algorithm is more efficient compared with the

tradi-tional Arnoldi method, which takes 527 s of CPU time to obtain all of the desired eigenpairs. 5. Conclusions

In this paper, we have modeled leaky surface acoustic wave propagation on a simple resonator with an interdigital

trans-ducer (IDT), where electrodes are arranged periodically on piezoelectric substrates (PZT) such as 64° YX-LiNbO3and 36°

YX-LiTaO3. The energy conservation Eq.(5)is discretized by a finite element method (FEM) applied to a single cell domain with

proper periodic boundary conditions, as shown inFig. 2. Eq.(5)is discretized on locally refined meshes to increase the

accu-racy of our numerical solutions. Our FEM simulation for predicting the center frequency of the stopping bands of the

reso-nator is convergent and accurate with an error of approximately 1% compared with experimental data, as shown inTables 2

and 3. For computing the dispersion diagram near the center of the stopping band of the resonator, we transform the GEP

(12)into the TPQEP(3)to reveal the important reciprocal relationship of the eigenvalues in which the eigenvalues appear in

2.14 2.16 2.18 2.2 2.22 2.24 2.26 2.28 2.3 2.32 x 109 0 1 2 3 4 5 6 7 x 10 −15 frequencyω

max. relative residual

Fig. 6. Maximal relative residuals of the 10 desired eigenvalues produced by Algorithm 1.

2.14 2.16 2.18 2.2 2.22 2.24 2.26 2.28 2.3 2.32 x 109 470 480 490 500 510 520 530 540 550 frequencyω

CPU times (sec.)

Arnoldi Algorithm 1

reciprocal pairs ðk; 1=kÞ. The TPQEP(3)is then solved by GTSHIRA so that the reciprocal relationship of the eigenvalues can be automatically preserved. Our numerical results show that the traditional Arnoldi method converges slowly and fails to pre-serve the reciprocal property of the eigenvalues near the unit circle. In constrast, the proposed structure-preserving method in Algorithm 1 not only converges to those eigenpairs faster than does the traditional Arnoldi method, but it also possesses high accuracy in terms of the relative residual error. Furthermore, the reciprocal property of the eigenpairs is kept nicely under machine precision. Recently, some studies have been shown that the wave propagation speed in PZT can be affected

by stress[26,27]and temperature[28]. With the proposed structure-preserving algorithm, we would like to investigate how

the center and range of the SAW frequency changes under initial stress and thermal effect in our future studies. Acknowledgments

This work is partially supported by the National Science Council and the National Center for Theoretical Sciences in Tai-wan. The author C.T. Wu thanks the support from National Science Council in Taiwan under the Grant No. 99-2115-M-009-001.

References

[1]H. Allik, T. Hughes, Finite element method for piezoelectric vibration, Int. J. Numer. Methods Eng. 2 (1970) 151–157.

[2]M.B. Angel, M.I. Rocha-Gaso, M.I. Carmen, A.V. Antonio, Surface generated acoustic wave biosensors for detection of pathogens: a review, Sensors 9 (2009) 5740–5769.

[3]S.N. Atluri, M. Nakagaki, Computational methods for plane problems of fracture, in: Computational Method in Mechanics of Fracture, North-Holland Publishing Co., Amsterdam, 1986, pp. 170–227.

[4]B.A. Auld, Acoustic Fields and Waves in Solids, Krieger Publishing Company, Malabar, Florida, 1990. [5]W. Bond, The mathematics of the physical properties of crystals, BSTJ 22 (1943) 1–72.

[6] M. Buchner, W. Ruile, A. Dietz, R. Dill, FEM analysis of the reflection coefficient of SAWS in an infinite periodic array, in: Proceedings of the IEEE Ultrasonics Symposium, 1991, pp. 371–375.

[7] C. Cai, H. Zheng, M.S. Khan, K.C. Hung, Modeling of material damping properties in ANSYS, in: International ANSYS Conference Proceedings, 2002. [8]C.K. Campbell, Surface Acoustic Wave Devices for Mobile and Wireless Communications, Academic Press, INC., 1998.

[9]D.P. Chen, H.A. Haus, Analysis of metal-strip SAW gratings and transducers, IEEE Trans. Sonics Ultrason. 32 (3) (1985) 395–408.

[10]L.C. Chin, V.V. Vardan, V.K. Vardan, Hybrid finite element formulation for periodic piezoelectric arrays subjected to fluid loading, Int. J. Numer. Methods Eng. 37 (1994) 2987–3003.

[11]E.K.-W. Chu, T.-M. Hwang, W.-W. Lin, C.-T. Wu, Vibration of fast trains, palindromic eigenvalue problems and structure-preserving doubling algorithms, J. Comput. Appl. Math. 219 (2008) 237–252.

[12]P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman Publishing INC., 1985.

[13]M. Hofer, N. Finger, G. Kovacs, J. Schöberl, S. Zaglmayr, U. Langer, R. Lerch, Finite-element simulation of wave propagation in periodic piezoelectric SAW structure, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 53 (6) (2006) 1192–1201.

[14]T.-M. Huang, W.-W. Lin, J. Qian, Structure-preserving algorithms for palindromic quadratic eigenvalue problems arising from vibration on fast trains, SIAM J. Matrix Anal. Appl. 30 (2008) 1566–1592.

[15]C.F. Ipsen, Accurate eigenvalues for fast trains, SIAM News 37 (2004).

[16]M. Koshiba, S. Mitobe, M. Suzuki, Finite-element solution of periodic waveguides for acoustic waves, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 34 (4) (1987) 472–477.

[17]J. Kushibiki, I. Takanaga, T. Sannomiya, Accurate measurements of the acoustical physical constants of LiNbO3and LiTaO3single crystals, IEEE Trans.

Ultrason. Ferroelectr. Freq. Control 46 (1999) 1315–1323.

[18]R. Lerch, Simulation of piezoelectric devices by two- and three-dimensional finite elements, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 37 (2) (1990) 1990.

[19]W.-W. Lin, A new method for computing the closed-loop eigenvalues of a discrete-time algebraic Riccatic equation, Linear Algebra Appl. 96 (1987) 157–180.

[20]J.E. San Miguel Medina, F. Buiochi, J.C. Adamowski, Numerical modeling of a circular piezoelectric ultrasonic transducer radiating in water, ABCM Symp. Ser. Mechatron. 2 (2006) 458–464.

[21]S. Mitobe, M. Koshiba, M. Suzuki, Theoretical determination of equivalent circuit parameters for reflective SAW metallic gratings, Electron. Commun. Jpn. 70 (6) (1987) 37–45.

[22]M. Mohamed, EL Gowini, W.A. Moussa, A finite element model of a MEMS-based surface acoustic wave hydrogen sensor, Sensors 10 (2010) 1232– 1250.

[23]G. Nadar, E.C. Nelli Silva, J.C. Adamowski, Effective damping value of piezoelectric transducer determined by experimental techniques and numerical analysis, ABCM Symp. Ser. Mechatron. 1 (2004) 271–279.

[24]R.V. Patel, On computing the eigenvalues of a symplectic pencil, Linear Algebra Appl. 188 (1993) 591–611.

[25] U. Rösler, D. Cohrs, A. Dietz, G. Fisherauer, W. Ruile, P. Russer, R. Weigel. Determination of leaky saw propagation, reflection and coupling on LiTaO3, in:

IEEE Ultrasonics Symposium, 1995, pp. 247–250.

[26]A.M. Abd-Alla, S.M. Ahmed, Propagation of love waves in a non-homogeneous orthotropic elastic layer under initial stress overlying semi-infinite medium, Appl. Math. Comput. 106 (1999) 265–275.

[27]M.S. Son, Y.J. Kang, The effect of initial stress on the propagation behavior of SH waves in piezoelectric coupled plates, Ultrasonics 51 (4) (2011) 489– 495.