Giáo trình Matlab cơ bản - Chương 2: Ma trận

seholder.m:  
clear all, clc 
global c d 
a = [ 1  2  3  4; 2  7  6  5; 3  6  9  0;4  5  0  3]; 
[eval, evec] = symmateig(a) 
n = size(a, 1); 
§31. THUẬT TOÁN RUTISHAUSER 
  Cho ma trận [A] đối xứng, xác định dương, ta dùng thuật toán Crout để 
phân tích nó thành: 
  [A] = [L][R]                   (1) 
Sau đó ta xây dựng ma trận: 
  [A1] = [R][L]                  (2) 
 125
Ma trận [A1] sẽ có cùng giá trị riêng với ma trận [A]. Ta lại tiếp tục phân tích 
ma trận [A1] theo thuật toán Crout: 
  [A1] = [L1][R1]                  (3) 
Và sau đó ta xây dựng ma trận: 
  [A2] = [R1][L1] 
Tiếp tục quá trình ta nhận được ma trận: 
  [An] = [Ln][Rn] 
  [An+1] = [Rn][Ln] 
Khi n → ∞,  [An]  trở  thành ma  trận  tam giác  trên có các phần  tử  trên đường 
chéo chính là các giá trị riêng của ma trận [A]. Phép lặp này hội tụ khi [A] là 
ma trận đố xứng, xác định dương. Khi [A] không xác định dương, phép  lặp 
không ổn định. Phân tích Crout cũng không thực hiện được khi có một phần 
tử rii = 0. Ta xây dựng hàm rutishauser() để thực hiện thuật toán trên: 
function [lambda, evec] = rutishauser(a) 
% Tinh cac gia tri rieng bang thuat toan phan tich LR 
for i = 1: 500 
    s1 = diag(a); 
    [l, r] = crout(a); 
    a = r*l; 
    s2 = diag(a); 
    if abs(max(s2 ‐ s1)) < 1e‐8 
       break; 
   end 
end 
lambda = diag(a); 
n = size(lambda, 1); 
evec = zeros(n); 
c = evec; 
e = eye(n); 
for i = 1:n 
    b = a ‐ lambda(i)*e; 
    c = nulld(b); 
    evec(:, i) = c(:, 1);   
end 
 126
lambda = diag(lambda); 
Để tính các giá trị riêng của ma trận ta dùng chương trình ctrutishauser.m: 
clear all, clc 
a = [ 1 2 3 4; 5 6 7 8; 9 0 1 2; 3 4 5 6]; 
[lambda, evec] = rutishauser(a) 
§32. THUẬT TOÁN FADDEEV – LEVERIER ‐ SOURIAU 
  Cho ma trận [A]: 
  [ ]
⎡ ⎤⎢ ⎥⎢ ⎥= ⎢ ⎥⎢ ⎥⎣ ⎦
L
L
M M L M
L
11 11 1n
21 22 2n
n1 n2 nn
a a a
a a a
A
a a a
Các giá trị riêng của ma trận [A] là nghiệm của đa thức đặc trưng: 
  [ ] [ ]( )− λ =det A E 0                 (1) 
Khi khai triển định thức trên, ta được đa thức cấp n: 
  − −λ = λ − λ − λ − − =Ln n 1 n 1n 1 1 nP ( ) p p p 0           (2) 
Gọi vết của ma trận là tổng các phần tử trên đường chéo chính: 
  + + +L11 22 nntrace([A]) = a a a             (3) 
Các hệ số của đa thức (2) được xác định theo: 
  [ ] [ ] [ ]= =1 1 1p trace( B ) B A  
  [ ] [ ] [ ] [ ] [ ]( )= = −2 1 2 1 11p trace( B ) B A B p E2  
[ ] [ ] [ ] [ ] [ ]( )= = −3 3 3 2 21p trace( B ) B A B p E3  
[ ] [ ] [ ] [ ] [ ]( )− −= = −
L
n n n n 1 n 1
1p trace( B ) B A B p E
n
Giải hệ (3) ta có các giá trị riêng của [A]. 
Ta xây dựng hàm fadlev() để thực hiện thuật toán trên: 
function [lambda, evec] = fadlev(a) 
%Tim cac gia tri rieng va vec to rieng bang phuong phap Faddeev ‐ Leverrier 
n = size(a, 1); 
b = a; 
 127
