Luận văn Thạc sĩ Toán học: Bài toán giá trị riêng bậc hai

Mục lục

Mở đầu 1

Danh sách ký hiệu 2

1 Bài toán giá trị riêng bậc hai 3

1.1 Bài toán giá trị riêng tiêu chuẩn . . . . . . . . . 3 . . . . . . .

1.1.1 Khái niệm . . . . . . . . . . . . . . . . 3 . . . . . . .

1.1.2 Một số thuật toán tìm giá trị riêng . . . 4 . . . . . . .

1.1.3 Bài toán giá trị riêng suy rộng . . . . . 9 . . . . . . .

1.2 Bài toán giá trị riêng bậc hai . . . . . . . . . . 13 . . . . . . .

1.2.1 Khái niệm . . . . . . . . . . . . . . . . 13 . . . . . . .

1.2.2 Tuyến tính hóa bài toán giá trị riêng bậc hai 15 . . . . .

1.2.3 Bộ ba Jordan của Q(λ ) . . . . . . . . . . . 18 . . . . .

1.2.4 Một số tính chất của bài toán giá trị riêng bậc hai 19 . .

1.3 Một số ứng dụng khác của bài toán giá trị riêng bậc hai 21 . . .

1.3.1 Biểu diễn nghiệm của phương trình vi phân tuyến

tính cấp hai . . . . . . . . . . . . . . . 21 . . . . . . .

1.3.2 Bài toán hạn chế bình phương nhỏ nhất 22 . . . . . . .

1.3.3 Một vài ví dụ . . . . . . . . . . . . . . 23 . . . . . . .

2 Giải số bài toán giá trị riêng bậc hai 26

2.1 Phương pháp số cho bài toán đặc . . . . . . . . 26 . . . . . . .

2.1.1 Phương pháp Newton . . . . . . . . . . 26 . . . . . . .

2.1.2 Phân tích Schur thực suy rộng . . . . . . . . . . . . 28

2.2 Phương pháp số cho bài toán thưa . . . . . . . . . . . . . . 29

2.3 Ví dụ số . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

Kết luận 35

Tài liệu tham khảo 36

Mở đầu

Trong chương trình đại học sinh viên chỉ được giới thiệu bài toán giá trị

riêng bậc một tiêu chuẩn. Trong khi đó, có rất nhiều bài toán, đặc biệt là trong

lĩnh vực cơ học, được qui về bài toán giá trị riêng bậc hai, ta có thể đưa về bài

toán giá trị riêng suy rộng bậc một, nhưng mặt khác có thể nghiên cứu độc

lập. Trong luận văn này chúng tôi sẽ tìm hiểu và trình bày "Bài toán giá trị

riêng bậc hai". Ngoài phần mở đầu và kết luận, luận văn này gồm hai chương.

Chương 1: Bài toán giá trị riêng bậc hai. Chương này chúng tôi trình bày

khái niệm về bài toán giá trị riêng bậc hai, tính chất và ứng dụng của bài toán

giá trị riêng bậc hai.

Chương 2: Giải số bài toán giá trị riêng bậc hai. Chương này trình bày một

vài phương pháp giải bài toán giá trị riêng bậc hai. Chúng tôi chia bài toán

ra hai loại dựa trên kích thước bài toán và dạng dữ liệu. Bài toán đặc (thông

thường) và cỡ của bài toán nhỏ. Còn bài toán thưa là bài toán có kích cỡ lớn

nhưng dữ liệu dạng thưa. Căn cứ vào đặc điểm của từng bài toán, phương

pháp giải cũng có nhiều khác biệt. Chúng tôi trình bày phương pháp Newton

và phương pháp phân tích Schur cho bài toán đặc và phương pháp dựa trên

không gian con Krylov cho bài toán thưa. Ngoài ra có thêm một vài ví dụ số

để minh họa cho các phương pháp trên.

Luận văn được hoàn thành tại Trường Đại học Khoa học - Đại học Thái

Nguyên. Em muốn gửi lời biết ơn sâu sắc nhất tới thầy giáo TS. Nguyễn

Thanh Sơn đã giúp đỡ, hướng dẫn tận tình và đầy trách nhiệm để em hoàn

thành luận văn này. Tôi cũng xin được gửi lời cảm ơn chân thành tới các thầy

giáo, cô giáo của Trường Đại học Khoa học - Đại học Thái Nguyên, gia đình

tôi và các bạn lớp cao học toán K9Y đã tạo điều kiện thuận lợi cho tôi trong

quá trình học cao học và thực hiện bản luận văn này. Trong quá trình viết luận

văn không tránh khỏi sai sót rất mong nhận được sự góp ý chân thành của độc

giả.

Thái Nguyên, ngày 10 tháng 11 năm 2017

Tác giả luận văn

Phí Thị Nho

Danh sách ký hiệu

Trong toàn luận văn, ta dùng những ký hiệu với các ý nghĩa xác định trong

bảng dưới đây:

ma trận A.

chuyển vị của ma trận thực A.

A AT A∗ = (A)T liên hợp phức của ma trận A.

A liên hợp của số phức của ma trận A.

ker(A) nhân của ma trận A.

span(A) không gian con sinh bởi các cột của ma trận A.

A > 0 (A ≥ 0) ma trận A xác định dương (nửa xác định dương).

deg(P) bậc của đa thức P.

det(A) định thức của ma trận A.

rank(B) hạng của ma trận B.

QEP bài toán giá trị riêng bậc hai.

||x|| chuẩn Ơclit.

gradien của P.

∇P K j(x, A) không gian con Krylov .

Chương 1

Bài toán giá trị riêng bậc hai

Nội dung chính của chương này là định nghĩa, tính chất và một số ứng

dụng của bài toán giá trị riêng bậc hai. Tuy nhiên, chúng tôi cũng dành một

thời lượng đáng kể cho việc trình bày bài toán giá trị riêng tiêu chuẩn và bài

toán giá trị riêng suy rộng. Lí do chính là ta có thể chuyển một bài toán giá

trị riêng bậc hai về bài toán giá trị riêng bậc một để giải. Thêm vào đó, nhiều

phương pháp giải số bài toán bậc hai cũng xuất phát từ những ý tưởng tương

tự cho bài toán bậc một. Khi viết chương này chúng tôi đã tham khảo các tài

liệu [1–3].

1.1 Bài toán giá trị riêng tiêu chuẩn

1.1.1 Khái niệm

Định nghĩa 1.1.1. Cho một ma trận vuông A ∈ Rn×n. Tìm đại lượng vô hướng λ ∈ R và véc tơ x ∈ Rn, x (cid:54)= 0, sao cho:

Ax = λ x (1.1)

hay

(A − λ I)x = 0 (1.2)

có một nghiệm không tầm thường. Cho cặp (λ , x) là một nghiệm của (1.1)

hoặc (1.2) tương ứng thì

(1) λ được gọi là một giá trị riêng của A;

(2) x được gọi là một véc tơ riêng của A;

(3) (λ , x) được gọi là cặp riêng của A;

(4) Đặt σ (A) của tất cả các giá trị riêng của A được gọi là phổ của A;

(5) Tất cả các véc tơ riêng có cùng giá trị riêng λ cùng với véc tơ 0 tạo thành một không gian con tuyến tính của R gọi là không gian riêng của λ ;

(6) Một nghiệm không tầm thường y của y∗A = λ y∗ được gọi là véc tơ riêng trái tương ứng với λ . Một véc tơ riêng trái của A là một véc tơ riêng phải của AT tương ứng với các giá trị riêng λ , ta viết : AT y = λ y.

