• 沒有找到結果。

# Jacobi Davidson method and its applications

N/A
N/A
Protected

Academic year: 2022

Share "Jacobi Davidson method and its applications"

Copied!
74
0
0

(1)

### Jacobi Davidson method and its applications

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University, Taiwan

October 31, 2008

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 1 / 71

(2)

### Outline

1 Some basic theorems

2 Jacobi-Davidson method for linear eigenvalue problems

3 Polynomial eigenvalue problems

4 Locking method for computing λ2, . . . , λ`

5 Non-equivalence deflation of quadratic eigenproblems

6 Non-equivalence deflation for cubic polynomial eigenproblems

7 Applications

8 Numerical experiments for linear eigenproblems

9 Numerical experiments for polynomial eigenproblems

(3)

### Some basic theorems

Eigenproblems

Standard eigenproblems: Ax = λx Generalized eigenproblems: Ax = λBx

Higher order poly. eigenproblems: (A0+ λA1+ .... + λnAn)x = 0 Eigenproblems of λ-matrices: F (λ)x = 0

What do we care ?

(i) In theory: eigenstructure, spectrum decomposition, canonical form, . . . , etc.

(ii) In computation: eigenvalues, eigenvectors, invariant subspaces, . . . , etc.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 3 / 71

(4)

Theorem (Fischer)

Let the Hermitian matrix A have eigenvalues λ1 ≥ λ2 ≥ · · · ≥ λn. Then

λi = min dim(W)=n−i +1



w ∈W,kw kmax2=1wHAw



and

λi = max dim(W)=i



min

w ∈W,kw k2=1wHAw

 .

(5)

Theorem

Let X be an eigenspace of A and let X be a basis for X . Then there is a unique matrix L such that

AX = XL.

The matrix L is given by

L = XIAX , where XI is a matrix satisfying XIX = I .

If (λ, x ) is an eigenpair of A with x ∈ X , then (λ, XIx ) is an eigenpair of L.

Conversely, if (λ, u) is an eigenpair of L, then (λ, Xu) is an eigenpair of A.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 5 / 71

(6)

### Proof:

Let

X = [x1· · · xk] and Y = AX = [y1· · · yk] .

Since yi ∈ X and X is a basis for X , there is a unique vector `i such that yi = X `i.

If we set L = [`1· · · `k], then AX = XL and L = XIXL = XIAX .

Now let (λ, x ) be an eigenpair of A with x ∈ X . Then there is a unique vector u such that x = Xu. However, u = XIx . Hence

λx = Ax = AXu = XLu ⇒ λu = λXIx = Lu.

Conversely, if Lu = λu, then

A(Xu) = (AX )u = (XL)u = X (Lu) = λ(Xu), so that (λ, Xu) is an eigenpair of A.

(7)

Theorem (Optimal residuals) Let [X X] be unitary. Let

R = AX − XL and SH = XHA − LXH. Then kRk and kS k are minimized when

L = XHAX , in which case

(a) kRk = kXHAX k, (b) kSk = kXHAXk, (c) XHR = 0.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 7 / 71

(8)

### Proof:

Set

 XH XH

 A

X X  =

 ˆL H

G M

 . Then

 XH XH

 R =

 ˆL H

G M

  XH XH

 X −

 XH XH

 XL =

 L − Lˆ G

 . It implies that

kRk =

 XH XH

 R

=

 ˆL − L G

 , which is minimized when L = ˆL = XHAX and

min kRk = kG k = kXHAX k.

The proof for S is similar. If L = XHAX , then

XHR = XHAX − XHXL = XHAX − L = 0.

(9)

Definition

Let X be of full column rank and let XI be a left inverse of X . Then XIAX is a Rayleigh quotient of A.

Theorem

Let X be orthonormal, A be Hermitian and R = AX − XL.

If `1, . . . , `k are the eigenvalues of L, then there are eigenvalues λj1, . . . , λjk of A such that

|`i − λji| ≤ kRk2 and v u u t

k

X

i =1