for i = 1:n 
    s = 0; 
    for k = 1:n 
        s = s + b(k, k); 
    end 
    p = s/i; 
    b = a*(b ‐ p*eye(n)); 
    r(i+1) = ‐p; 
end 
r(1) = 1; 
lambda = roots(r); 
n = size(lambda, 1); 
evec = zeros(n); 
c = evec; 
e = eye(n); 
for i = 1:n 
    b = a ‐ lambda(i)*e; 
    c = nulld(b); 
    evec(:, i) = c(:, 1);   
end 
lambda = diag(lambda); 
Để tìm các giá trị riêng của ma trận ta dùng chương trình ctfaclev.m: 
clear all, clc 
a = [11 2 3 1 4; 2 9 3 5 2; 4 5 6 7 8; 3 2 1 6 7; 4 9 6 8 2]; 
[lambda, evec] = fadlev(a) 
§33. THUẬT TOÁN KRYLOV 
Đa thức đặc trưng của ma trận [A] có dạng: 
  − −λ = λ − λ − λ − − =Ln n 1 n 2n 1 1 nP ( ) p p p 0           (1) 
Ta viết lại (1) dưới dạng: 
−
=
λ = λ + λ∑n 1n in i
i 0
P ( ) p                 (2) 
Theo định lý Hamilton ‐ Kelly ta có P([A]) = 0. Xét dãy lặp: 
 128
[ ]+ =(i 1) (i)v A v   i = 0, 1,..., n‐1            (3) 
Do P([A]) = 0 nên: 
  [ ]( ) −
=
ρ = + ρ =∑n 1(0) (n) (i)i
i 0
A v v v 0              (4) 
Đặt  ⎡ ⎤= ⎣ ⎦
T(i) (i) (i)
1 nv v ,...,v  ta có: 
−
=
+ ρ =∑n 1(n) (i)j i
i 0
v v 0                  (5) 
hay: 
−
−
−
⎡ ⎤ ⎡ ⎤ρ⎡ ⎤⎢ ⎥ ⎢ ⎥⎢ ⎥ = −⎢ ⎥ ⎢ ⎥⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ρ⎣ ⎦⎣ ⎦ ⎣ ⎦
L
M L M M M
L
(n 1) (0) (n)
1 1 n 1 1
(n 1) (0) (n)
n n 0 n
v v v
v v v
            (6) 
Vì  [ ]+ =(k 1) (k)v A v  
nên: 
[ ]( )+
=
= =∑n(k 1) (k) (k)i i ,j ji
j 1
v A v a v   i = 1, 2,..., n    k = 0, 1,..., n‐1    (7) 
Tóm lại quá trình tính toán theo thuật toán Krylov như sau: 
  ‐ Chọn v(0) tuỳ ý, tính  (k)iv  theo (7) 
  ‐ Giải hệ (6) để tìm ρ 
Khi đã tìm được ρ  ta có đa thức đặc trưng.  
Ta xây dựng hàm krylov() để thực hiện thuật toán trên: 
function [lambda, evec] = krylov(a) 
% Tim gia tri rieng bang thuat toan Krylov 
n = size(a, 1); 
v = zeros(n, 1); 
v(1) = 1; 
b = zeros(n, n); 
b(:, n) = v; 
for k = 2:n 
    v = a*v; 
    b(:, n‐k+1) = v; 
end 
c = a*v; 
 129
rho = ‐b\c; 
rho = [1 rhoʹ]; 
lambda = roots(rho); 
n = size(lambda, 1); 
evec = zeros(n); 
c = evec; 
e = eye(n); 
for i = 1:n 
    b = a ‐ lambda(i)*e; 
    c = nulld(b); 
    evec(:, i) = c(:, 1);   
end 
lambda = diag(lambda); 
Để tìm giá trị riêng của ma trận ta dùng chương trình ctkrylov.m: 
clear all, clc 
a = [17 24 30 17; 8 13 20 7; 2 10 8 6; ‐23 ‐43 ‐54 ‐26]; 
[lambda, evec] = krylov(a) 
§34. THUẬT TOÁN HYMAN 
  Cho ma trận Hessenberg [B] và λ là một giá trị riêng của [B]. Ta tìm: 
  α(λ), x1(λ), x2(λ),...,xn‐1(λ) 
sao cho [X] = [x1(λ), x2(λ),...,xn‐1(λ), xn]T với xn = 1 là nghiệm của bài toán: 
  [ ] [ ]( )[ ] [ ]1B E X e− λ = α  