1.1.2 Một số thuật toán tìm giá trị riêng

• Cơ sở trực giao cho không gian Krylov

Cho không gian con Krylov K j(x) = K j(x, A) ta có thể lấy

{x, Ax, ..., A( j−1)x},

làm hệ sinh. Tuy nhiên các vectơ Akx hội tụ đến vectơ riêng tương ứng với giá trị riêng có modul lớn nhất của A nên chúng sẽ sớm có xu hướng phụ thuộc

tuyến tính. Do đó quá trình trực giao hoá Gram-Schmidt được áp dụng cho

các vectơ trong cơ sở này để tìm một cơ sở trực chuẩn của không gian Krylov. Giả sử {q1, q2, ..., qi} là cơ sở trực chuẩn cho K i(x), ở đây i ≤ j. Chúng

ta xây dựng vectơ q j+1 bằng cách trực chuẩn hóa A jx với q1, q2, ..., q j

i A jx,

j ∑ i=1

y j := A jx − qiq∗

và sau đó chuẩn hóa các vectơ kết quả

q j+1 = y j/||y j||.

Ta được {q1, q2, ..., q j+1} là một cơ sở trực chuẩn của K j+1(x), gọi chung là cơ sở Arnoldi. Các vectơ được gọi là vectơ Arnoldi tương ứng. Các vectơ qi

có thể được tính như sau

K j+1(x, A) = ℜ([x, Ax, ..., A jx]), (q1 = x/||x||),

(Aq1 = αq1 + β q2, β (cid:54)= 0),

= ℜ([q1, Aq1, ..., A jq1]) = ℜ([q1, αq1 + β q2, A(αq1 + β q2), ..., A j−1(αq1 + β q2)]), = ℜ([q1, q2, Aq2, ..., A j−1q2]), ... = ℜ([q1, q2, ..., q j−1, Aq j]).

Vì vậy, thay vì trực chuẩn hoá A jq1 với q1, q2, ..., q j, chúng ta có thể trực chuẩn hoá Aq j với q1, q2, ..., q j để có được q j+1. Điều này có lợi về mặt tính toán vì sẽ giúp giảm số phép tính cần thực hiện. Các thành phần r j của Aq j trực chuẩn với q1, q2, ..., q j được cho bởi

i Aqi).

j ∑ i=1

(1.3) qi(q∗ r j = Aq j −

Nếu r j = 0 thì thủ tục dừng lại. Điều đó có nghĩa là chúng ta đã tìm thấy một không gian con bất biến, cụ thể là span{q1, q2, ..., q j}. Nếu ||r j|| > 0 ta có được q j+1

. q j+1 = r j ||r j||

Do q j+1 và r j cùng phương, ta có

j+1Aq j

j+1r j = ||r j|| = q∗ q∗

(1.4)

i Aq j thì (1.3) - (1.4) có thể được viết

Phương trình cuối cùng được suy ra từ việc q j+1 trực chuẩn với tất cả các vectơ Arnoldi trước đó. Đặt hij = q∗

j+1 ∑ i=1

(1.5) Aq j = qihij.

Tìm một cơ sở trực giao của không gian Krylov

Algorithm 1 Thuật toán Arnoldi Require: A ∈ R. Ensure: Một cơ sở trực giao cho K x x 1: Tính q1 = (cid:107)x2(cid:107) 2: for j = 1, . . . do r := Aq j for i = 1, . . . , j do

i r; r := r − qihi j;

hi j := q∗

end for h j+1, j := (cid:107)r(cid:107); if h j+1 j = 0 then

10:

11:

(cid:1) , H ∈ R( j+1)× j return (cid:0)q1, q2, ..., q j

12: end for 13: return (q1, q2, ..., qk+1) , H ∈ R(k+1)×k 14: Thuật toán Arnoldi dừng lại nếu h j+1 j = 0.

; end if q j+1 = r h j+1 j

Định nghĩa Qk = [q1, q2, ..., qk] thì phương trình (1.5) được viết như sau

(1.6) AQk = QkHk + [0, . . . 0, qk+1hk+1k] .

Phương trình (1.6) được gọi là quan hệ Arnoldi.

• Cơ sở Lanczos

Cơ sở Lanczos được xây dựng giống như cơ sở Arnoldi nhưng áp dụng cho ma trận Hermite hoặc đối xứng thực. Bằng cách nhân trái (1.6) với Q∗ k

chúng ta được

kAQk = Q∗

kQkHk = Hk.

Q∗

Nếu A là Hermite thì Hk cũng là Hermite và do đó là ba đường chéo. Chúng ta ký hiệu ma trận ba đường chéo này là Tk. Do tính đối xứng phương trình

(1.3) được đơn giản hoá như sau

jAq j) (cid:124) (cid:123)(cid:122) (cid:125) α j∈R

j−1Aq j) (cid:125) (cid:123)(cid:122) β j−1∈F

