Bài giảng phần Ma trận

            pp = max(size(p));            v = nulld(aa, small);                     jv(:,p) = v*(rand(size(v, 2), pp) ‐ 0.5)*16;            k = k + pp;        end        clear v;        for k = 1:dns(1)             v(:,1) = jv(:, k);            for p = 2:ec(k)                v(:, p) = a1*v(:, p‐1);            end            vv = fliplr(v(:, 1:ec(k)));            M = [M vv];  84       end     end  end  k = abs(det(M))^(‐1/n);  M = k*M;   Mi = inv(M);  D = Mi*a*M;  d0 = diag(D);  d1 = diag(D, 1);  D = diag(d0) + diag(d1, 1);  function Z = nulld(A, small)  norma = sqrt(mean(mean(abs(A))));  tiny = norma*eps;  if nargin<2     small = (1 + norma)*sqrt(eps);  end  [m, n] = size(A);  if  m~= n     error(ʹMa tran phai vuong!ʹ)  end  p = find(abs(A)<tiny);  if ~isempty(p)     A(p) = 0;  end  [U,S,V] = svd(A, 0);  S = diag(S);  s = S;  norma = max(s);  smax = max(s);  if smax == 0;      smax = 1;  end  s = s/smax;  snorm = s;  85 t = find(s>0);  if isempty(t);      Z = V;      return;  end  p = find(s<tiny);  if ~isempty(p)     s(p) = tiny;  end  p = find(s == 0);  if ~isempty(p)     s(p) = small*min(s(t));  end  logs = log10(s);  sdifflog = ‐diff(logs);  smax = max(sdifflog);  r = find(sdifflog == smax);  if min(size(r))>0     r = r(end);  else     r = 0;  end  Z = V(:,r+1:n);  if snorm == ones(n, 1);        Z = zeros(n, 0);   end  if max(S) <= small;       Z = V;   end  §13. PHÂN TÍCH MA TRẬN THEO CÁC GIÁ TRỊ KÌ DỊ     Phân tích ma trận theo các giá trị kì dị (kì dị value) được thực hiện trên  các ma trận vuông hay chữ nhật. Ta có:    [Anp] = [Unn][Snp][Vpp]  Trong đó:  86   [U]T[U] = [Enn]    [V]T[V] = [Epp]  nghĩa là các ma trận [U] và [V] là trực giao.  Các cột của  [U]  là các vec  tơ kì dị  trái,  [S] có các giá  trị kì dị và  là ma  trận  đường chéo và [V]T có các hàng là các vec tơ kì dị phải. Để tính các ma trận  [U], [S] và [V] ta tìm các giá trị riêng của [A][A]T và [A]T[A]. Các vec tơ riêng  của [A]T[A] tạo nên các cột của [V]. Các vec tơ riêng của [A][A]T tạo nên các  cột của  [U]. Các giá  trị kì dị của  [S]  là căn bậc hai của các giá  trị  riêng của  [A][A]T hay  [A]T[A]. Các giá  trị riêng  trên đường chéo của  [S] được sắp xếp  theo thứ tự giảm dần. Để hiểu được thuật toán này ta xét ví dụ sau:  Ta xây dựng hàm svddecomp() để thực hiện thuật toán này:  function [U, S , V] = svddecomp(A)  [m, n] = size(A);  if (m > n)      % ta can timcac vec to kì dị phai truoc      %[V D] = eigs(Aʹ*A)      [V, D] = eigs(Aʹ*A);      [dummy, perm] = sort(‐diag(D));      S = diag(sqrt(diag(D(perm, perm))));      V = V(:, perm);      sinv = diag(1./sqrt(diag(D)));      U = (A*V)*sinv;  else      % ta can tim cac vec to kì dị trai truoc      % [U D] = eigs(A*Aʹ)      [U, D] = eigs(A*Aʹ);      [dummy, perm] = sort(‐diag(D));      S = diag(sqrt(diag(D(perm, perm))));      U = U(:, perm);      sinv = diag(1./sqrt(diag(D)));      V = sinv*(Uʹ*A);      V = Vʹ;  end  87 Để phân tích một ma trận ta dùng chương trình ctsvddecomp.m:  clear all, clc  %a = [ 1 2 3 4 5; 6 7 8 9 0; 3 4 5 6 7; 8 9 0 1 2; 2 4 6 8 1];  a = [ 1 1; 0 1; 1 0];  [u, s, v] = svddecomp(a)  §14. PHÂN TÍCH SCHUR    Cho ma trận vuông [A], cấp n ta phân tích nó thành:    [A] = [T][U][T]*  Trong đó:  [T]  ‐ ma  trận  unita  và  [T]*  là ma  trận  chuyển  vị  liên  hợp  của  [T](lâý  chuyển vị của [T] rồi lấy liên hợp của các phần tử).  [U] = [Λ] + [N]  [Λ] ‐ là ma trận đường chéo có các phần tử là các giá trị riêng của [A]  [N]  ‐ ma  trận  tam giác phải,  có  các phần  tử  thuộc  đường  chéo  chính  bằng zero.  Mọi ma trận vuông đều có thể phân tích Schur. Tuy nhiên phân tích này  không duy nhất. Nếu [A] là ma trận trực giao thì [U] là ma trận đường chéo  và các cột của [T] là các vec tơ riêng của [A]. Phân tích Schur khi này được gọi  là phân tích phổ. Nếu [A] xác định dương, phân tích Schur chính là phân tích  SVD.  Để phân tích ma trận theo thuật toán Schur ta thực hiện các bước sau:  - Tìm giá trị riêng λ1 của [A] và vec tơ riêng tương ứng [v1]  - Chọn n‐ 1 vec tơ [w1],...,[wn‐1] độc lập tuyến tính và trực giao với [v1]  - Tạo [V1] là ma trận có các cột là [v1], [w1],...,[wn‐1] và tính:  [ ] [ ][ ] [ ]11 1 1V A V 0 A ∗ λ ∗⎡ ⎤= ⎢ ⎥⎣ ⎦ Trong đó [A1] là ma trận (n‐1)×(n‐1)  - Lặp lại quá trình với ma trận [A1] ta có:       [ ] [ ][ ] [ ]22 1 2 2V A V 0 A ∗ λ ∗⎡ ⎤= ⎢ ⎥⎣ ⎦ Trong đó [A2] là ma trận (n‐2)×(n‐2)  88 - Do [ ] [ ] [ ] [ ] 1 2 2 1 2 T A T 0 0 0 A ∗ ∗ λ ∗ ∗⎡ ⎤⎢ ⎥= λ ∗⎢ ⎥⎢ ⎥⎣ ⎦ Trong đó [ ]2 1 2ˆT V V⎡ ⎤ ⎡ ⎤= ⎣ ⎦ ⎣ ⎦  với  [ ]2 2 1 0 Vˆ 0 V ⎡ ⎤⎡ ⎤ = ⎢ ⎥⎣ ⎦ ⎣ ⎦ - Tiếp tục quá trình ta tìm được [V1],...,[Vn]. Cuối cùng [U] = [T]*[A][T]   2 1 2 n ˆ ˆT V V V⎡ ⎤ ⎡ ⎤ ⎡ ⎤⎤⎡ =⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦L   Ta xây dựng hàm schurdecom() thực hiện thuật toán trên:  function [T, U] = schurdecom(a)  % Phan tich Schur ma tran A   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;  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);       89 end  p = grams(evec);  T = conj(transpose(p))*a*p;  U = p;  Để phân tích ma trận ta dùng chương trình ctschur.m:  clear all, clc  a = [ 1 2 3 5; 4 5 6 2; 4 6 8 9; 9 3 6 7];  [t, u] = schurdecom(a)  §15. ĐỊNH THỨC CỦA MA TRẬN  Cho một ma trận vuông cấp n. Ta cần tìm định thức của nó. Trước hết  chúng ta nhắc lại một số tính chất quan trọng của định thức:  - nếu nhân  tất cả các phần  tử của một hàng  (hay cột) với k  thì  định  thức được nhân với k  - định  thức không đổi nếu  ta cộng  thêm vào một hàng  tổ hợp  tuyến  tính của các hàng còn lại.  - nếu đổi chỗ hai hàng cho nhau thì định thức đổi dấu  Trước khi đi đến định nghĩa về định thức ta tìm hiểu khái niệm về hoán  vị và phép thế.   Cho một dãy số, nếu ta đổi chỗ các số trong dãy cho nhau thì ta đã thực  hiện một phép hoán vị. Ví dụ 123, 132,..  là các hoán vị của dãy số  {1, 2, 3}.  Trong hoán vị α1α2…αi…αj…αn  ta nói αi  làm một nghịch  thế với αj nếu  i <  j  mà αi > αj. Ví dụ trong hoán vị 1432 số 4 làm với số  3 một nghịch thế  , số 4  làm với số 2 một nghịch thế, số 3 làm với số 2 một nghịch thế. Một hoán vị gọi  là chẵn nếu tổng số nghịch thế trong hoán vị đó là một số chẵn; một hoán vị  gọi là lẻ trong trường hợp ngược lại. Như vậy  1432 là một hoán vị lẻ.  Cho một dãy  số, nếu  ta  tạo  ra một dãy  số mới bằng cách  đổi chỗ các  phần tử cho nhau thì ta đã thực hiện một phép thế.    Ví dụ  2 1 4 3 p 1 4 2 3 ⎛ ⎞= ⎜ ⎟⎝ ⎠  là phép thế biến 2 thành 1, 1 thành 4,  4 thành 2 và 3  thành 3.  Một phép thế gọi là chẵn nếu tính chẵn lẻ của dòng trên và dòng dưới  như nhau và lẻ trong trường hợp ngược lại. Phép thế trên là phép thể lẻ.  90 Cho ma  trận  vuông  [A]  cấp  n.  Các  phần  tử  của  hàng  thứ  i  là  ai,1,  ai,2,…,ai,n. Các phần tử của cột thứ j là a1,j, a2,j ,…, an,j. Ta xem hàng thứ i là một  vec  tơ, kí hiệu  là Ai* và cột  thứ  j cũng  là một vec  tơ, kí hiệu  là A*j. Với mỗi  phép thế:  1 2 n 1 1 n i i i p j j j ⎛ ⎞= ⎜ ⎟⎝ ⎠ L L                 (1)  ta lập tích:  1 1 2 2 n ni j i j i j a a aK                   (2)  Trước mỗi tích (2) ta đặt dấu + nếu  và dấu ‐ nếu phép thế (1) lẻ. Sau đó ta lập  tổng của n! tích có dấu như vậy, nghĩa là tổng:  1 1 2 2 n n t(p) i j i j i j p ( 1) a a a−∑ K                 (3)  trong đó:    t(p) = 1 nếu phép thế p lẻ    t(p) = 0 nếu phép thế p chẵn  Tổng (4) được gọi là định thức của ma trận vuông [A], cấp n.  Ta  xây  dựng  hàm  determinant()  để  tính  định  thức  của ma  trận  theo  định  nghĩa:  function d = determinant(A)  % DETERMINANT tinh dinh thuc theo dinh nghia.  [m, n] = size(A);  if ( m ~= n )      fprintf ( ʹ\nʹ );        fprintf ( ʹ Chi ma tran vuong moi co dinh thuc!\nʹ );       return  end  p = zeros(1, n);  nf = prod([1:n]);  d = 0.0;  for i = 1:nf      p = nextperm(p);      s = permsign(p);      x = diag(A([1:n],p));  91     d = d + s*prod(x);  end  function psign = permsign(p)  % PERMSIGN tra ve dau phep the   .  %    +1, neu phep the chan,  %    ‐1, neu phep the le.  n = length ( p );  psign = 1;  for i = 1:n‐1      j = i;      while (p(j) ~= i)          j = j + 1;      end      if ( j ~= i )          temp = p(i);          p(i) = p(j);          p(j) = temp;            psign = ‐ psign;      end  end  function q = nextperm(p)  n = length(p);  q = p;  if(n == 1)       q = 1;  elseif (q == 0)       q = [1:n];  else      i = n ‐ 1;      while (q(i) > q(i+1))          i = i ‐ 1;          if (i == 0)               break;  92         end      end          if (i == 0)           q = [1:n];             else          j = n;          while (q(j) < q(i))              j = j ‐ 1;          end                  t = q(j);          q(j) = q(i);          q(i) = t;                 q(i+1:n) = q(n:‐1:i+1);      end  end  Để tính định thức ta dùng chương trình ctdeterminant.m:  clear all, clc  %a = [1  2; 3  5];  a = [1  3  5; 3  4  6; 4  6  3];  d = determinant(a)  §16. TÍNH ĐỊNH THỨC BẰNG CÁCH PHÂN TÍCH MA TRẬN  Cho ma trận [A]. Nếu  [A] = [B]×[C]  thì   det[A] = det[B]×det[C]  Mặt khác với một ma trận tam giác, ví dụ:    [ ] 11 12 13 14 22 23 24 33 34 44 b b b b 0 b b b B 0 0 b b 0 0 0 b ⎡ ⎤⎢ ⎥⎢ ⎥= ⎢ ⎥⎢ ⎥⎣ ⎦ thì  93   [ ] 11 22 33 44det B b b b b=   nghĩa là đối với ma trận tam giác, định thức bằng tích các phần tử trên đường  chéo chính.     Khi  phân  tích ma  trận  [A]  theo  thuật  toán Doolitte,  ta  dùng  chương  trình ctdoodecmp.m để tính định thức của nó:  clear all, clc  a = [1 2 3 4; 3 4 5 7; 2 3 8 5; 4 9 1 4];  [l, r] = doolittle(a);  d = prod(diag(l))*prod(diag(r))             Khi phân  tích  ma trận [A] theo thuật toán Crout, ta dùng chương trình  ctcrotdecmp.m để tính định thức của nó:  clear all, clc  a = [1 2 3 4; 3 4 5 7; 2 3 8 5; 4 9 1 4];  [l, r] = crout(a);  d = prod(diag(l))*prod(diag(r))  Khi phân  tích ma  trận  [A]  theo  thuật  toán Choleski,  ta dùng  chương  trình ctcholdecmp.m để tính định thức của nó:  clear all, clc  a = [4 ‐2 2;‐2 2 ‐4;2 ‐4  11];  a = pascal(5);  l = choleski(a);  d = prod(diag(l))*prod(diag(lʹ))  §17. THUẬT TOÁN TRỤ CHIÓ  Cho ma trận [A] có  1,1a 0≠ . Ta xây dựng ma trận [B] có các phần tử   i ,j 1,1 i ,j i ,n n ,jb a a a a= −    thì:    [ ] [ ]2 n1,1A a B−=   nghĩa là:  94 [ ] 11 12 11 13 11 1n 21 22 21 23 21 2n 11 12 11 13 11 1n n 2 31 32 31 33 31 3n11 11 12 11 13 11 1n n1 n2 n1 n3 n a a a a a a det det det a a a a a a a a a a a a det det det a a a a a adet A a det a a a a a a det det det a a a a a − ⎡ ⎤ ⎡ ⎤ ⎡ ⎤⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤⎢ ⎥ ⎢ ⎥ ⎢ ⎥= ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ⎡ ⎤⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎣ ⎦ L L M M L M L 1 nna ⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎡ ⎤⎢ ⎥⎢ ⎥⎣ ⎦⎣ ⎦ Ta xây dựng hàm chiopivot() để thực hiện thuật toán trên:  function d = chiopivot(a)  %tinh dinh thuc bang thuat toan Chio pivotal condensation  if a(1, 1) == 0      error(ʹKhong dung phuong phap  nay  de tinh dinh thuc duoc !ʹ);      return;  end  c = a(1, 1);  n = size(a, 1);  if (n <= 2)      d = a(1, 1)*a(2, 2) ‐ a(2, 1)*a(1, 2);      return  end  m = n ‐ 1;  while (m >= 1)      for i = 1:m%hang          b(i, 1:m) = a(1, 1)*a(i+1, 2:m+1) ‐ a(i+1, 1)*a(1, 2:m+1);      end      if (m > 2)          a = b;          c = c*a(1,1);          clear b;      end            m = m ‐ 1;  end  d = b(1, 1)*b(2, 2) ‐ b(2, 1)*b(1, 2);  95 d = d/c;  Để tính định thức ta dùng chương trình ctchiopivot.m:  clear all, clc  a = [1 2 3 4; 3 4 5 7; 2 3 8 5; 4 9 1 4];  d = chiopivot(a)  §18. THUẬT TOÁN LAPLACE  Để tính định thức theo thuật toán Laplace, ta khai triển định thức theo  hàng hay cột. Cho ma trận vuông [A] cấp n. Nếu bỏ đi hàng i và cột j (tức xoá  hàng và cột chứa phần tử aij) thì ta có một ma trận cấp (n ‐ 1), định thức của  nó gọi là định thức con cấp (n ‐ 1) ứng với phần tử aij (minor) , ký hiệu là Mij.  Ta chú ý đến hàng  thứ  i và aij  là một phần  tử của hàng đó. Trong det[A]  ta  gộp những số hạng chứa aij  lại và đặt aij  làm thừa số chung, hệ số của nó kí  hiệu là Aij và gọi là phần bù đại số (cofactor) của phần tử aij. Cofactor Aij của  ma trận [A] là:    i jij ijA ( 1) M += −   Định thức của [A] khi khai triển theo hàng là:    [ ] n ij ij i 1 det A a A = =∑   Ta xây dựng hàm cofactor() để tính các phần bù đại số:  function c = cofactor(A, i, j)  % cac phan bu dai so cua ma tran  % dung de nghich dao mt  % C = cofactor(A, i, j) tra ve phan bu dai so cua  %ng i, cot j cua A.  if nargin == 3      M = A;      M(i,:) = [];      M(:,j) = [];      c = (‐1)^(i+j) * det(M);  else  96     [n, n] = size(A);      for i = 1:n          for j = 1:n              c(i,j) = cofactor(A, i, j);          end      end  end  Sau khi phần bù đại số, ta xây dựng hàm cofactordet() để tính định thức của  [A]  function d = cofactordet(A)  d = 0;  for i = 1:size(A, 1)      c = cofactor(A, i, 1);      d = d + A(i, 1)*c;  end  Để tính định thức ta dùng chương trình ctcofactordet.m:  clear all, clc  a = [1 2 3 4; 3 4 5 7; 2 3 8 5; 4 9 1 4];  det = cofactordet(a)  §19. THUẬT TOÁN DODGSON  Thuật toán cô đặc Dodgson ( Dodgson condensation) dùng để tính định  thức  của ma  trận vuông. Nó  được Charles Ludwidge Dodgson  đưa  ra.  Để  tính định thức của ma trận cấp n × n, ta xây dựng các ma trận cấp (n ‐ 1) × (n ‐  1) cho đến ma trận cấp 1 × 1 là định thức của ma trận cần tìm.    Bước đầu tiên ta xây dựng ma trận cấp (n ‐ 1)×(n ‐ 1)  từ các định thức  của các ma trận con 2×2. Ví dụ với ma trận   5 1 0 2 8 5 0 6 7 ⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦ 97 ta có các ma trận con là:  5 1 1 0 2 8 8 5 2 8 8 5 0 6 6 7 ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦   Các ma trận này sẽ tạo ra các phân tử của ma trận 2×2. Phần tử ử hàng r, cột c  là định thức của ma trận con 2×2 của ma trận ban đầu với hàng r và cột c ở  góc trên trái. Như vậy ma trận mới là:  38 5 12 26 ⎡ ⎤⎢ ⎥⎣ ⎦   Ma trận k×k được tạo ra bằng cách lấy định thức của ma trận con 2×2 của ma  trận  (k+1)×(k+1) như ở  trên và chia nó cho phần  tử  tương ứng của ma  trận  trung tâm, nghĩa là ma trận đã bỏ đi hàng trên cùng, hàng dưới cùng, cột bên  phải và cột bên trái của ma trận (k+2)×(k+2)  Ta xây dựng hàm dodgson() để thực hiện thuật toán trên:  function dt = dodgson(a)  if size(a, 1) ~= size(a, 2)      error(ʹMa tran A phai la mt tran vuongʹ);  end;  n = size(a, 1);  if n == 2      dt = a(1, 1)*a(2, 2) ‐ a(2, 1)*a(1, 2);      return  end;  if n == 3;      for i = 1:n‐1          b(i, 1:n‐1) = a(i, 1:n‐1).*a(i+1, 2:n) ‐ a(i+1, 1:n‐1).*a(i, 2:n);             end      dt = (b(1, 1)*b(2, 2) ‐ b(2, 1)*b(1, 2))/a(2,2);      return  end  c = a;  c(:, 1) = [];  c(:, n‐1) = [];  c(1, :) = [];  98 c(n ‐ 1, :) = [];  for i = 1:n‐1      b(i, 1:n‐1) = a(i, 1:n‐1).*a(i+1, 2:n) ‐ a(i+1, 1:n‐1).*a(i, 2:n);       end  m = size(b, 1);  while m >= 2      for i = 1:m‐1          for j = 1:m‐1              d(i, j) = (b(i, j)*b(i+1, j+1) ‐ b(i+1, j)*b(i, j+1))/c(i, j);                    end      end      if m > 3          c = b;          c(:, 1) = [];          c(:, m‐1) = [];          c(1, :) = [];          c(m ‐ 1, :) = [];          b = d;      end      m = m ‐ 1;  end  dt = (d(1, 1)*d(2, 2) ‐ d(2, 1)*d(1, 2))/b(2,2);  Để tính định thức ta dùng chương trình ctdodgson.m:  clear all, clc;  a = [1 2 3 4; 5 1 3 2; 4 9 2 2; 6 3 4 1];  dt = dodgson(a)  §20. NGHỊCH ĐẢO MA TRẬN BẰNG CÁCH DÙNG MINOR    Cho ma trận [A], ta có:    ( ) [ ]− = j,i1 i ,j A a det A Trong đó:  99   ( )−1 i ,ja  là phần tử ở hàng i, cột j của ma trận [A]‐1    Ai,j là phần bù đại số của phần tử ai,j của ma trận [A]  Ta xây dựng hàm minorinv() để thực hiện thuật toán trên:  function c = minorinv(a)  % Tim ma tran nghich dao bang thuat toan minor  n = size(a, 1);  ms = det(a);  for i = 1:n      for k = 1:n         b = cofactor(a, i, k);         c(i, k) = b/ms;      end     end  c = transpose(c);  Để tìm ma trận nghịch đảo ta dùng chương trình ctminorinv.m:  clear all, clc;  a = [1  3  5;  3  4  9; 5  9  6];  b = minorinv(a)  §21. NGHỊCH ĐẢO BẰNG CÁCH PHÂN TÍCH MA TRẬN    Cho ma trận [A[, ta có thể phân tích nó thành ma trận tam giác phải [R]  và tam giác trái [L]:    [A] = [L][R]  Do định nghĩa ma trận nghịch đảo:    [L]‐1[L] = [L][L]‐1 = [E]  nên:    [R] = [L]‐1[L][R] = [L]‐1[A]  và:    [L] = [L][R][R]‐1 = [A][R]‐1  Do vậy:    [A] = [L][R] = ([A][R]‐1)([L]‐1[A]) = [A][R]‐1[L]‐1[A]   100   [A][R]‐1[L]‐1 = [E]   Kết quả là: [R]‐1[L]‐1 = [R]‐1   Với ma trận tam giác phải [R], các hàng khi nghịch đảo là l1,..,ln được tính theo  cách sau:    ‐ Cho [e1],.., [ei],..,[en] là các cột của ma trận đơn vị [E] cấp n    ‐  [ ]= nn n,n e l a   ‐  [ ]( )+ += − −Li i i ,i 1 i 1 i ,n n n,n 1l e a l a l a Ta xây dựng hàm luinverse() để thực hiên thuật toán tính định thức của ma  trận [R]:  function r = luinverse(a)  % Nghich dao ma tran tam giac phai  n = size(a, 1);  e = zeros(n, n);  c = zeros(n, 1);  l = e;  b = e;  for i = 1:n      c(i) = 1;      e(:, i) = c;      c(i) = 0;  end  l(:, n) = e(:, n)/a(n, n);  for i = n‐1:‐1:1      c = zeros(n, 1);      for k = i+1:n          c = c + a(i, k)*l(:, k);      end      l(:, i) = (e(:, i) ‐ c)/a(i, i);  end  r = lʹ;  101 Để nghich đảo ma trận tam giác ta dùng chương trình ctluinv.m:  clear all, clc  a = [1 2 3; 0 4 5;  0 0 2];  r = luinverse(a)  Để nghịch đảo ma trận tam giác trái ta dùng các quan hệ:    [L]‐1 = ([LT]‐1)T    vì [L]T sẽ  là một ma trận tam giác phải. Ta dùng chương trình ctluinverse.m  để nghịch đảo ma trận thông thường:   clear all, clc  a = [ 1  3  5; 3  4  9;  5  9  6];  [l, r] = doolittle(a);  b = luinverse(r)*(luinverse(lʹ))ʹ  §22. NGHỊCH ĐẢO MA TRẬN BẰNG THUẬT TOÁN GAUSS ‐ JORDAN  Cho ma trận [A]  và ma trận đơn vị [E] tương ứng. Dạng của ma trận [E]  cấp 4, là:  ⎟⎟ ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜ ⎝ ⎛ = 1000 0100 0010 0001 E    Như vậy, vấn đề  là ta cần tìm ma trận [A]‐1. Phương pháp  loại trừ để nhận  được ma trận nghịch đảo [A]‐1 được thực hiện qua n giai đoạn, mỗi một giai  đoạn gồm hai bước. Đối với giai đoạn thứ k:     ‐ chuẩn hoá phần tử akk bằng cách nhân hàng với nghịch đảo của nó    ‐ làm cho bằng không các phần tử phía trên và phía dưới đường chéo  cho đến cột thứ k. Khi k = n thì [A](k) sẽ trở thành ma trận đơn vị và [E]  trở thành [A]‐1  Ta xây dựng một hàm nghịch đảo invmat():  function x = invmat(a)  % Nghich dao ma tran a  102 %Cu phap: x = invmat(a)  k = size(a, 1);  n = k;  b = eye(n);  a = [a, b];  i = 1;  while i<=n      if a(i, i) ~= 0          c = a(i, i);          a(i, i:2*n) = a(i, i:2*n)/c;      end    for k = 1:n         if k~=i              c = a(k, i);              a(k, i:2*n) = a(k, i:2*n)‐ a(i, i:2*n)*c;          end      end    i = i+1;  end  x(:, 1:k) = a(:, (1+k):(2*k));  Để nghịch đảo ma trận     [ ] 2 1 1 A 1 2 1 1 1 2 ⎛ ⎞⎜ ⎟= ⎜ ⎟⎜ ⎟⎝ ⎠ ta dùng chương trình ctinvmat.m  clear all, clc  a = [ 2  1  1;  1  2  1;  1  1  2];  b = invmat(a)  §23. NGHỊCH ĐẢO BẰNG THUẬT TOÁN MOORE ‐ PENROSE    Cho ma trận [A] cấp m×n. Ma trận nghịch đảo tổng quát hoá Moore –  Pensore là ma trận n×m được Moore đưa ra năm 1920 và sau đó Pensore đưa  ra năm 1955. Ma trận Moore – Pensore [B] thoả mãn các điều kiện:  103   [A][B][A] = [A]    [B][A][B] = [B]  Từ các điều kiên này ta có:    [B] = [A]T([M][M])‐1  Để tính ma trận [B] ta dùng chương trình ctpseudoinv.m:  clear all, clc  a = [ 3 4 5 6 7; 3 6 7 8 9; 2 3 4 5 4];  b = aʹ*invmat(a*aʹ)  §24. GIÁ TRỊ RIÊNG VÀ VEC TƠ RIÊNG CỦA MA TRẬN  Cho ma trận [A] vuông, cấp n. Ta gọi số vô hướng λ là giá trị riêng của  ma trận [A] và [V] là vec tơ riêng phải của [A] tương ứng với λ nếu:    [A][V] = λ[V]                  (1)  [V] gọi là vec tơ riêng trái của [A] tương ứng với λ nếu:    [V]T[A] = λ[V]T  Ta có thể viết lại (1) dưới dạng:      [ ] [ ]( )[ ]− λ =A E V 0                   Từ đó ta có hệ n phương trình thuần nhất. Nghiệm tầm thường của hệ chính  là [X] = 0. Hệ sẽ có nghiệm không tầm thường nếu:    [ ] [ ]( )− λ =det A E 0                 (2)  Khai triển (2) ta nhận được phương trình đặc tính:    − +λ + λ + + λ + =Ln n 11 2 n n 1a a a a 0               Các nghiệm λi của phương trình đặc tính là các giá trị riêng của [A]. Nghiệm  [Vi] của:     [ ] [ ]( )[ ]− λ =iA E V 0   là các vec tơ riêng của [A]. Ta xây  dựng hàm nulld() để giải hệ phương trình  tuyến tính thuần nhất này:  function Z = nulld(A, small)  norma = sqrt(mean(mean(abs(A))));  tiny = norma*eps;  if nargin<2     small = (1 + norma)*sqrt(eps);  end  104 [m, n] = size(A);  if  m~= n     error(ʹMa tran phai vuong!ʹ)  end  p = find(abs(A)<tiny);  if ~isempty(p)     A(p) = 0;  end  [U,S,V] = svd(A,0);  S = diag(S);  s = S;  norma = max(s);  smax = max(s);  if smax == 0;      smax = 1;  end  s = s/smax;  snorm = s;  t = find(s>0);  if isempty(t);      Z = V;      return;  end  p = find(s<tiny);  if ~isempty(p)     s(p) = tiny;  end  p = find(s == 0);  if ~isempty(p)     s(p) = small*min(s(t));  end  logs = log10(s);  sdifflog = ‐diff(logs);  smax = max(sdifflog);  r = find(sdifflog == smax);  105 if min(size(r))>0     r = r(end);  else     r = 0;  end  Z = V(:,r+1:n);  if snorm == ones(n,1);        Z = zeros(n,0);   end  if max(S) <= small;       Z = V;   end  §25. BIẾN ĐỔI ĐỒNG DẠNG  Ta khảo sát bài toán về giá trị riêng của ma trận:    [A][X] = λ[X]                  (1)  với [A] là ma trận đối xứng. Bây giờ ta áp dụng phép biến đổi:    [X] = [P][X]*                   (2)  với [P] là ma trận không suy biến.  Thay (2) vào (1) và nhân hai vế với [P]‐1 ta có:    [ ] [ ][ ][ ] [ ] [ ][ ]− − ∗= λ1 1P A P X P P X   hay:  [ ] [ ] [ ]∗ ∗ ∗= λA X X                   (3)  Trong đó   [ ] [ ][ ] [ ] [ ][ ]∗ −= 1A P X P A P   Do λ không đổi khi thực hiện phép biến đổi nên giá trị riêng của ma trận [A]  cũng chính là giá trị riêng của ma trận [A]*. Các ma trận có cùng giá trị riêng  được gọi là ma trận đồng dạng và phép biến đổi giữa chúng gọi là phép biến  đổi đồng dạng.   Phép biến đổi đồng dạng thường dùng để đưa bài toán tìm giá trị riêng  về dạng dễ giải hơn. Giả sử ta có cách nào đó để tìm ma trận [P] mà nó đường  chéo hoá ma trận [A], lúc đó (3) có dạng:  106 ∗ ∗ ∗ ∗ ∗ ∗ ⎡ ⎤ ⎡ ⎤− λ ⎡ ⎤⎢ ⎥ ⎢ ⎥ ⎢ ⎥− λ⎢ ⎥ ⎢ ⎥ ⎢ ⎥=⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥− λ ⎣ ⎦⎣ ⎦ ⎣ ⎦ L L MM M L M L 11 1 22 1 nn 1 0A 0 0 x 00 A 0 x 0 00 0 A x Nghiệm của phương trình trên là:    ∗ ∗ ∗λ = λ = λ =L1 11 2 22 n nnA A A             (4)  [ ] [ ] [ ]∗ ∗ ∗ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥= = =⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎣ ⎦ ⎣ ⎦ LM M M1 2 n 1 0 0 0 1 0 X X X 0 0 1 hay:  [ ] [ ]∗ ∗ ∗ ∗⎡ ⎤= =⎣ ⎦L1 2 nX X X X E   Theo (2), vec tơ riêng của [A] là:    [X] = [P][X]* = [P][E] = [P]              (5)  Như vậy ma trận biến đổi [P] là ma trận các vec tơ riêng của [A] và các  giá trị riêng của [A] là các số hạng trên đường chéo của [A]*.  §26. TÌM GIÁ TRỊ RIÊNG BẰNG CÁC PHƯƠNG PHÁP LUỸ THỪA  1. Phương pháp luỹ thừa vô hướng: Giả sử ma trận vuông [A] cấp n có n giá  trị riêng phân biệt:    λ > λ ≥ λ ≥ ≥ λL1 2 3 n                (1)  và các vec tơ riêng tương ứng [V1], [V2],...,[Vn]  Ta chọn một vec tơ [X] bất kì có các thành phần khác zero. Khi đó [X] sẽ được  biểu diễn bằng:    [ ] [ ] [ ] [ ]= α + α + + αL1 1 2 2 n nX V V V             (2)  Do [ ][ ] [ ]= λn n nA V V  nên:  [ ][ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ][ ] [ ][ ] [ ][ ] [ ] [ ] [ ] [ ] [ ] [ ] = α + α + + α = α + α + + α = α λ + α λ + + α λ ⎛ ⎞λ λ= λ α + α + + α⎜ ⎟λ λ⎝ ⎠ L L L L 1 1 2 2 n n 1 1 2 2 n n 1 1 1 2 2 2 n n n 2 n 1 1 1 2 2 n 2 1 1 A X A V A V A V A V A V A V V V V V V V Lặp lại lần thứ k ta có:  107   [ ] [ ] [ ] [ ] [ ] [ ]⎧ ⎫⎛ ⎞ ⎛ ⎞λ λ⎪ ⎪= = λ α + α + + α⎨ ⎬⎜ ⎟ ⎜ ⎟λ λ⎝ ⎠ ⎝ ⎠⎪ ⎪⎩ ⎭L k k k k 2 n k 0 1 1 1 2 2 n 2 1 1 X A X V V V   (3)  Khi k → ∞, do (1) nên (3) hội tụ về [V1] khi α1 ≠ 0. Vậy là giá trị riêng λ1, có trị  tuyệt đối lớn nhất, và vec tơ riêng tương ứng [V1], có thể tính bằng cách cho  trước một vec tơ [X0] có các thành phần khác zero theo hướng [V1] và lặp theo  thủ tục sau:    [ ] [ ] [ ][ ] [ ]+ ∞= →λ k k 1 1 1 k X X A V X               (4)  Trong đó:    [ ] [ ]{ }∞= nX max X   Tốc độ hội tụ phụ thuộc vào độ lớn của  λ λ2 1 .  Thuật toán tóm lại gồm các bước:    ‐ cho vec tơ ban đầu [X0]    ‐ tính [X1] = [A][X0]  ‐ tính  [ ]X   ‐ cho [X] = [ ] [ ] [ ]=X X X  và lặp lại các bước 2 ‐ 4 cho đến khi hội tụ  Ta xây dựng hàm eigpower() để thực hiện thuật toán trên:  function [lambda, x, iter] = eigpower(A, tol, maxiter, x0)  % Tim gia tri rieng l