nghĩa là: 
1,1 1 1,2 2 1,n
2,1 1 2,2 2 2,n
n,n 1 n 1 n,n
(b )x b x b 1
b x (b )x b 1 0
b x (b )1 0− −
− λ + + + = α⎧⎪ + − λ + + =⎪⎨⎪⎪ + − λ =⎩
L
L
M  
Thay thế ngược ta có được các nghiệm xn‐1, xn‐2,..., x1 và α(λ).  
Theo quy tắc Cramer ta có: 
 130
[ ] [ ]( )
1,1
2,1
n 1,n 1
n,n 1
n
b
b 0
det
0 b 0
0 b 0
1 x
det B E
− −
−
− λ ∗ α⎡ ⎤⎢ ⎥∗⎢ ⎥− λ⎢ ⎥⎢ ⎥⎣ ⎦= = − λ  
            [ ] [ ]( )
n 1
2,1 n ,n 1( )( 1) b b
det B E
−
−α λ −= − λ
L
Do đó: 
  [ ] [ ]( )n 1
2,1 n ,n 1
( 1)( ) det B E
b b
−
−
−α λ = − λL  
Đa thức α(λ) chỉ sai khác đa thức đặc trưng p(λ) một nhân tử hằng. Do vậy ta 
có thể dùng α(λ) để tìm các giá trị riêng thay cho p(λ). Hàm alffa() thực hiện 
công việc này:  
function r = alfa(a); 
n = size(a, 1); 
b = a; 
l = 1; 
for i = 2:n 
    l = l*a(i, i‐1); 
end 
sign = (‐1)^(n‐1); 
for i = 1:n 
    s = 0; 
    for k = 1:n 
        s = s + b(k, k); 
    end 
    p = s/i; 
    b = a*(b ‐ p*eye(n)); 
    r(i+1) = ‐p; 
end 
r(1) = 1; 
r = sign*r/l; 
 131
Nếu cho một ma trận bất kì [A], ta có thể biến đổi Householder nó thành ma 
trận Hessenberg [B] đồng dạng với [A]. Sau khi có các giá trị riêng, ta tìm các 
vec tơ riêng tương ứng. Hàm hyman() thực hiện thuật toán này: 
function [lambda, evec] = hyman(a) 
%Tim cac gia tri rieng va vec to rieng bang phuong phap Hyman 
b = hessenberg(a); 
r = alfa(b); 
lambda = roots(r); 
n = size(lambda, 1); 
evec = zeros(n); 
c = evec; 
e = eye(n); 
for i = 1:n 
    b = a ‐ lambda(i)*e; 
    c = nulld(b); 
    evec(:, i) = c(:,1);   
end 
lambda = diag(lambda); 
Để  tìm các giá trị riêng và vec tơ riêng tương ứng của một ma trận ta dùng 
chương trình cthyman.m: 
clear all, clc 
a = [ 1 2 3 4; 2 5 2 3;7 1 4 1; 3 2 3 7]; 
[lambda, evec] = hyman(a) 
§35. TRỰC GIAO HOÁ ARNOLDI 
  Trong đại số tuyến tính, phương pháp Arnoldi được W. Arnoldi đưa ra 
năm 1951. Phương pháp lặp Arnoldi sử dụng trực giao hoá Gram – Schmidt 
để tạo ra dãy các vec tơ trực giao [q1], .., [qn] gọi là các vec tơ Arnoldi. Thuật 
toán cụ thể gồm các bước: 
- Cho vec tơ bất kì [q1] có  [ ]1q 1=  
- Lặp từ k = 2,3,... 
 132
• [A][qk‐1] = [qk] 
• for j = 1:k‐1 
  ∗  Tj k j,k 1q q h −⎡ ⎤⎡ ⎤ =⎣ ⎦ ⎣ ⎦  
∗  k j,k 1 jq h q−⎡ ⎤ ⎡ ⎤= ⎣ ⎦⎣ ⎦  
•  [ ]k k,k 1q h −=  
• [ ] kk
k,k 1
qq
h −
=  
Ta gọi [Qn] là ma trận tạo từ n vec tơ Arnoldi đầu tiên [q1],..,[qn] đầu tiên và 
[Hn] là ma trận (Hessenberg trên) tạo bởi n hàng đầu tiên của [H] và có: 
  [Hn] = [Q]T[A][Q]  