(`i − λji)2 ≤√

2kRkF.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 9 / 71

(10)

### Jacobi’s orthogonal component correction (JOCC), 1846

Consider the eigenvalue problem A

 1 z



 α cT b F

  1 z



= λ

 1 z



, (1)

where A is diagonal dominant and α is the largest diagonal element.

(1) is equivalent to

 λ = α + cTz, (F − λI )z = −b.

Jacobi iteration : (with z1 = 0)

 θk = α + cTzk,

(D − θkI )zk+1 = (D − F )zk − b (2) where D = diag(F ).

(11)

### Davidson’s method (1975)

Algorithm (Davidson’s method) Given unit vector v , set V = [v ] Iterate until convergence

Compute desired eigenpair (θ, s) ofVTAV. Computeu = Vs andr = Au − θu.

If (k r k2< ε), stop.

Solve (DA− θI )t = r.

Orthog. t ⊥ V → v , V = [V , v ] end

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 11 / 71

(12)

Let uk = (1, zkT)T. Then

rk = (A − θkI )uk =

 α − θk+ cTzk (F − θkI )zk + b



Substituting the residual vector rk into linear systems (DA− θkI )tk = −rk, where DA =

 α 0

0 D

 , we get

(D − θkI )yk = −(F − θkI )zk − b

= (D − F )zk − (D − θkI )zk− b From (2) and above equality, we see that

(D − θkI )(zk+ yk) = (D − F )zk − b = (D − θkI )zk+1

This implies that zk+1 = zk+ yk as one step of JOCC starting with zk.

(13)

### Polynomial eigenvalue problems

Polynomial eigenvalue problems:

A(λ)x ≡ λτAτ + · · · + +λ2A2+ λA1+ A0 x = 0. (3) Enlarged linear eigenvalue problem: (e.g. cubic polynomial)

0 I 0

0 0 I

A0 A1 A2

 x λx λ2x

= λ

I 0 0

0 I 0

0 0 −A3

 x λx λ2x

. Disadvantages:

I The order of the larger matrices are tripled

I The condition number of eigenvalues and eigenvectors may increase

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 13 / 71

(14)

Consider the quadratic eigenvalue problem1 Q(λ)x ≡ λ2M + λC + K x = 0 with

M = 1

2In, K = diag1≤j ≤n j2π2(j2π2+ τ − κv2)/2 , C = −CT = (cij) with cij =

( 4ij

j2−i2v if i + j is odd 0 otherwise.

Take v = 10, κ = 0.8 and τ = 77.9.

Enlarged linear eigenvalue problem:

 K 0 0 I

  x λx



= λ

 −C −M

I 0

  x λx



1F. Tisseur and K. Meerbergen, SIAM Rev. Vol. 43, pp. 235-286, 2001

(15)

−1 0 1

x 10−9

−4

−3

−2

−1 0 1 2 3

4x 104 n = 60

polyeig exact sol

Figure: The spectrum of Q(λ)

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 15 / 71

(16)

### Polynomial Jacobi-Davidson method

k, uk): approx. eigenpair of A(λ), θk ≈ λ, with

uk = Vksk, VkTA(λ)Vksk = 0 and k sk k2= 1.

Let

rk = A(θk)uk. Then

ukTrk = uTkA(θk)uk = skTVkTA(θk)Vksk = 0 ⇒ rk ⊥ uk Find the correction t such that

A(λ)(uk+ t) = 0.

That is

A(λ)t = −A(λ)uk = −rk + (A(θk) − A(λ))uk.

(17)

Let

pk = A0k)uk

τ

X

i =1

i θi −1k Ai

! uk.

A(λ) = A − λI :

pk = −uk,

(A(θk) − A(λ))uk = (λ − θk)uk = (θk − λk)pk A(λ) = A − λB:

pk = −Buk,

(A(θk) − A(λ))uk = (λ − θk)Buk = (θk− λ)pk A(λ) =Pτ

i =0λiAi with τ ≥ 2:

(A(θk) − A(λ))uk =



k− λ)A0k) −1

2(θk − λ)2A00k)

 uk

= (θk − λ)pk −1

2(θk − λ)2A00k)uk

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 17 / 71

(18)

Hence

A(λ)t = −rk+ (θk − λ)pk− 1

2(θk− λ)2A00k)uk

Since rk ⊥ uk, we have



I −pkukT ukTpk



A(λ)t = −rk −1

2(θk− λ)2



I − pkukT ukTpk



A00k)uk. Correction equation:



I − pkuTk ukTpk



A(θk)(I − ukukT)t = −rk and t ⊥ uk,

or



I − pkukT ukTpk



(A − θkB)



I − ukpTk pkTuk



t = −rk and t ⊥B uk,

with symmetric positive definite matrix B.

(19)

### Solving correction vector t

Correction equation:



I − pkuTk ukTpk



A(θk)(I − ukukT)t = −rk. (4) Method I:

Use preconditioning iterative approximations, e.g., GMRES, to solve (4).

Use a preconditioner Mp



I − pkuTk ukTpk

 M

I − ukuTk





I − pkukT ukTpk

 A(θk)



I − ukukT

 , where M is an approximation of A(θk).

In each of the iterative steps, it needs to solve the linear system

Mpt = y , t ⊥ uk (5)

for a given y .

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 19 / 71

(20)

Since t ⊥ uk, Eq. (5) can be rewritten as



I − pkuTk ukTpk



Mt = y ⇒ Mt = ukTMt

ukTpk pk+ y ≡ ηkpk + y . Hence

t = M−1y + ηkM−1pk, where

ηk = − ukTM−1y ukTM−1pk.

SSOR preconditioner: Let A(θk) = L + D + U. Then M = (D + ωL)D−1(D + ωU).

(21)

Method II: Since t ⊥ uk, Eq. (4) can be rewritten as A(θk)t = uTkA(θk)t

ukTpk pk − rk ≡ εpk − rk. (6) Let t1 and t2 be approximated solutions of the following linear

systems:

A(θk)t = −rk and A(θk)t = pk, respectively. Then the approximated solution ˜t for (6) is

˜t = t1+ εt2 for ε = −ukTt1 ukTt2

. The approximated solution ˜t for (6) is

˜t = −M−1rk+ εM−1pk for ε = uTkM−1rk ukTM−1pk, where M is an approximation of A(θk).

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 21 / 71

(22)

Method III:

Eq. (6) implies that

t = εA(θk)−1pk − A(θk)−1rk = εA(θk)−1pk − uk. Let t1 be approximated solution of the following linear system:

A(θk)t = pk. Then the approximated solution ˜t for (6) is

˜t = εt1− uk for ε =

ukTt1−1

.

(23)

Algorithm (Jacobi-Davidson Algorithm for solving A(λ)x = 0) Choose an n-by-m orthonormal matrix V0

Do k = 0, 1, 2, · · ·

Compute all the eigenpairs of VkTA(λ)Vk = 0.

Select the desired (target) eigenpair (θk, sk) with kskk2= 1.

Computeuk = Vksk,rk = A(θk)uk andpk = A0k)uk. If (krkk2 < ε), λ = θk, x = uk, Stop

Solve (approximately) a tk ⊥ uk from (I − pku

T k

uTkpk)A(θk)(I − ukuTk)t = −rk. Orthogonalize tk ⊥ Vk → vk+1, Vk+1 = [Vk, vk+1]

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 23 / 71

(24)

### Restarting

A double precision real variable needs to 8 bytes to save in Memory.

125 double precision real variables ≈ 1 KB 125, 000 double precision real variables ≈ 1 MB

Keep the locked Schur vectors as well as the Schur vectors of interest in the subspace and throw away those we are not interested.

(25)

(v) Solve correct equation (approximately) to obtain a t ⊥ uk by the method determined below.

If ( krkk2> 0.1 and k ≤ 9 ) then

Use {BiCGSTAB, No precond., 7, 10−3} else

Use {GMRES, SSOR, 30, 10−3} End if

Figure: The heuristic strategy for computing the first target eigenvalue.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 25 / 71

(26)

### Locking for Ax = λx

Vk with VkV = Ik are convergentSchur vectors, i.e., AVk = VkTk

for someupper triangular Tk. SetV = [Vk, Vq]with VV = Ik+q in k + 1-th iteration of Jacobi-Davidson Algorithm. Then

VAV =

 VkAVk VkAVq VqAVk VqAVq



=

 Tk VkAVq VqVkTk VqAVq



=

 Tk VkAVq 0 VqAVq

 .

(27)

Locking Polynomial Jacobi–Davidson Method (0) Given A(λ) =Pτ

i =0λiAi and the number of desired eigenvalues σ.

(1) Initialize V = [Vini] as an orthonormal matrix and Vx = [].

(2) For j = 1, 2, ..., σ

(2.1) Iterate until convergence

(i) Compute thej th desired eigenpairs (θ, s) of VTA(θ)V . (ii) Compute u = Vs, p = A0(θ)u, and r = A(θ)u.

(iii) If (||r ||2 < ε), Set λj = θ, xj = u, Goto Step (2.2).

(iv) Solve (approximately) a t ⊥ u from (I − puuTTp)A(θ)(I − uuT)t = −r . (v) Orthogonalize t ⊥ V → v , V = [V , v ] (2.2) Orthogonalizexj ⊥ Vx → xj; Vx = [Vx, xj]

(2.3) Choose an orthonormal matrixVini ⊥ Vx; Set V = [Vx, Vini] Ref: G. L. G. Sleijpen, G. L. Booten, D. R. Fokkema and H. A. van der Vorst, BIT, 36:595-633, 1996

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 27 / 71

(28)

(v) Solve correct equation (approximately) to obtain a t ⊥ uk by the method determined below.

If ( k rk k2> 0.1 and k < 10 ) then

Use {BiCGSTAB, No precond., 7, 10−3} else if ( k rk k2≥ 0.1 and k > 14 ) then

Use {GMRES, SSOR, 30, 10−3}

else if ( k rk k2< 0.1 and k rk−1 k2/ k rk k2< 4) then Set j = min(30, j + 2) and use {GMRES, SSOR, j, 10−3} else

Use {GMRES, SSOR, j, 10−3} End if

Figure: The heuristic strategy for computing eigenvalues other than the first convergent eigenvalue.

(29)

### Non-equivalence deflation of quadratic eigenproblems

Let λ1 be a real eigenvalue of Q(λ) and x1, z1 be the associated right and left eigenvectors, respectively, with z1TKx1 = 1. Let

θ1 = (z1TMx1)−1. We introduce a deflated quadratic eigenproblem

Q(λ)x ≡e h

λ2M + λ ee C + eKi x = 0, where

Me = M − θ1Mx1z1TM, Ce = C + θ1

λ1

(Mx1z1TK + Kx1z1TM), Ke = K − θ1

λ21Kx1z1TK .

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 29 / 71

(30)

### Complex deflation

Let λ1 = α1+ i β1 be a complex eigenvalue of Q(λ) and x1 = x1R+ ix1I, z1 = z1R+ iz1I be the associated right and left eigenvectors, respectively, such that

Z1TKX1 = I2, where X1 = [x1R, x1I] and Z1 = [z1R, z1I]. Let

Θ1 = (Z1TMX1)−1.

Then we introduce a deflated quadratic eigenproblem with Me = M − MX1Θ1Z1TM,

Ce = C + MX1Θ1Λ−T1 Z1TK + KX1Λ−11 ΘT1Z1TM, Ke = K − KX1Λ−11 Θ1Λ−T1 Z1TK

in which Λ1 =

 α1 β1

−β1 α1

 .

(31)

Theorem

(i) Let λ1 be a simple real eigenvalue of Q(λ). Then the spectrum of Q(λ) is given bye

(σ(Q(λ)){λ1}) ∪ {∞}

provided that λ21 6= θ1.

(ii) Let λ1 be a simple complex eigenvalue of Q(λ). Then the spectrum of eQ(λ) is given by

σ(Q(λ)){λ1, ¯λ1} ∪ {∞, ∞}

provided that Λ1ΛT1 6= Θ1.

Furthermore, in both cases (i) and (ii), if λ2 6= λ1 and (λ2, x2) is an eigenpair of Q(λ) then the pair (λ2, x2) is also an eigenpair of eQ(λ).

Ref: T.-M. Hwang, W.-W. Lin and V. Mehrmann, SIAM J. Sci. Comput.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 31 / 71

(32)

Suppose that M, C , K are symmetric. Given an eigenmatrix pair

1, X1) ∈ Rr ×r × Rn×r of Q(λ), where Λ1 is nonsingular and X1 satisfies X1TKX1 = Ir, Θ1:= (X1TMX1)−1.

We define ˜Q(λ) := λ2M + λ ˜˜ C + ˜K , where M˜ := M − MX1Θ1X1TM,

C˜ := C + MX1Θ1Λ−T1 X1TK + KX1Λ−11 Θ1X1TM, K˜ := K − KX1Λ−11 Θ1Λ−T1 X1TK .

Theorem

Suppose that Θ1− Λ1ΛT1 is nonsingular. Then the eigenvalues of the real symmetric quadratic pencil ˜Q(λ) are the same as those of Q(λ) except that the eigenvalues of Λ1, which are closed under complex conjugation, are replaced by r infinities.

(33)

### Proof:

Since (Λ1, X1) is an eigenmatrix pair of Q(λ), i.e., MX1Λ21+ CX1Λ1+ KX1 = 0, we have

Q(λ)˜ = Q(λ) + [MX1(λIr + Λ1) + CX1] Θ1Λ−T1 (X1TK − λΛT1X1TM)

= Q(λ) + Q(λ)X1(λIr − Λ1)−1Θ1Λ−T1 (X1TK − λΛT1X1TM).

By using the identity

det(In+ RS ) = det(Im+ SR), where R, ST ∈ Rn×m, we have

det[ ˜Q(λ)]

= det[Q(λ)]det[I + X1(λIr − Λ1)−1Θ1Λ−T1 (X1TK − λΛT1X1TM)]

= det[Q(λ)]det[Ir + (λIr − Λ1)−1Θ1Λ−T1 (Ir − λΛT1Θ−11 )]

= det[Q(λ)]

det(λIr − Λ1)det(Θ1Λ−T1 − Λ1).

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 33 / 71

(34)

Since (Θ1− Λ1ΛT1) ∈ Rr ×r is nonsingular, we have det(Θ1Λ−T1 − Λ1) 6= 0.

Therefore, ˜Q(λ) has the same eigenvalues as Q(λ) except that r eigenvalues of Λ1 are replaced by r infinities.

(35)

### Non-equivalence deflation for cubic polynomial eigenproblems

Let (Λ, Vu) ∈ Rr ×r × Rn×r be an eigenmatrix pair of A(λ) with VuTVu= Ir and 0 /∈ σ(Λ), i.e.,

A3VuΛ3+ A2VuΛ2+ A1VuΛ + A0Vu= 0. (7) Define a new deflated cubic eigenvalue problem by

A(λ)u = (λe 33+ λ22+ λ˜A1+ ˜A0)u = 0, (8) where





0 = A0,

1 = A1− (A1VuVuT+ A2VuΛVuT + A3VuΛ2VuT), A˜2 = A2− (A2VuVuT+ A3VuΛVuT),

3 = A3− A3VuVuT.

(9)

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 35 / 71

(36)

Lemma

Let A(λ) and eA(λ) be cubic pencils given by (3) and (8), respectively.

Then it holds

A(λ) = A(λ)e 

In− λVu(λIr − Λ)−1VuT

. (10)

Theorem

Let (Λ, Vu) be an eigenmatrix pair of A(λ) with VuTVu= Ir. Then (i) (σ(A(λ))σ(Λ)) ∪ {∞} = σ(eA(λ)).

(ii) Let (µ, z) be an eigenpair of A(λ) with k z k2= 1 and µ /∈ σ(Λ).

Define

˜

z = (In− µVuΛ−1VuT)z ≡ T (µ)z. (11) Then (µ, ˜z) is an eigenpair of eA(λ).

(37)

### Proof of Lemma:

Using (9) and (7), and the fundamental matrix calculation, we have A(λ)e = A(λ) − λ



λ2A3VFVFT + λA2VFVFT + λA3VFΛVFT+ A1VFVFT +A2VFΛVFT + A3VFΛ2VFT



= A(λ) − λ

A3VF(λIr − Λ)3(λIr − Λ)−1VFT +3A3VFΛ(λIr − Λ)2(λIr − Λ)−1VFT +3A3VFΛ2(λIr − Λ)(λIr − Λ)−1VFT +A2VF(λIr − Λ)2(λIr − Λ)−1VFT +2A2VFΛ(λIr − Λ)(λIr − Λ)−1VFT

+A1VF(λIr − Λ)(λIr − Λ)−1VFT

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 37 / 71

(38)

A(λ)e = A(λ) − λA3VF3Ir − Λ3) + A2VF2Ir − Λ2) +A1VF(λIr − Λ) + A0VF − A0VF] (λIr − Λ)−1VFT

o

= A(λ) − λh

A(λ)VF(λIr − Λ)−1VFTi

= A(λ)h

In− λVF(λIr − Λ)−1VFTi .

(39)

### Proof of Theorem

: (i) Using the identity

det(In+ RS ) = det(Im+ SR) and Lemma 8, we have

det( eA(λ)) = det(A(λ)) det



In− λVF(λIr − Λ)−1VFT



= det(A(λ)) det In− λ(λIr − Λ)−1

= det(A(λ)) det(λIr − Λ)−1det(−Λ).

Since 0 /∈ σ(Λ), det(−Λ) 6= 0. Thus, eA(λ) and A(λ) have the same finite spectrum except the eigenvalues in σ(Λ). Furthermore, dividing Eq. (8) by λ3 and using the fact that

Ae3VF = (A3− A3VFVFT)VF = 0,

we see that (diagr{∞, · · · , ∞}, VF) is an eigenmatrix pair of eA(λ) corresponding to infinite eigenvalues.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 39 / 71

(40)

(ii) Since µ /∈ σ(Λ), the matrix T (µ) = (I − µVFΛ−1VFT) in (11) is invertible with the inverse

T (µ)−1= In− µVF(µIr − Λ)−1VFT. (12) From Lemma 8, we have

A(µ)˜e z = A(µ) h

In− µVF(µIr − Λ)−1VFT i h

In− µVFΛ−1VFT i

z = 0.

This completes the proof.

(41)

### Applications

Quantum well (2 dim.)

Quantum wire (1 dim.)

Quantum dot (0 dim.)

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 41 / 71

(42)

### Nanometer

1nm = 10−9m

Nano-scale ≈ 1–100 nm A semiconductor QD ≈ 10 nm QD : hair ≈ 1 : 10000

Why consider quantum effects?

Small devices imply significant quantum effect

(43)

### Quantum dot

Cross-sections of hetero-structure InAs/GaAs QDs by Transmission Electron Microscope [Schoenfeld, 00]

150 Å

300 Å 150 Å

300 Å 150 Å

300 Å 150 Å

300 Å

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 43 / 71

(44)

(45)

### Quantum dot

Cross-sections of hetero-structure InAs/GaAs QDs by Transmission Electron Microscope [Schoenfeld, 00]

150 Å

300 Å 150 Å

300 Å 150 Å

300 Å 150 Å

300 Å

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 45 / 71

(46)

### Energy levels (eigenvalues)

Confined Discrete Energy Levels Confined Discrete

Energy Levels

(47)

### Numerical experiments for linear eigenproblems

TheSchr¨odinger equationfor Semiconductor:

−∇ · (α∇u) + Vu = λu,

where α =

( α2m~2

1 inside, α+2m~2

2 outside, V =

 V= V1 inside, V+= V2 outside

~: Plank constant m`: parabolic effective mass V`: confinement potential λ: total energy

Interface condition:

α∂u

∂n ∂D

= α+∂u

∂n ∂D+

Dirichlet boundary conditions

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 47 / 71

(48)

Symmetric eigenvalue problem

Ax = λx where A is a symmetric matrix.

Reference:

Tsung-Min Hwang, Wen-Wei Lin, Wei-Cheng Wang and Weichung Wang, Numerical simulation of three dimensional pyramid quantum dot, Journal of Computational Physics, Vol 196, pp. 208-232, 2004.

Figure: PRB, 54, 8743, (1996)

(49)

### Three dimensional pyramid quantum dot

H

W W

Pyramid Dot Cuboid Matrix i

2 1. 5 1 0. 5 0 0.5 1 1.5 2

1 0. 5 0 0.5 1 1.5 2

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 2 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 2 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 2 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 1 2 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 1 1 2 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 1 1 1 5 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 1 1 3 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 1 3 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 3 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 3 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 3 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 49 / 71

(50)

Finite volume discretized scheme Symmetric eigenvalue problems:

Ax = λx

Second order convergent rate Correction vector with SSOR preconditioner:

t = −MA−1r + εMA−1p

with

ε = uTMA−1r uTMA−1p,

whereMA= (D + ωL)D−1(D + ωU).

H

W W

Pyramid Dot Cuboid Matrix i

Figure: Structure schema of a pyramid quantum dot.

(51)

The width of the QD (matrix) base is 12.4 nm (24.8 nm); the height of the QD (matrix) is 6.2 nm (18.6 nm).

InAs QD: α1 = 0.024me and V1 = 0.0 GaAs matrix: α2= 0.067me and V2= 0.70.

Stopping criteria: residual < 10−10 Convergent rate

(L,M,N) Mtx. dim. λ1 Rate λ2 Rate

(16, 16, 12) 2,475 0.4226 - 0.6527 - (32, 32, 24) 22,103 0.4001 - 0.6423 - (64, 64, 48) 186,543 0.3934 1.744 0.6391 1.708 (128,128, 96) 1,532,255 0.3916 1.905 0.6383 1.866 (256,256,192) 12,419,775 0.3911 1.954 0.6380 1.912 Table: Convergent rate =log2 (4h)− λ(2h))/(λ(2h)− λ(h))

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 51 / 71

(52)

dim(A) = 1, 532, 255on PC with P4 1.8 GHz CPUand1 GB of main memory.

Value Ite. no. CPU time (sec.)

λ1 0.3916 72 278.0

λ2 0.6383 72 284.3

λ3 0.6383 125 521.4

dim(A) = 32, 401, 863on HP workstation with a 1.0 GHz Intel Itanium II CPUand 12 GBsof main memory.

Value Ite. no. CPU time (sec.)

λ1 0.3910 138 5,852

λ2 0.6380 133 5,354

λ3 0.6380 220 8,511

(53)

Unsymmetric eigenvalue problem

Ax = λx where A is a unsymmetric matrix.

Reference:

Tsung-Min Hwang, Wei-Hua Wang and Weichung Wang, Efficient numerical schemes for electronic states in coupled quantum dots, accepted for publication in Journal of Nanoscience and

Nanotechnology.

JAP, 90-12, (2001)-

Figure: JAP, 90-12, (2001)

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 53 / 71

(54)

### Vertically aligned quantum dot array

(a) Nonuniform mes h (b) Uniform mesh

(55)

Reduce three-dimensional systems to two-dimensional systems by Fourier transformation.

Finite volume discretized scheme Unsymmetric eigenvalue problems:

Ax = λx

Second order convergent rate

Method II with SSOR preconditioner:

MA= (D + ωL)D−1(D + ωU).

Figure: Structure schema of a cylindrical vertically aligned quantum dot array.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 55 / 71

(56)

Reduce three-dimensional systems to two-dimensional systems by Fourier transformation.

Finite volume discretized scheme Unsymmetric eigenvalue problems:

Ax = λx

Second order convergent rate

Method II with SSOR preconditioner:

MA= (D + ωL)D−1(D + ωU). Figure: Structure schema of a cylindrical vertically aligned quantum dot array.

(57)

Reduce three-dimensional systems to two-dimensional systems by Fourier transformation.

Finite volume discretized scheme Unsymmetric eigenvalue problems:

Ax = λx

Second order convergent rate

Method II with SSOR preconditioner:

MA= (D + ωL)D−1(D + ωU).

Figure: Structure schema of a cylindrical vertically aligned quantum dot array.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 55 / 71

(58)

gap = 6nm and dim(A) = 12, 288, 000on HP workstation with a 1.0 GHz Intel Itanium II CPU and12 GBsof main memory.

r2 = 3.198 r2 = 3.223

Value Ite. no. Time (sec.) Value Ite. no. Time (sec.)

λ1 0.1587 83 20,287 0.1587 84 23,840

λ2 0.35553 109 26,030 0.3526 91 31,612

λ3 0.3558 53 13,831 0.3555 61 25,706

gap = 3nm and dim(A) = 11, 059, 200on HP workstation with a 1.0 GHz Intel Itanium II CPU and12 GBsof main memory.

r2 = 3.198 r2 = 3.223

Value Ite. no. Time (sec.) Value Ite. no. Time (sec.)

λ1 0.1586 85 20,072 0.1586 83 18,455

λ2 0.3540 257 56,811 0.3515 130 29,506

λ3 0.3564 53 12,257 0.3558 54 13,082

(59)

Generalized eigenvalue problem

Ax = λBx

where A is symmetric positive definiteand B is a positive diagonalmatrix.

Reference:

Tsung-Min Hwang, Wei-Cheng Wang and Weichung Wang, Numerical schemes for three dimensional irregular shape quantum dots over curvilinear coordinate systems, accepted for publication in Journal of Computational Physics.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 57 / 71

(60)

### Three dimensional arbitrary shape quantum dots

Quantum dot Matrix

ZA ZB

ZC

1 M

Rmtx

0 H

N1

N2

N

H HC

HRmtx

R adial coordinate

Axial coordinate

(61)

Curvilinear coordinate system

Reduce three-dimensional systems to two-dimensional systems by Fourier transformation.

Jump condition capturing scheme (nine points)

Generalized eigenvalue problems:

Ax = λBx , A issymmetric positive definite and B is apositive diagonal matrix.

Second order convergent rate

Method I with SSOR preconditioner:

MA= (D + ωL)D−1(D + ωU).

Quantum dot Matrix

Figure: Structure schema of the quantum dot model.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 59 / 71

Matrix factorization (MF) and its extensions are now widely used in recommender systems.. In this talk I will briefly discuss three research works related to

BAL 1000 Brown almost-linear func, nonconvex, dense Hessian1. BT 1000 Broyden tridiagonal func, nonconvex,

Wang, Unique continuation for the elasticity sys- tem and a counterexample for second order elliptic systems, Harmonic Analysis, Partial Differential Equations, Complex Analysis,

Polynomial Jacobi Davidson Method for Large/Sparse Eigenvalue Problems..

Department of Mathematics – NTNU Tsung-Min Hwang November 16, 2003... We say that the rate of

Arbenz et al.[1] proposed a hybrid preconditioner combining a hierarchical basis preconditioner and an algebraic multigrid preconditioner for the correc- tion equation in the

We give some numerical results to illustrate that the first pass of Algorithm RRLU(r) fails but the second pass succeeds in revealing the nearly rank

We will give a quasi-spectral characterization of a connected bipartite weighted 2-punctually distance-regular graph whose halved graphs are distance-regular.. In the case the