Giáo trình tính toán khoa học - Chương 3

Chia sẻ: Nguyen Thanhvan | Ngày: | Loại File: PDF | Số trang:31

Thêm vào BST

Báo xấu

310
lượt xem 91
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Để các bạn làm quen với các công cụ của Matlab trong đại số tuyến tính, trước hết chúng tôi cần nhắc lại một số khái niệm về ma trận và các phép toán trên ma trận. 3.1.1 Chuyển vị ma trận Cho ma trận A=(aij)mxn. Ma trận chuyển vị của A là ma trận A' =(a'ij)nxm sao cho a'ij=aji. Nếu A’=A thì A được gọi là một ma trận đối xứng.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Giáo trình tính toán khoa học - Chương 3

Chương 3 GIẢI TÍCH MA TRẬN VÀ ĐẠI SỐ TUYẾN TÍNH 3.1 GIẢI TÍCH MA TRẬN Để các bạn làm quen với các công cụ của Matlab trong đại số tuyến tính, trước hết chúng tôi cần nhắc lại một số khái niệm về ma trận và các phép toán trên ma trận. 3.1.1 Chuyển vị ma trận Cho ma trận A=(aij)mxn. Ma trận chuyển vị của A là ma trận A' =(a'ij)nxm sao cho a 'ij=aji. Nếu A’=A thì A được gọi là một ma trận đối xứng. Thí dụ 1. 1 5 9 1 2 3 4  2 6 10  Nếu A   5 6 7 8  thì A '   . 7 11  3  9 10 11 12      4 8 12  3.1.2 Định thức của ma trận vuông Cho ma trận vuông A cấp n . Ta gọi ma trận con của ma trận A tương ứng với phần tử a ij là ma trận vuông cấp n-1 suy từ A b ằng cách bỏ đi các phần tử hàng i và cột j . Thí dụ 2. Sau đây là ma trận A và các ma trận con M23 và M12 tương ứng với các phần tử a23 và a12 của nó:  1 2 3  1 2  4 6 A   4 5 6  , M 23    , M12   7 9  . 7 8  7 8 9     Định thức của ma trận A vuông cấp n , gọi là định thức cấp n, được định nghĩa theo phương pháp qui n ạp như sau: - Nếu A là ma trận cấp 1 hay A=(a11), thì det(A)=a11; 62
a11 a12  - Nếu A là ma trận cấp 2 hay A   thì  a21 a22    d et(A) = a 11 a22 - a12 a 21=a 11det(M11) - a12 d et(M12 ); - Tổng quát: Nếu A là ma trận cấp n≥2 thì det( A)  a11det(M11)-a12 det(M12 )+...+ (-1)1+n a1n d et(M1n ), n  (1) j 1 a1 j det(M1 j ) . hay (3.1) det( A)  j 1 Để tính định thức của ma trận cấp n theo công thức qui nạp, ta phải tính n định thức cấp n -1. Do đó số phép tính để tính định thức của ma trận cấp n sẽ tăng theo tốc độ của n!. Trong thực tế, khi cần giải các b ài toán về đại số tuyến tính cỡ lớn, người ta cố gắng tránh các phương pháp có sử dụng định thức. 3.1.3 Ma trận trận nghịch đảo Giả sử A là ma trận vuông cấp n. Nếu tồn tại ma trận B vuông cấp n sao cho AB=BA=E (E là m a trận đơn vị cấp n) th ì A được gọi là ma trận khả nghịch. Khi đó B được gọi là ma trận nghịch đảo của A và được kí hiệu là B =A-1 . Chỉ ma trận có định thức khác không (hay ma trận không suy biến) mới khả nghịch đảo và ma trận nghịch đảo của A có thể tính đ ược bằng công thức:  c11 c21 ... nn1  1  c12 c22 ... cn2  1 A1  CT  . (3.2)  det( A)  ... ... ... ...  det( A) c ... cnn   1n c2 n  trong đó C=(cij)nxn, với cij=(-1)i+jdet(Mij) gọi là phần phụ đại số của phần tử a ij của ma trận A. Chú ý rằng với khái niệm về phần phụ đại số ta có: n n  aij cij   aij cij . det( A)  j 1 i 1 Với công thức trên, đ ể tính ma trận nghịch đảo của ma trận cấp n thì cần phải tính n 2 định thức cấp n-1 . Rõ ràng đ ây là phương pháp rất kém hiệu quả. Vì vậy ta cần phải tìm kiếm các phương pháp tính khác, hiệu quả hơn. Sau đây là một số các hàm ma trận và vector trong Matlab: 63
Bảng 3-1 Một số hàm ma trận và vector trong Matlab Hàm Ý nghĩa Tính m a trận nghịch đảo của ma trận vuông A. inv(A) Tính định thức của ma trận vuông A. det(A) A' hoặc A.' Tạo ma trận chuyển vị của ma trận A. Hàm vết của ma trận hay tổng các phần tử trên đường chéo trace(A) của ma trận A. Tính hạng của ma trận A. rank(A) [m,k]=min(x) Tính giá trị nhỏ nhất m trong các to ạ độ của vector x và vị trí k đạt min. Nếu x là ma trận thì kết quả là vector hàng gồm giá trị min của các cột. Tính giá trị lớn nhất M trong các toạ độ của vector x và vị trí [M,k]=max( k đạt max. Nếu x là ma trận thì kết quả là vector hàng gồm giá x) trị max của các cột. Tính giá trị trung b ình của các phần tử của vector x. Nếu x là mean(x) ma trận thì kết quả là vector hàng gồm giá trị trung bình của các cột. [y,k]=sort(x) Sắp xếp lại các phần tử của x theo thứ tự tăng dần, kết quả trả cho vector y. Vector k là vector số thứ tự cũ trong x của các ph ần tử trong y. Nếu x là ma trận thì các cột của x được sắp xếp tăng dần. Tính tổng các phần tử của vector x. Nếu x là một ma trận thì sum(x) kết quả là 1 vector hàng, m à mỗi phần tử của vector là tổng các ph ần tử của một cột tương ứng. Cộng dồn các phần tử của vector x. cumsum(x) Tính tích các phần tử của vector x. Nếu x là một ma trận thì prod(x) kết quả là 1 vector hàng, m à mỗi phần tử của vector là tích các ph ần tử của một cột tương ứng. Nhân dồn các phần tử của vector x. cumprod 64
Thí dụ 3. >> A=[ 1 2 3 151 3 2 1]; >> det(A) ans = -32 >> inv(A) ans = -3/32 -1/8 13/32 -1/16 1/4 -1/16 13/32 -1/8 -3/32 >> max(A) ans = 3 5 3 >> min(A) ans = 1 2 1 >> sum(A) ans = 5 9 5 >> cumsum(A) ans = 1 2 3 2 7 4 5 9 5 >> x= [ 0 0 1 -2 3 4 5 -6 -7 -8]; >> [M,k]=max(x) M= 65
5 k= 7 Chú ý k là vị trí đầu tiên của phần tử đạt max trong x.Tương tự ta có: >> [m,p]=min(x) m= -8 p= 9 >> mean(x) ans = -10/9 >> sort(x) ans = -8 -7 -6 -2 0 0 1 3 4 5 >> [y,k]=sort(x) y= -8 -7 -6 -2 0 0 1 3 4 5 k= 10 9 8 4 1 2 3 5 6 7 >> D= sort(A) D= 1 2 1 1 2 1 3 5 3 >> prod(A) ans = 3 20 3 66
>> cumprod(A) ans = 1 2 3 1 10 3 3 20 3 >>rank (D) ans = 2 3.1.4 Tích vô hướng Trong không gian vector V, tích vô hướng của 2 vector x và y, được ký hiệu là , là một số thực ho àn toàn xác định đối với mọi cặp vector x và y, thoả mãn các tính chất: - = < y, x >; - = < x, z > + ; - = = k , với mọi kR; -  0 ; = 0  x=0. Trong mỗi không gian vector có th ể định nghĩa nhiều tích vô hướng khác nhau. Tuy nhiên, tích vô hướng Euclide là loại tích vô h ướng được sử dụng nhiều nhất trong không gian Rn. Tích vô hướng Euclide của hai vector x và y được xác định như sau: n  xi yi . (3.3)  x, y  x1 y 1 x2 y2  ...  xn yn  i 1 3.1.5 Chuẩn của vector Cho một hàm vector hàm f(x) xác định trên không gian vector V, ký hiệu là x . Nếu nó thoả mãn các tính chất sau đây với xV: x  0, x  0  x  0 , 1) x   . x , R, 2) x  y  x  y (gọi là bất đẳng thức tam giác), 3) thì x đ ược gọi là chuẩn của vector x. 67
Từ đ ịnh nghĩa về chuẩn vector, với mỗi số p>0 có một chuẩn loại p tương ứng trong Rn như sau : 1/ p p n     xi  x với 0 < p  +. p  i1  Trong đó có 3 loại chuẩn vector thường được sử dụng nhiều nh ất là:  max  xi  ; - x  i - x 1 = x1  x2  ...  xn ; 1   2 2 2 - x 2  x1  x2  ...  xn 2 . Từ đó có 3 loại chuẩn của ma trận tương thích với chuẩn vector. Cho A là một ma trận vuông A =( a ij )nn, khi đó:     - A  = max  aij  ; i j        - A 1 = max   aij  ; j i      j  , với  là các trị riêng của ma trận đối xứng ATA. - A 2 = max j j j còn được gọi là các trị kì dị (single value) của ma trận A. Các chuẩn của ma trận cũng thỏa m ãn các tính chất 1)-3) của chuẩn vector. Ngoài ra, giữa chuẩn loại p=0,1,2 của vector và chuẩn loại tươn g ứng của ma trận tho ả m ãn tính chất sau: Ax p  A p x p . Bất đẳng thức trên gọi là tính tương thích của chuẩn ma trận đối với chuẩn vector. Chu ẩn loại 2 cúa vector trong không gian Rn: 2 2 2 x 2  x1  x2  ...  xn còn được gọi là chuẩn Euclide. Thí dụ 4. Để tính các tích vô hướng và các chuẩn m a trận và chu ẩn vector trong Matlab có th ể làm như sau: 68
>> x = [1 2 3 4]; >> y = [ 5 6 7 8]; %% Tích vô hướng Euclide >> tvh = x*y' tvh= 70 >> chuanVC=max(abs(x)) chuanVC = 4 >> chuan1=sum(abs(x)) chuan1 = 10 >> chuan2 = sqrt(x*x') chuan2 = 5.4772 %% Lỗi : sai kích thước >> chuan =sqrt(x' *x); >> A=[ 1 2 3; 4 5 6; 7 8 9]; >>ChuanVC =max(sum(abs(A'))) ChuanVC = 24 >> Chuan1 =max(sum(abs(A))) Chuan1 = 18  Hàm NORM Cú pháp: norm(V,p) Giải thích. Hàm NORM tính chuẩn loại p của ma trận và vector V. - Nếu V là một vector p > 0 hoặc  inf thì : Chu ẩn loại p; norm (V,p) = sum(abs(V).^p)^(1/p) : Chu ẩn Euclide; norm (V) = norm(V,2) 69
: Chu ẩn loại vô cùng; norm (V,inf) = max(abs(V)) norm (V,-inf) = min(abs(V)). - Nếu V là một ma trận th ì p ch ỉ có thể 1, 2, inf hoặc 'fro' và : Trị kì dị lớn nhất của V; norm (V) = max(svd(V)) norm (V,2) = norm (V); norm (V,1) = max(sum(abs(V))); norm (V,inf) = max(sum(abs(V'))); norm (V,'fro')= sqrt(sum(diag(V'*V))) : chu ẩn Frobenius. 3.1.6 Phân loại ma trận Khi nghiên cứu các bài toán trong đại số tuyến tính, người ta thường chia ma trận th ành 2 loại: - Ma trận lưu trữ đ ược: Các ma trận mà các phần tử của chúng có thể lưu trữ và xử lí đ ược trong bộ nhớ của MTĐT. Kích thước của các ma trận này thường không lớn lắm n = 0(103). - Ma trận thưa: Ma trận có kích thước rất lớn, nhưng phần lớn các phần tử của chúng đều bằng 0 hoặc ma trận có phần tử không cần lưu trữ m à có th ể tính được bằng qui tắc nào đó. Thí dụ 5. Một ma trận thưa có dạng 3 đường chéo: 4 1 0 0 ... 0 0 0 1 0 0 0 4 1 0 ... 0 0 0 0 1 4 1 ...   0 0 0. A0 0 1 4 ...  ... ... ... ... ... 0 ... ...  0 1 4 1 0 0 0 ... 0 0 1 4 0 0 0 ...   3.1.7 Số điều kiện của ma trận . là một chuẩn nào đó của vector và ma trận và A là một ma trận Giả sử vuông. Nếu đặt : Ax Ax M = sup và m = inf x x x 0 x 0 70
1 thì dễ dàng chứng minh rằng M = A và nếu A không suy biến thì m = A -1 . M  A . A -1 .được gọi là số điều kiện của ma trận A và được ký Khi đó, tỷ số m hiệu là cond(A). Giả sử x là nghiệm đúng của hệ phương trình tuyến tính Ax = b và x+x là nghiệm của hệ phương trình xấp xỉ A( x  x)  b  b . Khi đó b  Ax  M x và b  Ax  m x . Từ đó suy ra: x M b b . (3.4)   cond ( A) x mb b Biểu thức trên cho thấy sai số tương đối của nghiệm có thể lớn bằng sai số tương đối của vế phải nhân với cond(A). Ma trận A được gọi là ma trận có số điều kiện xấu (ill-conditioned -matrix ) n ếu cond(A)>>1. Vì vậy việc giải gần đúng một hệ phương trình với ma trận hệ số A có số điều kiện xấu rất kém ổn định; Vế phải b chỉ cần thay đổi nhỏ cũng có thể gây sai số khá lớn cho lời giải của hệ phương trình. Số điều kiện của ma trận có các tính chất: i) cond(A)  1; ii) Nếu A là ma trận trực giao (tức là A'=A-1) thì con(A)=1; iii) Với mọi c  0 R : ta đều có cond(cA) = cond(A); max di n iv) Nếu D = diag d i 1 thì cond(D) = . min di 100 Thí dụ 6. Giả sử A=diag 0,11 , khi đó det(A)=10-100 nên A là ma trận gần suy biến. Tuy nhiên do A=0,1E nên A là ma trận có số điều kiện rất tốt và hệ phương trình Ax=1 dễ dàng giải đ ược một cách chính xác là: xi =10 , i  1,100 .  Hàm COND Cú pháp: cond(A,p) Giải thích. Hàm COND tính số điều kiện của ma trận A theo chu ẩn loại p. cond (A,p)= norm(A,p)*norm(inv(A),p), nếu p=1,2, inf hoặc 'fro'; cond (A) = cond(A,2). 71
Thí dụ 7. >> A=hilb(10); >> C = cond(A) C= 1.6025e+013 Ta thấy ma trận Hilbert là ma trận có số điều kiện rất xấu. 3.2 MỘT SỐ MA TRẬN ĐẶC BIỆT TRONG MATLAB Bảng 3-2 Danh sách một số ma trận đặc biệt có trong Matlab Ma trận Ý nghĩa Ma trận ma phương cấp n magic(n) Ma trận Pascal cấp n pascal(n) Ma trận Hadamard cấp n , với n>2. Nó chỉ tồn tại khi n hadamard(n) chia hết cho 4. Ma trận Hilber cấp n hilb(n) Ma trận nghịch đảo của ma trận Hilbert cấp n invhilb(n) Ma trận tích tensor Kronecker của ma trận x và y. kron(x,y) Ma trận Vandermonde của vector x vander(x) Ma trận ma phương là ma trận vuông cấp n sinh bởi các số 1, 2, 3,.., n2 có các tổng của cột, các tổng h àng và các tổng đường chéo đều bằng nhau. Thí dụ 8. Ma trận magic cấp 3: 8 1 6 magic(3) =  3 7. 5 4 2 9   Ma trận Pascal là ma trận vuông cấp n chứa số của tam giác Pascal. Thí dụ 9. Ma trận Pascal cấp 4: 72
1 1 1 1 1 2 3 4 . pascal (4)=  3 6 10  1   1 4 10 20  Ma trận Hilber là ma trận vuông cấp n, có phần tử a ij = 1/(i+j-1). Thí dụ 10. %% Ma trận Hilbert cấp 4 >> A= hilb(4) A= 1 .0000 0.5000 0.3333 0.2500 0 .5000 0.3333 0.2500 0.2000 0 .3333 0.2500 0.2000 0.1667 0 .2500 0.2000 0.1667 0.1429 >> B=hilb(10); >> C = cond(B) C= 1.6025e+013 Kết quả trên cho th ấy ma trận Hilbert là ma trận có số điều kiện rất xấu. Đây là thí dụ điển hình về ma trận có số điều kiện xấu. Do đó ma trận Hilbert thường được dùng đ ể thử nghiệm về tính ổn định nghiệm của một phương pháp giải hệ phương trình tuyến tính. Ma trận Hadamard là ma trận vuông H cấp n gồm các phần tử 1 và -1 sao cho H'*H= n*eye(n). Nó chỉ tồn tại khi n chia hết cho 4. Matlab chỉ tính được ma trận Hadamard khi n, n/12 hay n/20 là lũy thừa của 2. Thí dụ 10. >> H=hadamard(4) H= 1 1 1 1 1 -1 1 -1 1 1 -1 -1 1 -1 -1 1 Ma trận Vandermonde của vector x=(x1,x2,.., xn) là ma trận: 73
 x1n1 ... x12 x1 1   n1  2 vander(x) =  x 2 ... x2 x2 1  .  ... ... ... ...  x n1 ... xn 1  2 xn n  Ma trận Vandermonde được sử dụng trong nội suy bằng đa thức. Tích tensor Kronecker của 2 ma trận Amxn và B là ma trận có dạng:  a11B a12 B ... a1n B    a B a22 B ... a2 n B  kron(A,B) =  21 .  ... ...  ... ...    am1B am 2 B ... amn B  Bảng 3-3 Các hàm trích phần tử từ một ma trận. Hàm Ý nghĩa diag(A,p) Nếu A là 1 ma trận thì hàm sẽ tạo ra 1 vector cột từ đường chéo thứ p của ma trận A. Nếu A là 1vector thì hàm sẽ tạo ra 1 ma trận vuông có đường chéo thứ p là vector A, còn các phần tử khác đều bằng 0. Mặc định của p là 0 (đường chéo chính). Hàm sẽ tạo ra một ma trận cùng cỡ với ma trận A, có các phần tử từ triu(A,p) đường chéo thứ p trở lên được lấy từ A, còn các ph ần tử khác (phía dưới đường chéo thứ p) đều bằng 0. Mặc định của p là 0. Hàm sẽ tạo ra một ma trận cùng cỡ với ma trận A, có các phần tử từ tril(A,p) đường chéo thứ p trở xuống đ ược lấy từ A, còn các phần tử khác (phía trên đường chéo thứ p) đều bằng 0. Mặc định của p là 0. Thí dụ 11. >> A =[ 1 2 34 56 7 8 9 10 11 12]; >> B = diag(A) B= 74
1 6 11 >> C = diag(A, 2) C= 3 8 >> D = diag(A,-3) D= Empty matrix: 0 -by-1 >> E = tril(A) E= 1 0 0 0 5 6 0 0 9 10 11 0 >> F=diag(B) F= 1 0 0 0 6 0 0 0 11 >> G = triu(A) G= 1 2 3 4 0 6 7 8 0 0 11 12 >> H = triu(A,2) H= 00 3 4 75
00 0 8 00 0 0 >> C = tril(A,-2) C= 0 0 0 0 0 0 0 0 9 0 0 0 3.3 CÁC PHƯƠNG PHÁP GIẢI HỆ PHƯƠNG TRÌNH TUYẾN TÍNH Xét h ệ phương trình đ ại số tuyến tính có dạng Ax=b, trong đó A là m ột ma trận vuông cấp n và b là một vector cột n chiều. Để giải hệ ph ương trình Ax=b ta có các phương pháp giải như sau: 3.3.1 Phương pháp Cramer Nếu A là ma trận không suy biến th ì hệ phương trình có duy nh ất nghiệm xác đ ịnh bởi biểu thức: det( A j ) (3.5) xj  ; j  1, n det( A) trong đó Aj là ma trận vuông cấp n suy từ A bằng cách thay cột thứ j của A bởi vector vế phải b. Xét về hiệu quả tính toán thì phương pháp này rất tồi. Người ta tính được số phép tính số học (+, -,  và / ) của phương pháp này là: NC(n)=(n+1)!n . Bạn h ãy tư ởng tượng, n ếu phải giải hệ phương trình 15 ẩn (n=15) theo phương pháp Cramer trên MTĐT có tốc độ tính toán là 1 tỷ phép tính/giây thì phải mất 3,6 ngày tính toán liên tục; Còn nếu số ẩn là n =20 thì thời gian tính toán là 32 nghìn năm. Những con số thật khủng khiếp. Tuy nhiên công thức Cramer có ý ngh ĩa khi ta nghiên cứu một số tính chất của lời giải. 76
3.3.2 Phương pháp ma trận nghịch đảo Đây là một hình thức th ể hiện khác của phương pháp Cramer. Nếu A là ma trận không suy biến thì từ biểu thức Ax=b suy ra hệ phương trình có duy nh ất nghiệm xác định bởi biểu thức x=A-1b. Vì vậy bạn dễ dàng có thể sử dụng Matlab để tính nghịch đảo của ma trận A và giải hệ phương trình như sau: >> x= inv(A)*b Tho ạt nhìn công thức thì có vẻ đơn giản và rõ ràng, nhưng phương pháp trên lại rất kém hiệu quả khi tính toán. Vì vậy việc giải hệ phương trình tuyến tính trên MTĐT rất hiếm khi n gười ta sử dụng ma trận nghịch đảo A-1. Tuy nhiên việc tính ma trận nghịch đảo rõ ràng là vẫn rất cần thiết bởi vì biểu diễn trên là công thức “chính xác” để giải b ài toán đại số tuyến tính. Nhiều khi chúng ta cần quan tâm vector nghiệm x và các tính chất của nó. 3.3.3 Phương pháp Gauss Xét h ệ phương trình Ax =b trong trường hợp A là ma trận tam giác trên và các hệ số trên đường chéo chính đều khác 0:  a11 a12 ... a1,n 1 a1n   x   b   1   1    0 a22 ... a2, n 1 a 2 n  x2   b2    ...   ...    ...  ... ... ... ...     an 1,n 1 an 1, n   xn 1   bn 1  0 0 ...     0 ann   xn   bn  0 ... 0   Khi đó hệ phương trình Ax=b có thể giải một cách đơn giản bằng phương pháp thế ngư ợc từ dưới lên theo công thức sau đây mà không phải tính một định thức n ào: b xn  n , ann n  bk  akj x j bk - ak , k 1xk 1 - ak ,k  2 xk  2 - ...- ak n xn j  k 1 , xk   akk akk với k  n, n  1,...,1 . Tương tự, nếu A là ma trận tam giác dưới với các phần tử trên đường chéo chính khác 0, thì có thể giải hệ phương trình bằng phương pháp thế xuôi từ trên 77
xuống như sau: b x1  1 , a11 k 1  akj x j bk  bk - ak1x1 - ak 2 x2 - ...- ak , k 1xk 1 j 1 , với k  2, n xk   akk akk  Phương pháp Gauss Nếu A không phải một ma trận tam giác, có thể dùng các phép biến đổi sơ cấp về h àng của ma trận để biến đổi dần ma trận mở rộng A =[A,b] của ma trận A về dạng ma trận tam giác trên, rồi dùng phương pháp thế n gược nói trên đ ể tính vector x. Công việc đó được gọi là quá trình xuôi của phương pháp Gauss: - Trước hết đặt A(0)=[A,b]. - Quá trình xuôi gồm n -1 bước. Tại bư ớc thứ k=1, n  1 thực hiện : + Đặt A(k)=A(k-1) ;  k 1  k 1  k   a  k 1  aik .akj , j  k , n  1 . + Tại hàng i >k tính (3.6) aij ij  k 1 a kk k1 Nếu trong quá trình tính toán các phần trử trụ akk đều khác không thì khi (n-1) kết thúc quá trình xuôi ta được ma trận A có dạng tam giác trên.  Phương pháp Gauss -Jordan (Còn gọi là phương pháp Gauss cải tiến hay p hương pháp Chọn trụ tối đại). Trong quá trình tính toán của ph ương pháp Gauss, nếu gặp phần tử trụ  k 1  0 thì không thể tiến hành bước tiếp theo được; Hoặc nếu phần tử trụ có akk giá trị quá nhỏ th ì kết quả của phép chia cũng sẽ kém chính xác vì vậy cần phải cải tiến phương pháp Gauss như sau : Đặt A(0)= [A,b ]. - Quá trình xuôi gồm n bước. Tại b ước thứ k=1, n thực hiện : -  k 1 + Tìm r  arg max  aik  : i  k ,n  ;    78
+ Nếu ark=0 thì ma trận A suy biến, nên phải ngừng giải; Ngược lại, nếu r  k thì tiến h ành đổi chỗ 2 h àng r và k của ma trận A(k-1); + Đặt A(k)=A(k-1) ;  k 1  k   akj + Tại h àng k tính akj , j  k, n  1;  k 1 a kk  k 1  k 1  k   a  k 1  aik .akj , j  k , n  1, + Tại h àng i >k tính aij ij  k 1 a kk k   k 1  a k 1 .a k  , j  k , n  1 . hay aij  aij ik kj Kết quả của thủ tục là thu được ma trận tam giác trên A(n) có các phần tử trên đường chéo chính bằng 1. Có thể tính đ ược số phép tính số học cần sử dụng cho phương pháp Gauss (kể cả 2 quá trình xuôi và ngược) là: 1 (4n3  9 n 2  7 n) . NG  n   6 3.3.4 Phương pháp phân tích QR và phương pháp phân tích LU Để giải hệ phương trình tuyến tính trên MTĐT người ta thường phân tích ma trận A thành tích của 2 ma trận, mỗi ma trận đó dễ nghịch đảo hơn ma trận A. Do đó sẽ giảm nhẹ khối lượng tính toán. Có 2 phương pháp phân tích ma trận:  Phương pháp phân tích QR Phân tích ma trận A thành tích của ma trận trực giao Q (ma trận unitar) và ma trận tam giác trên R. Ta có thể viết lại hệ phương trình gốc nh ư sau: Ax = QRx = b Việc tính nghịch đ ảo của ma trận Q là tầm thường; Và khi đó hệ phương trình trở th ành: Rx = QTb Hệ phương trình mới có dạng tam giác trên nên có thể dễ d àng giải bằng quá trình th ế ngược. Như vậy việc giải hệ không cần p hải tính nghịch đảo của ma trận A. Để xác định các ma trận Q và R có thể sử dụng hàm QR của Matlab có cú pháp như sau: [Q,R,P] =qr(A) 79
Giải thích: - A là một ma trận vuông cần phân tích; - Q là ma trận trực giao và R là ma trận tam giác trên; - P là ma trận giao hoán thỏa m ãn Q*R = P*A. Sở dĩ cần phải sử dụng thêm ma trận giao hoán vì trong quá trình tính toán, giống như phương pháp Gauss, khi gặp trường hợp phải chia cho 0 thì cần phải hoán vị các h àng của ma trận.  Phương pháp phân tích LU Hãy xem xét lại phương pháp khử Gauss. Tại b ước thứ k=1, n  1 của phương pháp Gauss, ta giữ lại k hàng đầu của ma trận A(k-1) và biến đổi phần tử aikk 1 (i>k) thành 0. Mỗi bước như vậy tương đương với việc nhân bên trái ma  trận A(k-1) với một ma trận sơ cấp có dạng tam giác dưới Ek. Quá trình đó đưa dần dần m a trận A về dạng ma trận tam giác trên U có các ph ần tử trên đường chéo bằng 1. Nói cách khác EnEn-1...E1A= MA=U hay A =M-1U=LU. Phép phân tích này th ực hiện được nếu nh ư tất cả các phần tử trụ  k 1  0, k  1,n tức là các đ ịnh thức con chính của ma trận A đều phải khác akk không. Tuy nhiên trong khi tính toán ta có thể gặp trường hợp akk   0 , giống  k 1 như phương pháp Gauss cải tiến, ta cần đổi chỗ hàng k với một hàng r nào đó ở dưới. Vì vậy cũng cần lưu trữ các phép đổi chỗ bằng một ma trận giao hoán để sử dụng cho việc biến đổi tương ứng vector vế phải của hệ phương trình . Để tìm các ma trận L và U ta có th ể dùng thủ tục Crout gồm n-1 bước. Vì L là ma trận tam giác dưới có các phần tử đường chéo bằng 1 và U là ma trận tam giác trên nên trong khi tính toán ta có thể lưu trữ chúng chung vào ma trận kí hiệu là B có dạng:  u11 u12 ... u1n    l u22 ... u2 n  B   21 .  ... ... ...  ...    ln1 ln 2 ... unn  Quá trình tính toán đ ược thực hiện theo thủ tục gồm n bước sau đây. 80
 Thủ tục Crout Bước 1 . Đặt D1 = A. Bước k=2,3...n. Tại bước k-1 ma trận A đã được biến đổi th ành:  u11 u12 ... u1k 1 ... u1n     l21 u22 ... u2 k 1 ... u2 n   ... ...  ... ... ... ...  ,  lk 1,1 lk 1,2 ... u k -1, k -1 ... uk 1, n   ...  ... ...    ln1  ln 2 ... ln ,k -1   trong đó Dk là ma trận vuông cấp n-k+1 (tức là có k-1 hàng và k-1 cột của B T T  u kk uk  u uk   thành  kk đã được tính toán xong). Tiếp tục biến đổi Dk=    sk H k 1   lk Dk 1  theo công thức: 1 T sk và Dk 1  H k 1  lk uk . lk  ukk Bằng phép phân tích trên ta đã biểu diễn A thành tích của hai ma trận tam giác dưới L và tam giác trên U. Vì vậy chúng ta có: Ax = LUx = b Đặt Ux = y, rồi giải hệ phương trình Ly = b. Vì L là tam giác dưới nên d ễ dàng giải bằng quá trình thế xuôi để tính y. Sau đó ta lại giải hệ phương trình Ux=y b ằng quá trình th ế ngược để tìm nghiệm x. Như vậy chúng ta cũng không cần phải tính ma trận ngh ịch đảo n ào cả. Để xác định các ma trận L và U có thể sử dụng lệnh của Matlab theo cú pháp như sau: [L,U,P] =lu(A) Giải thích: - A là một ma trận vuông cần phân tích ; - L là ma trận tam giác dưới và U là ma trận tam giác trên; - P là ma trận giao hoán thỏa m ãn L*U = P*A 81