Ta xây dựng hàm arnoldi() thực hiện thuật toán trên:  
function [Hn, Vn] = arnoldi(A); 
k = 50; 
m = size(A, 1); 
v = rand(m, 1); 
n = size(A, 1); 
H = zeros(m+1, m); 
V = zeros(n, m+1); 
V(:,1) = v/norm(v); 
for j = 1:m 
    w = A*V(:, j); 
    for i = 1:j 
        H(i, j) = V(:, i)ʹ*w; 
        w = w ‐ H(i, j)*V(:, i); 
    end 
    H(j+1, j) = norm(w); 
    if H(j+1, j) <= eps, 
       break; 
    end 
    V(:, j+1) = w/H(j+1, j); 
end 
Hn = H(1:m, :); 
 133
Vn = V(:, 1:m); 
Để phân tích ma trận ta dùng chương trình ctarnoldi.m: 
clear all, clc 
A = [ 7 2 3 ‐1; 2 8 5 1; 3 5 12 9; ‐1 1 9 7];  
[H, V] = arnoldi(A) 
§36. TRỰC GIAO HOÁ LANCZOS 
  Cho ma trận [A] đối xứng. Ta phân tích ma trận thành các ma trận [Q] 
và [T] sao cho: 
  [A] = [Q][T][Q]T 
với:  [Q][Q]T  =  [E],  nghĩa  là  [Q]  là ma  trận  trực  giao  và  [T]  là ma  trận  ba 
đường chéo đối xứng. Thuật toán Lanczos gồm các bước: 
- Cho vec tơ [v1] có  [ ]1v 1=  
- Cho v0 = 0, β0 = 0. [V] = [v1] 
- Lặp : 
•  j 1 j j 1 j 1v A v v+ − −⎡ ⎤ ⎡ ⎤ ⎡ ⎤⎡ ⎤= −β⎣ ⎦⎣ ⎦ ⎣ ⎦ ⎣ ⎦  
•  Tj j 1alfa v v +⎡ ⎤ ⎡ ⎤= ⎣ ⎦ ⎣ ⎦  
•  j 1 j 1 jv v v+ +⎡ ⎤ ⎡ ⎤ ⎡ ⎤= − α⎣ ⎦ ⎣ ⎦ ⎣ ⎦  
•  Tj 1 j 1 j 1v v V V v+ + +⎡ ⎤ ⎡ ⎤⎡ ⎤ ⎡ ⎤ ⎡ ⎤= −⎣ ⎦ ⎣ ⎦ ⎣ ⎦⎣ ⎦ ⎣ ⎦  
•  j j 1v +⎡ ⎤β = ⎣ ⎦  
•  j 1 j 1 jv v+ +⎡ ⎤ ⎡ ⎤= β⎣ ⎦ ⎣ ⎦  
• [ ] j 1V V,v +⎡ ⎤= ⎣ ⎦  
Ta xây dựng hàm lanczos() để thực hiện thuật toán Lanczos 
function [T, Q] = lanczos(A); 
% thuat toan Lanczos cho ma tran doi xung  
n = size(A, 1); 
Q(:, 1) = rand(n, 1); 
Q(:, 1) = Q(:, 1)./norm(Q(:, 1)); 
a(1) = Q(:, 1)ʹ*A*Q(:, 1); 
Q(:, 2) = A*Q(:, 1) ‐ a(1)*Q(:, 1); 
 134
b(1) = norm(Q(:, 2)); 
Q(:, 2) = 1./b(1)*Q(:,  2); 
for i = 2:n‐1 
    a(i) = Q(:, i)ʹ*A*Q(:, i); 
    Q(:, i+1) = A*Q(:, i) ‐ a(i)*Q(:, i) ‐ b(i‐1)*Q(:, i‐1); 
    b(i) = norm(Q(:, i+1)); 
    Q(:, i+1) = 1./b(i)*Q(:, i+1); 
end 
a(n) = Q(:, n)ʹ*A*Q(:, n); 
T = diag(a) + diag(b, 1) + diag(b, ‐1); 
Để phân tích ma trận ta dùng chương trình ctlanczos.m: 
clear all, clc 
A = [4 1 3 ‐2; 1 ‐2 4 1; 3 4 1 2; ‐2 1 2 3]; 
[T, Q] = lanczos(A);