(1.7) r j = Aq j − qi (q∗ −q j−1 (q∗ = Aq j − α jq j − β j−1q j−1. (cid:124)

j+1(Aq j − α jq j − β j−1q j−1)

Tương tự như trên, chúng ta nhân từ bên trái (1.7) với q j+1 để có được

j+1r j = q∗ j+1Aq j = β j.

||r j|| = q∗ = q∗

Từ đây suy ra β j ∈ R, do đó

(1.8) β jq j+1 = r j, β j = ||r j||.

Từ (1.7) - (1.8) ta nhận được

Aq j = β j−1q j−1 + α jq j + β jq j+1.

Thu thập các phương trình với j = 1, ..., k chúng ta có

  α1 β1

β1 α2 β2

(1.9) AQk = Qk +βk[0, ..., 0, qk+1].

β2 α3 . . .                   . . . . . . βk−1 βk−1 αk

(cid:125) (cid:124) (cid:123)(cid:122) Tk

Như vậy, Tk ∈ Rk×k là đối xứng thực. Đẳng thức (1.9) được gọi là quan hệ Lanczos.

Thuật toán Lanczos được trình bày trong Thuật toán 2 này có ba vectơ q, r

và v được sử dụng.

Cho (λi, si) là một cặp riêng của Tm, tức là

(1.10) T si = λisi.

Khi đó, AQmsi = QmTmsi = λiQmsi. Vì vậy, giá trị riêng của Tm cũng là giá trị riêng của A. Các vectơ riêng của A tương ứng với giá trị riêng λi là

(1.11) yi = Qmsi.

Ta có thuật toán sau

1: Cho A ∈ Fn×n là Hermitian. Thuật toán này sẽ tìm mối quan hệ Lanczos (1.9), nghĩa là một cơ sở trực giao Qm = [q1, .., qm] cho K m(x) trong đó m là chỉ số nhỏ nhất sao cho K m(x) = K m+1(x), và (các yếu tố không tầm thường của) các ma trận ba đường chéo Tm .

2: q := x/||x||; Q1 = [q]; 3: r := Aq; 4: α1 := q∗r; 5: r := r − α1q; 6: β1 := ||r||; 7: for j = 2, 3, ... do

10:

11:

12:

13:

Algorithm 2 Thuật toán Lanczos cơ bản để tìm một cơ sở trực giao cho không gian Krylov K m(x).

14:

v = q; q := r/β j−1; Q j := [Q j−1, q]; r := Aq − β j−1v; α j := q∗r; r := r − α jq; β j := ||r||; if β j = 0 then

15:

trả lại (Q ∈ Fn× j; α1, ..., α j; β1, ..., β j−1);

16: end for

end if

Sau khi tính được các cơ sở trực giao V của một không gian con Krylov

của ma trận A, thay vì tính giá trị riêng của A, ta tính giá trị riêng của ma trận

V T AV có cỡ bé hơn ma trận A ban đầu. Người ta đã chứng minh được rằng các giá trị riêng của V T AV xấp xỉ rất tốt các giá trị riêng của A theo thứ tự

giảm dần của mô đun. Cụ thể, xin xem trong tài liệu [1].

1.1.3 Bài toán giá trị riêng suy rộng

• Chùm ma trận chính qui và dạng chính tắc Weierstrass

Định nghĩa 1.1.2. Ma trận A − λ B, ở đó A, B là những ma trận cỡ m × n, được gọi là chùm ma trận hoặc tổng quát hơn là một λ ma trận bậc một, trong đó

λ là tham số.

Định nghĩa 1.1.3. Nếu A, B là ma trận vuông cỡ n × n và det(A − λ B) không đồng nhất bằng không, chùm A − λ B được gọi là chính qui. Nếu không, nó

được gọi là suy biến. Khi A − λ B là chính qui, P(λ ) = det(A − λ B) được gọi

là đa thức đặc trưng của A − λ B. Các giá trị riêng của A − λ B được xác định

bởi

(1) Nghiệm của P(λ ) = 0

(2) ∞ (với bội n − deg(P)) nếu deg(P) < n

Ví dụ 1.1.4. Cho

    1 0 0 2 0 0

A − λ B = − λ 0 1 0 0 0 0             0 0 0 0 0 1

Ta có

P(λ ) = det(A − λ B)

= (1 − 2λ )(1 − 0λ )(0 − λ )

= (2λ −1)λ .

, 0, và ∞. Vậy các giá trị riêng là: 1 2

Mệnh đề 1.1.5. Cho A − λ B là chùm ma trận chính qui cỡ n × n. Nếu B không suy biến, mọi giá trị riêng của A − λ B là hữu hạn và trùng với các giá trị riêng của AB−1 hoặc B−1A. Nếu B là suy biến, A − λ B có giá trị riêng ∞

bội n − rank(B). Nếu A không suy biến, giá trị riêng của A − λ B trùng với nghịch đảo của các giá trị riêng của A−1B hoặc BA−1, trong đó một giá trị riêng không của A−1B tương ứng với một giá trị riêng ∞ của A − λ B.

Chứng minh. Nếu B không suy biến thì B−1tồn tại và det(B), det(B−1) đều

khác không. Ta có

det(AB−1 − λ I) = det((A − λ B)B−1)

= det(A − λ B)det(B−1).

Theo đó

det(AB−1 − λ I) = 0 ⇔ det(A − λ B) = 0,

tức là tập giá trị riêng của A − λ B và AB−1 là trùng nhau. Khẳng định cho B−1A được chứng minh tương tự.

Trong trường hợp B suy biến, giả sử rank(B) = k < n. Kí hiệu phân tích

giá trị riêng kì dị của B là

B = UΣV T ,

trong đó Σ là ma trận chéo chỉ có k phần tử đường chéo chính khác không và

U,V là các ma trận trực giao. Khi đó

det(A − λ B) = det(A − λUΣV T )

= det(U(U T AV − λ Σ)V T ) = det(U)det(U T AV − λ Σ)det(V T ) = ±det(U T AV − λ Σ),

do định thức của ma trận trực giao chỉ có thể bằng ±1. Khi đó rõ ràng det(U T AV − λ Σ) là một đa thức bậc k nên ∞ là giá trị riêng bội n − k.

Cuối cùng, nếu A không suy biến thì A−1 tồn tại. Ta có ngay

det(A − λ B) = det(A(I − λ A−1B)) = det(A)det(I − λ A−1B).

, Nếu giả sử λ (cid:54)= 0, đặt µ = 1 λ

A−1B) det(A − λ B) = det(A)det(I − 1 µ

= −µdet(A)det(A−1B − µI).

là giá trị riêng của

1 Hay λ là giá trị riêng của A − λ B khi và chỉ khi µ = λ A−1B. Khẳng định cho BA−1 được chứng minh tương tự.

Định nghĩa 1.1.6. Cho λ (cid:48) là một giá trị riêng hữu hạn của chùm chính qui A − λ B. x (cid:54)= 0 là một véc tơ riêng phải nếu (A − λ (cid:48)B)x = 0 hoặc tương đương Ax = λ (cid:48)Bx. Nếu λ (cid:48) = ∞ là một giá trị riêng và Bx = 0 thì x là véc tơ riêng phải. Một véc tơ riêng trái của A − λ B là một véc tơ riêng phải của (AT − λ BT ).

Ví dụ 1.1.7. Xét dao động tắt dần của chất điểm trong hệ lò xo. Có hai chùm

ma trận phát sinh tự nhiên từ bài toán này. Thứ nhất, có bài toán giá trị riêng

(cid:35) (cid:34) −M−1B −M−1K Ax = x = λ x. I 0

Khi

(cid:35) (cid:34) (cid:35) (cid:34) −B −K M 0 x. x = λ I 0 0 I

Đây là một công thức tốt hơn nếu M là điều kiện rất xấu, do đó M−1B và

M−1K là khó tính chính xác .

Thứ hai, nó rất phổ biến để xem xét các trường hợp B = 0 (không tắt dần)

do đó phương trình vi phân ban đầu là

M ¨x(t) + Kx(t) = 0.

Các cách giải của dạng

xi(t) = eλitxi(0),

i eλitMxi(0) + eλitKxi(0) = 0, λ 2

ta có

hoặc

λ 2 i Mxi(0) + Kxi(0) = 0.

Nói cách khác, −λ 2 i là một giá trị riêng và xi(0) là véc tơ riêng phải của chùm K − λ M. Khi ta giả thiết M không suy biến, đây cũng là một giá trị riêng và véc tơ riêng phải của M−1K. Định nghĩa tiếp theo cho thấy làm thế nào để

tổng quát khái niệm về sự tương đồng với ma trận chùm.

Định nghĩa 1.1.8. Cho PL và PR là ma trận không suy biến. Khi đó, chùm A − λ B và PLAPR − λ PLBPR được gọi là tương đương.

Nhận xét 1.1.9. Đây là một mở rộng khái niệm ma trận tương đương.

Mệnh đề 1.1.10. Chùm ma trận chính qui tương đương A − λ B và PLAPR − λ PLBPR có cùng giá trị riêng, x là một véc tơ riêng phải của A − λ B nếu và chỉ nếu P−1 R x là một véc tơ riêng phải của PLAPR − λ PLBPR, y là véc tơ riêng trái của A − λ B nếu và chỉ nếu (PT L )−1y là một véc tơ riêng trái của

PLAPR − λ PLBPR.

Chứng minh. Mệnh đề được chứng minh nhờ ba lập luận sau đây

R x = 0.

det(A − λ B) = 0 nếu và chỉ nếu det(PL(A − λ B))PR = 0.

L y)T PL(A − λ B)PR = 0.

(A − λ B)x = 0 nếu và chỉ nếu PL(A − λ B)PRP−1 yT (A − λ B) = 0 nếu và chỉ nếu (P−T

Định lí sau sẽ khái quát công thức về dạng chính tắc Jordan cho ma trận

đơn. Đối với chùm ma trận chính qui, nó được gọi là dạng chính tắc Weier-

strass.

Định lí 1.1.11. Cho A − λ B là chính qui, khi đó tồn tại Pm ma trận khả nghịch sao cho

PL(A − λ B)PR = diag(Jn1(λ1) − λ In1, ..., Jnk(λnk) − λ Ink, Nm1, ..., Nmr).

Ở đó, Jniλi là một khối Jordan cỡ ni × ni với giá trị riêng λi,

  λi

1 . . . , Jni(λi) = . . . . . .             . . . 1 λi

và Nmi là một khối Jordan cho λ = ∞ với bội mi

 

1 λ . . . Nmi = = Imi − λ Jmi(0). . . . . . .             . . . λ 1

1.2 Bài toán giá trị riêng bậc hai

1.2.1 Khái niệm

Định nghĩa 1.2.1. Cho các ma trận thực vuông M,C, K cỡ n × n. Kí hiệu

Q(λ ) = Mλ 2 +Cλ + K.

và gọi đó là một λ - ma trận bậc hai (hay vắn tắt là λ - ma trận). Bài toán giá trị riêng bậc hai là bài toán tìm số λ ∈ C sao cho có véc tơ khác không x ∈ Cn

để cho

(Mλ 2 +Cλ + K)x = 0. (1.12)

Ngoài ra, ta cũng xét bài toán tìm y, λ sao cho

y∗(Mλ 2 +Cλ + K) = 0, (1.13)

trong đó ∗ là kí hiệu liên hợp phức. Các véc tơ x, y lần lượt là các véc tơ riêng phải và véc tơ riêng trái tương ứng với các giá trị riêng λ . Bài toán giá trị

riêng bậc hai được gọi là đối xứng nếu M,C, K là các ma trận đối xứng. Đặt m(x) = x∗Mx, c(x) = x∗Cx, k(x) = x∗Kx với x ∈ Cn khác không.

Nếu bài toán giá trị riêng bậc hai đối xứng thì đây là hàm có giá trị thực.

Cũng giống như bài toán giá trị riêng tiêu chuẩn, bội đại số của một giá trị riêng λ0 là bậc α của nghiệm λ0 của đa thức det(Q(λ )). Bội hình học của γ của λ0 là số chiều của hạt nhân ker(Q(λ0)). Giá trị riêng được gọi là đơn nếu α = γ = 1 và được gọi là nửa đơn nếu α = γ. Giá trị riêng được gọi là

khuyết nếu nó không phải nửa đơn.

Véc tơ x1 được gọi là một véc tơ riêng suy rộng tương ứng với giá trị riêng

λ0 nếu x1 là nghiệm của phương trình

(1.14) Q(λ0)x1 = −Q(cid:48)(λ0)x0,

với x0 là một véc tơ riêng ứng với λ0, Q(cid:48)(λ0) là đạo hàm của Q theo từng phần tại λ0.

Ta có thể chỉ ra một giá trị riêng nửa đơn thì sẽ không có véc tơ riêng suy

rộng đi kèm. Ta cũng biết rằng nếu giá trị riêng khuyết thì số véc tơ riêng

tương ứng sẽ không lấp đủ bội đại số của nó. Khái niệm véc tơ riêng suy rộng

được định nghĩa như sau.

Định nghĩa 1.2.2. Các véc tơ x0, x1, ..., xm−1 được gọi là một xích Jordan có độ dài m của Q(λ ) tương ứng giá trị riêng λ0 nếu m đẳng thức sau được thỏa

mãn

Q(λ0)x0 =0, Q(λ0)x1 + Q(cid:48)(λ0)x0 =0,

Q(cid:48)(cid:48)(λ0)x0 =0, Q(λ0)x2 + Q(cid:48)(λ0)x1 + 1 2

· · · · · ·

Q(λ0)xm−1 + Q(cid:48)(λ0)xm−2 + Q(cid:48)(cid:48)(λ0)xm−3 =0, 1 2

với x0 (cid:54)= 0 là một véc tơ riêng tương ứng với λ0. Khi đó, các véc tơ

x1, x1, ..., xm−1,

được gọi là các véc tơ riêng suy rộng.

• Có thể coi véc tơ riêng x0 là một xích Jordan có độ dài Nhận xét 1.2.3.

1. Vì thế khái niệm này mở rộng khái niệm véc tơ riêng thông thường.

Chú ý 1.2.4. Các véc tơ riêng suy rộng không nhất thiết độc lập tuyến tính. Nói riêng, có thể có vài véc tơ riêng suy rộng ứng λ0 bằng 0.

• Bài toán đối xứng (1.12) QEP được gọi là hyperbolic Định nghĩa 1.2.5.

nếu

c(x)2 > 4m(x)k(x), ∀x (cid:54)= 0, x ∈ Cn

• Bài toán QEP được gọi là tựa hyperbolic nếu

c(x)2 > 4m(x)k(x),

với mọi x là véc tơ riêng của Q(λ ).

• Bài toán QEP được gọi là elliptic nếu

c(x)2 < 4m(x)k(x),

với mọi giá trị riêng của Q(λ ).

1.2.2 Tuyến tính hóa bài toán giá trị riêng bậc hai

Định nghĩa 1.2.6. Cho λ - ma trận bậc hai cỡ n × n

Q(λ ) = Mλ 2 +Cλ + K. (1.15)

Một λ - ma trận bậc một A − λ B được gọi là một tuyến tính hóa của (1.15)

nếu tồn tại 2 λ - ma trận E(λ ), F(λ ) cỡ 2n × 2n có định thức là hằng số khác

không sao cho

(cid:34) (cid:35) Q(λ ) 0 =E(λ )(A − λ B)F(λ ). 0 In

Bổ đề 1.2.7. Tập các giá trị riêng của Q(λ ) và A − λ B là trùng nhau.

Chứng minh. Thật vậy, λ0 là một giá trị riêng của Q(λ ) khi và chỉ khi det(Q(λ0)) = 0. Ta có

det(Q(λ0)) = det(Q(λ0))det(In)

(cid:32) (cid:34) (cid:35) (cid:33)

= det 0 Q(λ 0) 0 In

= det(E(λ0)(A − λ0B)F(λ0))

= det(E(λ0))det(A − λ0B)det(F(λ0)),

do giả thiết det(E(λ0)) (cid:54)= 0, det(F(λ0)) (cid:54)= 0, det(Q(λ0)) = 0 khi và chỉ khi det(A − λ0B) = 0.

Tuyến tính hóa luôn tồn tại và không duy nhất. Người ta cố gắng chọn

những tuyến tính hóa sao cho nó bảo toàn các tính chất của λ -ma trận ban

đầu hoặc thuận lợi cho quá trình tính toán. Sau đây là hai dạng tuyến tính hóa

phổ biến.

(cid:34) (cid:35) (cid:34) (cid:35) N 0 N 0 − λ L1 : −K −C 0 M

và

(cid:34) (cid:35) (cid:34) (cid:35) C M −K 0 , − λ L2 : 0 N N 0

trong đó N là một ma trận không suy biến cỡ n × n. Ta có thể chỉ ra, chẳng hạn dạng L1, ma trận E(λ ) và F(λ ) như sau

(cid:35) (cid:34) (C + λ M)N−1 −I , E(λ ) = N−1 0

(cid:34) I (cid:35) 0 , F(λ ) = I λ I

Để thấy rõ hơn mối liên hệ giữa giá trị riêng, véc tơ riêng của bài toán

QEP và tuyến tính hóa của nó. Đặt u = λ x, khi đó bài toán QEP được đưa về

dạng

λ Mu + λCx + Kx = 0,

hay

(cid:34) (cid:35) (cid:34) (cid:35) (cid:34) (cid:35) (cid:34) (cid:35) I x x 0 I 0 = 0. − λ u u −K −C 0 M

Phát biểu sau đây sẽ khẳng định việc chọn bất kì N

Bổ đề 1.2.8. Bài toán giá trị riêng L1, L2 là bất biến đối với N

(cid:34) (cid:35) N 0 Chứng minh. Xét bài toán L1, do N khả nghịch nên ma trận 0 M

cũng khả nghịch và

(cid:34) (cid:35)−1 (cid:35) (cid:34) N−1 N 0 0 = . 0 M 0 M−1

Khi đó, bài toán giá trị riêng L1 tương đương với bài toán

(cid:34) (cid:35)−1 (cid:34) (cid:35) N N 0 0 − λ I, −K −C 0 M

hay

(cid:34) (cid:35) I 0 − λ I. (1.16) −M−1K −M−1C

Rõ ràng bài toán giá trị riêng (1.16) không phụ thuộc vào việc chon N. Việc chứng minh cho L2 được thực hiện tương tự.

1.2.3 Bộ ba Jordan của Q(λ )

Nhắc lại rằng một phân tích Jordan của ma trận A là

A = XJX −1. (1.17)

với J là ma trận chéo khối J = (J1, · · ·, Jk) và Ji, i = 1, k, là các khối Jordan

  λi

1 . . . ∈ Rmi×mi. Ji =

            . . . . . . 1 λi

Ma trận không suy biến X chứa các véc tơ riêng và các véc tơ riêng suy rộng

của A.

Mở rộng khái niệm này nâng lên cho một λ - ma trận bậc hai, ta thu được khái niệm bộ ba Jordan. Giả sử det(M) (cid:54)= 0. Giả sử J = diag(J1, · · ·, Jk) với Ji, i = 1, · · ·, k là các khối Jordan cỡ mi và m1 + · · · + mk = 2n. Ta gọi

mi−1],

0, · · ·, xi

Xi = [xi

là ma trận cỡ n × mi có các cột lập thành một xích Jordan độ dài mk ứng với giá trị riêng λi. Tiếp theo, ta đặt

X = [X1, · · ·, Xk] ∈ Rn×2n.

Khi đó, cặp (X, J) được gọi là một cặp Jordan của Q(λ ) và thỏa mãn ma trận

(cid:34) (cid:35) X , XJ

là không suy biến và

MXJ2 +CXJ + KX = 0. (1.18)

Xích Jordan trái được xây dựng tương tự cho ma trận Y cỡ 2n × n định nghĩa

như sau

(cid:34) X (cid:35)−1 (cid:34) (cid:35) 0 Y = M−1. XJ I

Nếu mối giá trị riêng đều là nửa đơn, chuyển vị của những hàng của Y lập

nên tập các véc tơ riêng trái của Q(λ ). Cuối cùng bộ ba (X, J,Y ) được gọi là

một bộ ba Jordan. Người ta đã chỉ ra

(1.19) XY M = 0 và XJY M = I.

1.2.4 Một số tính chất của bài toán giá trị riêng bậc hai

Mệnh đề 1.2.9. Giả sử các ma trận M,C, K là thực. Khi đó mọi nghiệm của bài toán giá trị riêng bậc hai đều là số thực hoặc là các cặp số phức liên hợp

nhau. Hơn nữa nếu x là véc tơ riêng phải tương ứng với λ thì x cũng là véc tơ

riêng phải ứng với λ .

Chứng minh. Theo định nghĩa, tồn tại véc tơ khác không x ∈ C, số λ ∈ C sao

cho

(Mλ 2 +Cλ + K)x = 0. (1.20)

Do đó định thức

P(λ ) := det(Mλ 2 +Cλ + K) = 0.

Do M,C, K thực nên P(λ ) là một đa thức của λ có hệ số thực. Khẳng định

thứ nhất được suy ra từ tính chất của nghiệm của đa thức hệ số thực.

Để chứng minh khẳng định thứ hai, ta lấy liên hợp hai vế của (1.20) với

lưu ý M,C, K là các ma trận thực, ta thu được

(Mλ +Cλ + K)x = 0.

Đây cũng là điều cần chứng minh.

Mệnh đề 1.2.10. Các bài toán giá trị riêng bậc hai đối xứng và đồng thời là hyperbolic hoặc tựa hyperbolic đều có giá trị riêng thực.

Chứng minh. Xét phương trình bậc hai với ẩn là λ Giả sử x là một véc tơ

riêng ứng với λ

x∗Mxλ 2 + x∗Cxλ + x∗Kx = 0. (1.21)

. λ = (1.22) Nghiệm của phương trình này có thể biểu thị dưới dạng −c(x) ± (cid:112)c(x)2 − 4m(x)k(x) 2m(x)

Khi đó, kết luận của mệnh đề được suy ra từ tính không âm của đại lượng c(x)2 − 4m(x)k(x).

Mệnh đề 1.2.11. Bài toán QEP đối xứng, elliptic chỉ có giá trị riêng phức.

Chứng minh. Kết luận này được suy ra trực tiếp từ định nghĩa và dạng biểu

diễn (1.22) của giá trị riêng.

Một trường hợp riêng thú vị của mệnh đề trên được phát biểu như sau

Hệ quả 1.2.12. Giả sử M, K là các ma trận đối xứng xác định dương, C là phản đối xứng, tức là CT = −C. Khi đó, bài toán QEP chỉ có giá trị riêng

thuần ảo.

Chứng minh. Ta có, với x ∈ C thì

c(x) = x∗Cx = x∗CT x = −x∗Cx = −c(x).

Điều đó suy ra c(x) là số thuần ảo. Hơn nữa, bình phương của một số thuần

ảo là một số âm. Do đó, kết hợp với giả thiết m(x), k(x) dương và dạng biểu

diễn (1.22) của giá trị riêng, ta suy ra điều cần chứng minh.

Định nghĩa 1.2.13. Cho ma trận M,C đối xứng xác định dương, K đối xứng nửa xác định dương. Điều kiện tắt dần (over damping) được định nghĩa là

((x∗Cx)2 − 4(x∗Mx)(x∗Kx)) > 0 . (1.23) min ||x||2=1

Một hệ cơ học được miêu tả dưới dạng hệ phương trình vi phân có các ma

trận hệ số M,C, K thỏa mãn điều kiện (1.23) được gọi là hệ tắt dần.

Mệnh đề 1.2.14. Hệ tắt dần chỉ có các giá trị riêng thực, âm.

Chứng minh. Dễ thấy, một hệ tắt dần là hyperbolic. Áp dụng mệnh đề 1.2.10

ta suy ra rằng các giá trị riêng của nó là thực. Dấu của nó là âm được suy ra

từ việc áp dụng định lí Viete cho phương trình bậc hai (1.2.13).

1.3 Một số ứng dụng khác của bài toán giá trị riêng bậc

hai

1.3.1 Biểu diễn nghiệm của phương trình vi phân tuyến tính cấp

hai

Xét phương trình vi phân tuyến tính cấp hai

M ¨q(t) +C ˙q(t) + Kq(t) = f (t). (1.24)

Phương trình này là mô hình toán học cho nhiều hiện tượng và đặc biệt được

sử dụng nhiều trong cơ học cấu trúc. Trong (1.24), M,C, K là các ma trận thực

cỡ n × n.

Giả sử M không suy biến và (X, J,Y ) là bộ ba Jordan của Q(λ ). Từ đẳng

thức (1.18), ta có thể kiểm tra được rằng

q(t) = XeJta, (1.25)

với a ∈ R2n là một véc tơ hằng bất kì là nghiệm tổng quát của phương trình

thuần nhất tương ứng với (1.24). Thật vậy

(1.26) ˙q(t) = XJeJta, ¨q(t) = XJ2eJta.

Thay (1.25), (1.26) vào (1.24) với f (t) ≡ 0,

M ¨q(t) +C ˙q(t) + Kq(t) = (MXJ2 +CXJ + KX)eJta = 0

(cid:90) t

Tiếp theo, ta chỉ ra

e−JsY f (s)ds qp(t) = XeJt

(cid:90) t

là một nghiệm riêng của (1.24). Thật vậy

0 (cid:90) t

e−JsY f (s)ds ˙qp(t) = XJeJt

e−JsY f (s)ds + XJY f (t). ¨qp(t) = XJ2eJt

(cid:90) t

Vì thế, kết hợp (1.25) và (1.26), ta thu được

e−JsY f (s)ds + MXJY f (t) M ¨qp(t) +C ˙q(t) + Kqp(t) = (MXJ2 +CXJ + KX)eJt

= f (t).

(cid:90) t

Do vậy, nghiệm tổng quát của (1.24) là

q(t) = XeJt(a + e−JsY f (s)ds).

1.3.2 Bài toán hạn chế bình phương nhỏ nhất

Cho A ∈ Rn×n là ma trận đối xứng, b ∈ Rn. Xét bài toán bình phương tối

thiểu có ràng buộc

min{xT Ax − 2bT x : xT x = α 2}. (1.27)

Bài toán này có thể được đưa về một bài toán giá trị riêng bậc hai bằng cách

áp dụng phương pháp nhân tử Lagrange. Cho

φ (x, λ ) = xT Ax − 2bT x − λ (xT x − α 2).

Đạo hàm ϕ theo x, và λ sinh ra các phương trình

Ax − λ x = b, α 2 = xT x. (1.28)

Ta cần tìm nghiệm λ nhỏ nhất của các phương trình này để giải bài toán

(1.27). Giả sử λ không là một giá trị riêng của A. Ta đặt

y = (A − λ I)−2b = (A − λ I)−1x.

Khi đó (1.28) tương đương với:

bT y − α 2 = 0, (1.29)

(A − λ I)2y = b. (1.30)

Từ ( 1.29) ta có :

bT y α 2 = 1. Bằng cách khai triển (1.30) ta thu được bài toán giá trị riêng bậc hai đối xứng

)y = 0. (λ 2I − 2λ A + A2 − α −2bbT (cid:17) (cid:16) (1.31)

Nghiệm của (1.27) là

x = (A − λ I)−1 b, (1.32)

trong đó λ là giá trị riêng nhỏ nhất của (1.32).

1.3.3 Một vài ví dụ

Ví dụ 1.3.1. Cho λ ma trận bậc hai

Q1(λ ) = Mλ 2 + Dλ + K,

với

      0.5 1.75 1 1 0.2

, K = . M = , D = 1.5 7.5 2 1 1                   2.5 1 5.0 0.2 1

Dễ thấy cả ba ma trận hệ số M, D, K là các ma trận đối xứng, xác định dương.

Ngoài ra, ta kiểm tra được nó thỏa mãn điều kiện tắt dần. Do vậy, hệ cơ học được miêu tả bởi Q1λ là một hệ tắt dần.

Sử dụng phần mềm MATLAB, câu lệnh

[X, λ ] = polyeig(K, D, M),

ta tính được các giá trị riêng, lấy xấp xỉ 4 chữ số thập phân sau dấu phẩy là

λ1 = −4.7586, λ2 = −2.6614, λ3 = −1.6266,

λ4 = −1.1556, λ5 = −0.2713, λ6 = −0.0264.

và ma trận các véc tơ riêng   0.2415 −0.9901 −0.6870 0.9141 0.6933 −0.5228

−0.9700 −0.1268 −0.2213 0.2623 −0.2549 0.6165       0.0273

0.0601 −0.6921 0.3092 −0.6740 −0.5887 Như vậy, điều này đã minh họa cho mệnh đề 1.2.14 về giá trị riêng của

một hệ tắt dần. Ngoài ra, có thể kiểm tra thấy 3 véc tơ riêng đầu tiên độc lập

tuyến tính. Các véc tơ riêng khác đều biểu thị tuyến tính qua 3 véc tơ này.

Ta sẽ xét một tuyến tính hóa bảo toàn tính đối xứng của Q(λ ) nói trên

được định nghĩa như sau

(cid:34) (cid:35) (cid:34) (cid:35) K 0 K . , E1 = A1 = K D −M

và một tuyến tính hóa khác không bảo toàn tính đối xứng

(cid:34) (cid:35) (cid:34) (cid:35) K K 0 . A2 = , E2 = M −K −D

Bằng cách sử dụng các lệnh eig(A1, E1), và eig(A2, E2), ta thu được tập các giá trị riêng λ1, λ2, λ3, λ4, λ5, λ6, như đề cập ở trên.

Ví dụ 1.3.2. Xét λ -ma trận bậc hai

Q2(λ ) = Mλ 2 + Dλ + K,

  −1.5 1 0

với M, D giống như ở ví dụ 1 và K cho bởi 1 1 3     .   0 1 −1

Dễ thấy cả 3 ma trận đều đối xứng và ta có thể kiểm tra Q2(λ ) thỏa mãn điều kiện hyperbolic. Tuy nhiên, do K không nửa xác định dương nên Q2(λ )

là một hệ hyperbolic nhưng không phải hệ tắt dần. Bằng cách tương tự như ví dụ 1.3.1, ta tính được các giá trị riêng của Q2(λ ) và các tuyến tuyến tính hóa

của chúng như sau:

λ1 = −4.6886, λ2 = −4.0610, λ3 = 0.7556,

λ4 = −0.5652, λ5 = 0.2132, λ6 = −2.1540.

Điều này minh họa cho mệnh đề 1.2.10

Chương 2

Giải số bài toán giá trị riêng bậc hai

Trong chương này, chúng tôi trình bày một số phương pháp số giải bài

toán giá trị riêng bậc hai. Chúng tôi chia bài toán ra hai loại dựa theo dạng dữ

liệu. Nếu bài toán có cỡ nhỏ đến vừa và dạng đặc, tức là số phần tử 0 trong

các ma trận không đáng kể, chúng tôi trình bày phương pháp dựa theo phân

tích Schur của ma trận và phương pháp Newton giải phương trình phi tuyến.

Để viết chương này, ngoài tài liệu chính [5], chúng tôi còn tham khảo các tài

liệu [2–4]

2.1 Phương pháp số cho bài toán đặc

2.1.1 Phương pháp Newton

Để áp dụng phương pháp Newton giải phương trình phi tuyến cho bài toán

QEP, ta đưa nó về một bài toán trong không gian n + 1 chiều như sau

(cid:34) (cid:35) (cid:34) (cid:35) x P := = 0. (2.1) Q(λ )x vT x − 1 λ

Đạo hàm Fre’chet của P được tính bằng

(cid:35)

. ∇P = (2.2) (cid:34) Q(λ ) Q(cid:48)(λ )x vT 0

Nhắc lại rằng phần chính của phương pháp Newton cho phương trình phi

tuyến (cid:34) (cid:35) x P = 0, λ

là dãy lặp (cid:34) (cid:35) (cid:34) (cid:35) (cid:34) (cid:35) xs xs = . − (∇P)−1P (2.3) λs

xs+1 λs+1 Sử dụng (2.1) và (2.2) và kí hiệu Qs = Q(λs) và Q(cid:48) λs s = Q(cid:48)(λs), (2.3) được đưa

về dạng

(cid:35) (cid:34) (cid:35) (cid:34) (cid:35)

= − . (2.4) xs+1 − xs λs+1 − λs (cid:34) Qs Q(cid:48) sxs vT 0 s Qsxs vT s xs − 1

s xs = 1. Khi đó (2.4) trở thành

Giả sử xs được chuẩn tắc hóa đối với vs, tức là vT

sxs,

Qsxs+1 = −(λs+1 − λs)Q(cid:48) (2.5)

vT s xs+1 = 1.

, ta viết lại (2.5) dưới dạng thuật toán như sau. Đặt us+1 = xs+1 (−λs+1 + λs)

1: k = 0;

2: while (k < N; ||Xk+1 − Xk|| > ε) do Giải Qsuk+1 = QsXk tìm uk+1;

Algorithm 3 Thuật toán Newton cho QEP. Require: M,C, K, x0, λ0, ε, N. Ensure: (X, λ ).

; Tính λk+1 = λk − V T k Xk (V T k uk+1)

7: end while 8: X = Xk; λ = λk;

Tính Xk+1 = Cuk+1; k = k + 1;

Trong thuật toán trên C là một ma trận chuẩn tắc hóa bất kì

Để kết thúc mục này, chúng tôi trình bày một số cách chọn dãy vs. Cách đơn giản nhất là chọn vs là một trong các véc tơ đơn vị ei. Nó có nghĩa là giữ thành phần i của véc tơ xs luôn là hằng số. Trong trường hợp ta cần tính một và i cặp riêng, ta nên chọn vs sao cho nó trực giao với các véc tơ đã tính trước

đó; điều này giúp thuật toán không hội tụ tới những véc tơ đã tính được.

2.1.2 Phân tích Schur thực suy rộng

Phân tích Schur được biết đến là một cách rất hữu hiệu để tìm các giá trị

riêng của một ma trận. Trường hợp đơn giản nhất là trên phân tích Schur phức

của một ma trận đơn A. Khi đó, các giá trị riêng của ma trận đơn A được hiển

thị trên đường chéo chính của dạng Schur T của A. Với T là ma trận phức, tam giác trên thỏa mãn T = QHAQ, và Q là ma trận phức, trong đó QH là liên hợp phức của Q.

Trong trường hợp A- ma trận thực bậc một, ta có kết quả sau đây

Mệnh đề 2.1.1. Giả sử A, B ∈ Rn×n. Khi đó, luôn tồn tại các ma trận trực giao Q, và Z sao cho QT AZ = φ là tựa tam giác trên và QT BZ = ψ là tam

giác trên.

Trong trường hợp det(B) (cid:54)= 0, các phần tử trên đường chéo chính của ψ

khác không. Khi đó, các giá trị riêng của λ - ma trận A − λ B được tính như

sau. Ta xem xét đường chéo của ma trận tựa chéo trên của φ là dạng khối,

φi, i = 1, ..., k, trong đó φi là một phần tử (hay ma trận cỡ 1 × 1) hoặc φi là ma trận cỡ 2 × 2 nếu có phần tử phía dưới đường chéo chính khác không. Ta

gọi ψi là ma trận khối con của ψ có cùng vị trí với φi trong φ . Theo đó, các

nếu ma trận đó là cỡ 1 × 1 và sẽ là hai giá trị giá trị riêng của A − λ B sẽ là φi ψi riêng của λ -ma trận cỡ 2 × 2 φi − λ ψi. Kĩ thuật này có tên là phân tích QZ.

Để áp dụng kĩ thuật này cho bài toán QEP, ta chỉ cần áp dụng nó cho một

dạng tuyến tính hóa A − λ B của nó.

2.2 Phương pháp số cho bài toán thưa

Khi cỡ của bài toán QEP lớn thì các phương pháp trình bày ở mục 2.1 sẽ

trở nên đắt đỏ. Trong mục này, chúng ta trình bày một phương pháp dựa trên

không gian con Krylov, phù hợp cho bài toán cỡ lớn. Mục này được viết dựa

trên việc tham kháo tài liệu [6].

Trước tiên, để bắt đầu chúng tôi nhắc lại sơ lược phương pháp Arnoldi

cho bài toán giá trị riêng suy rộng (λ A + B)x = 0. Bởi phần chính của phương pháp này là xây dựng một cơ sở trực chuẩn {q1, ..., qk} cho không gian con

Krylov

Kk(A−1B, q1) = span{q1, A−1Bq, ..., (A−1B)k−1qk}.

Kí hiệu Qk = [q1, ..., qk], khi đó hình chiếu của A−1B lên Kk(A−1B, q1) là ma trận Hk = QT k (A−1B)Qk = [hi j]. Đây là một ma trận Hessenberg, tức là hi j = 0 ∀i, j sao cho i − j ≥ 2. Khi đó, ta có mối quan hệ

(A−1B)Qk = QkHk + hk+1,kqk+1eT k ,

trong đó ek ∈ Rk là véc tơ đơn vị thứ k. Nếu (θ0, u) với ||u|| = 1 là một cặp riêng của −Hk thì (θ0, x) với x = Qku có thể được coi là một xấp xỉ của một cặp riêng của ma trận A−1B. Cặp (θ0, x) được gọi là cặp Ritz bao gồm giá trị Ritz và véc tơ Ritz. Thặng dư của việc xấp xỉ này là

(θ0A + B)x = hk+1,kukAqk+1.

Do hk+1,k sẽ giảm khi k tăng và hơn nữa uk cùng giảm nên thặng dư này sẽ dần về 0. Điều này dẫn đến sự hội tụ của phương pháp.

Ta sẽ áp dụng phương pháp này với bài toán QEP. Do bài toán liên quan

đến số hạng chứa ma trận K, nên ta sẽ áp dụng thuật toán Arnoldi để xây dựng ma trận Hessenberg Hk của M−1C. Sau đó ta tiếp tục tính hình chiếu của M−1K lên không gian Krylov này:

k M−1KQk.

Gk = QT

Khi đó, ta xét bài toán QEP chiếu

Lk(λ ) = Iλ 2 + Hkλ + Gk.

Nếu θ0, u là một cặp riêng của Lk(λ ), tức là

(Iθ 2 + Hkθ + Gk)u = 0,

thì ta có thể sử dụng (θ0, x = Qku) là một cặp riêng xấp xỉ của L(λ ). Ta cũng

gọi chúng là một cặp Ritz. Xét thặng dư tương ứng

2 +Cθ0 + K)x

rk = (Mθ0

k u + ∆ku

= A(Qkθ 2 + QkHkθ0 + QkGk)u + + θ0hk+1, kAq+1eT

= θ0hk+1, kukAqk+1 + ∆ku,

trong đó uk là thành phần thứ k của u và

∆k = KQk − MQkGk.

Trong khi số hạng thứ nhất của ∆k sẽ nhỏ dần khi k tăng, số hạng thứ hai

không nhất thiết nhỏ. Do đó ta có thể viết

||rk|| ≈ ||∆ku|| ≤ ||∆k|| ≤ ||M||||M−1K||.

Do đó, áp dụng trực tiếp phương pháp Arnoldi sẽ chỉ đạt được thặng dư ở mức ||∆ku||. Lí do là vì không gian Krylov không chứa thông tin về ma trận M−1K và do đó hình chiếu Gk sẽ không chứa đủ thông tin. Để khắc phục điều này, ta sẽ sử dụng một chiến lược trượt và nghịch đảo để giảm bớt ảnh hưởng

của ∆k.

Xét một giá trị riêng của Q(λ ) gần nhất với một giá trị σ nào đó. Ta chọn luôn λ0 = σ là xấp xỉ ban đầu của giá trị riêng cần tìm và x0 là véc tơ riêng

(cid:48)

tương ứng nào đó. Sử dụng phép biến đổi trượt và nghịch đảo

, = µ 1 λ (cid:48) = 1 λ − λ0

dẫn đến

Q(λ ) = M(λ (cid:48) + λ0)2 +C(λ (cid:48)+0) + K = Mλ (cid:48)2 + (2λ0M +C)λ (cid:48) + Q(λ0).

Tiếp theo, ta tìm giá trị lớn nhất theo mô đun của bài toán trượt và nghịch đảo

(cid:98)Q(µ) = µ 2Q(λ ) = (cid:98)Mµ 2 + (cid:98)Cµ + (cid:98)K,

với

(cid:98)M = Q(λ0), (cid:98)C = 2λ0M +C, (cid:98)K = M.

Giả sử λ0 (cid:54)= 0, ta có thể viết

. (cid:98)K = (cid:98)C −C 2λ0

Cuối cùng, áp dụng thuật toán Arnoldi trình bày ở phần trước vào (cid:98)Q(µ) và xấp xỉ (cid:98)Q(µ) bởi (cid:98)Qk(µ) = Iµ 2 + Hkµ + Gk với Hk là hình chiếu của (cid:98)M−1 (cid:98)C sinh bởi thuật toán Arnoldi và Gk là hình chiếu của (cid:98)M−1 (cid:98)K, tức là

k (cid:98)M−1

Gk = QT (cid:98)KQk.

k (cid:98)M−1CQk.

Gk = , với (cid:98)Gk = QT Khi đó, ta dễ dàng kiểm tra thấy (Hk − (cid:99)Gk) 2λ 0 Ta có kết quả sau đây về thặng dư của phép xấp xỉ.

0 . Khi đó

Mệnh đề 2.2.1. Giả sử µ0 là giá trị riêng có mô đun lớn nhất của (cid:98)Qk(µ) với u = (u1, ..., uk)T là một véc tơ riêng có độ dài đơn vị tương ứng. Đặt x1 = Qku và λ1 = λ0 + µ −1

1 +Cλ1 + K)x1|| ≤

1 − λ 2 λ 2 0 2|λ0|

||(Mλ 2 hk+1,kuk||Q(λ0)qk+1||

1 − λ 2 λ 2 0 2|λ0|

+ (2.6) ||∆ku||,

trong đó ∆k = CQk − Q(λ0)Qk (cid:99)Gk.

Chứng minh mệnh đề này được trình bày chi tiết trong [6].

Cũng giống như phần trước, số hạng thứ nhất của vế phải của (2.6) sẽ

giảm khi k tăng trong khi không có gì đảm bảo số hạng thứ hai sẽ giảm nếu

không có một chiến lược đặc biệt. Ở đây ta sẽ cập nhật biên độ của phép tịnh

tiến cho đến khi giá trị riêng tìm được thỏa mãn |λ − σ | ≤ η với η > 0 cho

trước. Điều đó dẫn đến thuật toán sau đây.

1: for l = 1, 2, · · · , đến khi hội tụ do

Algorithm 4 Lặp trượt và nghịch đảo Arnoldi. Require: σ , η, dự đoán x, số bước lặp tối đa m, đặt λ0 = σ ; Ensure: Một cặp riêng (µ, u) lớn nhất theo modul

x ||x||2 (cid:98)M = Q(λ0), (cid:98)C = 2λ0M +C, (cid:98)K = M, q1 = . for j = 1, 2, · · · , m do

(cid:98)q = (cid:98)M−1 (cid:98)Cq; for i = 1, · · · , j do

i (cid:98)q; (cid:98)q1 = (cid:98)q − qihi j; i ( (cid:98)M−1 (cid:98)Kq j); g ji = qT

i ( (cid:98)M−1 (cid:98)Kqi);

10:

11:

12:

hi j = qT gi j = qT

13:

; end for h j+1, j = ||(cid:98)q||2; if h j+1, j > 0, then q j+1 = (cid:98)q h j+1, j

14:

else

15:

16:

break.

end if Tính cặp riêng lớn nhất (µ0, u) của Iµ 2 + H jµ + G j = 0.

17:

18:

γ1 = h j+1, j|u j||| (cid:98)Mq j+1||;

19:

20:

|λ 2 (CQ j − (cid:98)MQ j(H j − 2λ0G j))|u|;

21:

λ1 = λ0 + µ −1 0 ; 1 − λ 2 |λ 2 0 | 2|λ0| 1 − λ 2 0 | γ2 = 2|λ0| if γ1 ≤ 0.1γ2 then

22:

break.

23:

24:

25:

end if

26:

end for x = Qku; if |λ1 − σ | < η then

27:

λ0 = λ1;

28: end for

end if

2.3 Ví dụ số

Trong mục này, chúng tôi sẽ sử dụng phương pháp Newton để tính cặp

riêng của bài toán giá trị riêng bậc hai. Để thực hiện, chúng tôi sử dụng phần

mềm MATLAB.

Để tiện so sánh, chúng tôi sử dụng lại bài toán ở ví dụ 1.3.1, tiểu mục

1.3.3.

Chúng tôi chọn x0 = [1 − 1 1]T , λ0 = −10, vs = [1 1 0]T với mọi s. Chúng

tôi sử dụng tiêu chuẩn dừng

|λs+1 − λs| < tol = 10−6

Và cuối cùng ma trận C được chọn là ma trận đơn vị.

Sau 8 bước lặp, chúng tôi thu được giá trị riêng

λ = −4.758617077184492.

So với giá trị riêng được tính trực tiếp bằng phần mềm MATLAP

λ1 = −4.758617077184495

thì đây là kết quả rất chính xác. Cũng với chương trình này, ta thu được véc tơ riêng tương ứng có giá trị rất lớn (khoảng 1021). Tuy nhiên, cũng bằng

MATLAB, chúng ta kiểm tra được rằng véc tơ này và véc tơ cột thứ nhất của

X trong ví dụ 1.3.1 ở tiểu mục 1.3.3 là phụ thuộc tuyến tính. Điều đó tức là

chương trình đã xác định đúng không gian riêng tương ứng. Để giảm giá trị của véc tơ, ta có thể sử dụng ma trận C với chuẩn khoảng 10−3. Khi đó, sau 8 bước lặp, chuẩn của véc tơ thu được sẽ giảm đi 10−24 lần và làm cho các

giá trị thành phần của véc tơ thuận mắt hơn. Cũng cần nói thêm rằng thuật

toán hội tụ đến giá trị riêng có trị tuyệt đối lớn nhất vì nó gần giá trị khởi động λ0 = −10 nhất. Khi chúng tôi chọn λ0 = −10−3, thuật toán hội tụ đến −0.026425971910232 là xấp xỉ của λ6 = −0.026425971912474 sau 4 bước

lặp.