Bài giảng Phương pháp số - Bài 3: Ma trận và Hệ phương trình tuyến tính - Nguyễn Thị Vinh

BÀI 3 MA TRẬN VÀ HỆ PHƯƠNG TRÌNH TUYẾN TÍNH HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (1) 1. HỆ PHƯƠNG TRÌNH TUYẾN TÍNH gồm m phương trình n ẩn là một hệ có dạng           mnmm22m11m 2nn2222121 1nn1212111 bxa ...xaxa ............................................ bxa ... xaxa bxa ...xaxa thì hệ trên còn có thể viết ở dạng vectơ cột x1v1 +x2v2 ++ xnvn = b hay dạng phương trình ma trận Ax = b PHƯƠNG PHÁP SỐ-Bài 3 2 và x x x x, b b b bn),1,2,...,(j, a . a a v n 2 1 m 2 1 mj 2j 1j j                                           , a...aa ............ a...aa a...aa A mnm2m1 2n2221 1n1211              Nếu đặt HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (2) 2. HỆ DẠNG TAM GIÁC TRÊN VÀ CÁCH GIẢI Hệ dạng tam giác trên là hệ có dạng a11x1 + a12x2 ++ a1nxn = b1 trong đó a11, a22,,ann ≠ 0 a22x2 ++ a2nxn = b2 Cách giải: giải ngược từ dưới lên annxn = bn void heTGiac(vector > a, vector &x) { unsigned n = a.size(); vector y(n, 0); // y co n phan tu 0 y[n-1] = a[n-1][n]/a[n-1][n-1]; for (int i = n-2; i >= 0; i--) { double tong = 0.; for(unsigned j = i+1; j <= n-1; j++) tong = tong + a[i][j] * y[j]; y[i] = (a[i][n] - tong) / a[i][i]; } x = y; } PHƯƠNG PHÁP SỐ-Bài 3 3 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (3) 3. MỘT SỐ ĐỊNH LÍ VỀ NGHIỆM Định lí 4.1: Hệ Ax = b có nhiều nhất một nghiệm (tức là, nghiệm là duy nhất nếu tồn tại ) nếu và chỉ nếu hệ thuần nhất tương ứng Ax = 0 chỉ có nghiệm “tầm thường” x= 0. PHƯƠNG PHÁP SỐ-Bài 3 4 Định lí 4.2: Bất kì hệ PTTT thuần nhất nào với số phương trình ít hơn số ẩn đều có nghiệm không tầm thường (khác 0) Định lí 4.3 Nếu A là một ma trận cấp m × n và hệ Ax = b có nghiệm với mọi vectơ m chiều b, thì m ≤ n.. Định lí 4.4 Cho A là ma trận cấp n × n. Các khẳng định sau đây là tương đương: (i) Hệ thuần nhất Ax = 0 chỉ có nghiệm tầm thường x = 0. (ii) Với mọi vế phải b, hệ Ax = b luôn có nghiệm. (iii) A là ma trận khả nghịch HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (4) 4. LỜI GIẢI BẰNG SỐ CỦA HỆ PTTT Ax = b Xét các hệ PTTT có số phương trình đúng bằng số ẩn.  Nếu A khả nghịch thì hệ có duy nhất nghiệm đối với mọi b Có hai loại phương pháp giải: PHƯƠNG PHÁP SỐ-Bài 3 5 Lặp: xuất phát với một xấp xỉ ban đầu và dùng một thuật toán đã được lựa chọn phù hợp, sẽ đưa ra các xấp xỉ liên tiếp ngày càng tốt hơn.  Chỉ nhận được nghiệm gần đúng. Trực tiếp: là các phương pháp đưa ra nghiệm chính xác bởi một số hữu hạn các phép toán số học sơ cấp (khử Gauss).  Các phương pháp trực tiếp thực hiện trên máy tính thường không đưa đến các nghiệm chính xác do các sai số làm tròn. HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (5) 5. CÁC PHÉP BIẾN ĐỔI TƯƠNG ĐƯƠNG HỆ PTTT Định lí 4.7 Cho Ax = b là hệ PTTT, các xử lí hệ này bằng các phép toán sau đây dẫn tới các hệ tương đương: (i) Nhân một phương trình với một hằng số khác 0 (ii) Cộng một phương trình đã nhân với một hằng số vào một phương trình khác (iii) Đổi chỗ hai phương trình Nhận xét: Nếu với k và i xác định mà akk ≠ 0 thì chúng ta có thể khử ẩn xk từ phương trình thứ i bất kì bằng cách cộng phương trình thứ k đã nhân với –(aik /akk) vào phương trình thứ i. Khi đó hệ phương trình hệ quả sau là tương đương với hệ ban đầu . bÃx ~  PHƯƠNG PHÁP SỐ-Bài 3 6 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (6) 6. PHÉP KHỬ GAUSS: đưa hệ Ax = b, về hệ tam giác trên: Gọi W là ma trận cấp n x (n + 1) chứa ma trận vuông A ở n cột đầu và vectơ b ở cột cuối cùng. Với k = 1, , n–1 lặp các công việc: Tìm hàng trụ i ≥ k gần hàng k nhất sao cho wik ≠ 0 Nếu không có hàng i như vậy thì dừng (A không khả nghịch) Nếu tìm được hàng i, thì đổi chỗ hàng i với hàng k trụ wkk Với i = k+1, ..., n, lặp các công việc sau (khử các pt dưới trụ) mi = wik / wkk là hệ số nhân cho hàng i Với j = k+1, , n+1, lặp các công việc sau (biến đổi hàng i) wij = w ij – mi x wkj Nếu wnn = 0 thì dừng (A không khả nghịch) Nếu trái lại  hệ tam giác trên Ux = y với và yj = wi,n+1 i = 1, , n          ji0 jiw u ij ij PHƯƠNG PHÁP SỐ-Bài 3 7 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (7) 6. PHÉP KHỬ GAUSS (tiếp) Ví dụ : Giải hệ PTTT sau: 2x – 3y = 3 4x – 5y + z = 7 2x – y – 3z = 5 Giải: Bước 1: Dùng trụ đầu tiên của hệ để khử những hệ số bên dưới trụ đó: m2 = 4/2 = 2, m3 = 2/2 = 1  lấy pt 2 trừ đi 2 lần pt 1, và pt 3 trừ đi 1 lần pt 1 ta được hệ: 2x – 3y = 3 1 y + z = 1 2 y – 3z = 2 Bước 2: Dùng trụ thứ hai để khử hệ số bên dưới nó: m3 = 2/1 = 2 lấy pt 3 trừ đi 2 lần pt 2 ta được hệ 2x – 3y = 3 1y + z = 1 –5z = 0 giải ngược từ dưới lên ta được: z = 0, y = 1, x = 3 PHƯƠNG PHÁP SỐ-Bài 3 8 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (8) void giaiGauss(vector > a) { unsigned i, j, k, n = a.size(); double t, m; for (k = 0; k < n–1; k++) { // BUOC KHU XUOI HE TAM GIAC TREN while (a[i] [k] == 0 && i < n) i++; // Tim hang i gan k nhat: a[i] [k]≠0 if (a[i] [k] != 0 && i != k) { for (j = i; j <= n; j++) { // Doi cho hang i va hang k t = a[i] [j]; a[i] [j] = a[k] [j]; a[k] [j] = t; } } else if(a[i] [k] == 0) { // Ma tran A suy bien cout<<"He khong co duy nhat nghiem"<< endl; exit(1); } for (i = k+1; i <= a.size()–1; i++) { // Tao ma tran dang tam giac tren m = a[i] [k] / a[k] [k]; for (j = k+1; j <= a.size(); j++) a[i] [j] = a[i] [j] – m * a[k] [j]; } } // BUOC QUET NGUOC - GIAI HE TAM GIAC TREN }  Độ phức tạp thời gian của thuật toán là O(n3) PHƯƠNG PHÁP SỐ-Bài 3 9 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (9) 7. GIẢI HỆ PTTT DẠNG 3 ĐƯỜNG CHÉO Ma trận A = (aij) cấp n có dạng ba đường chéo nếu aij = 0 bất cứ khi nào |i – j| > 1. Hệ PTTT 3 đường chéo có dạng Ví dụ: hệ PTTT bên phải là hệ ba đường chéo nnn1nn 1nn1n1n1n2n1n 3433323 2322212 12111 bxdxl bxuxdxl bxuxdxl bxuxdxl bxuxd         0 0 0 1 2100 1210 0121 0012                   PHƯƠNG PHÁP SỐ-Bài 3 10 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (10) 7. GIẢI HỆ PTTT DẠNG 3 ĐƯỜNG CHÉO Cho trước các hệ số li, di, ui và vế phải bi của hệ ba đường chéo lixi–1 + dixi + uixi+1 = bi, i = 1 , . . . , n (với l1 = un = 0) Bước 1: Với k = 2, , n lặp các công việc sau đây Nếu dk−1 = 0, đó là dấu hiệu sai và dừng Trái lại, m = lk / dk–1 và tiếp tục tính dk = dk – m*uk–1 bk = bk – m*bk–1 Nếu dn = 0, đó là dấu hiệu sai và dừng Trái lại: xn = bn / dn Bước 2: Với k = n–1, , 1 lặp các công việc sau Nhận xét: Độ phức tạp thời gian của thuật toán là O(n) k 1kkk k d *xub x    PHƯƠNG PHÁP SỐ-Bài 3 11 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (11) 7. GIẢI HỆ PTTT DẠNG 3 ĐƯỜNG CHÉO (tiếp) for (unsigned k=1; k<n; k++) { if (d[k-1] == 0) { cout << "Khong giai duoc" << endl; return; } else { double m = l[k] / d[k-1]; d[k] = d[k] - m * u[k-1]; b[k] = b[k] - m * b[k-1]; } } If (d[n-1] == 0) { cout << "He khong co nghiem duy nhat" << endl; return; } else { x[n-1] = b[n-1] / d[n-1]; for(int i = n - 2; i >= 0; i--) x[i] = (b[i] - u[i] * x[i+1]) / d[i]; } for (unsigned i = 0; i < n; i++) cout << "x[" << i + 1 << "]=" << x[i] << endl; PHƯƠNG PHÁP SỐ-Bài 3 12 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (12) 8. CHIẾN LƯỢC XOAY Nhận xét: Thuật toán khử Gauss đưa ra lời giải sai cho ví dụ sau 0.0003 x1 + 1.566 x2 = 1.569 0.0003 x1 + 1.566 x2 = 1.569 0.3454 x1 – 2.436 x2 = 1.018 – 1805.42 x2 = – 1807.46 PHƯƠNG PHÁP SỐ-Bài 3 13 Nếu lấy a12 làm trụ, ta có nghiệm đúng! Phân tích Trong bước khử Gauss đưa hệ Ax = b về hệ tam giác trên Ux = y sử dụng ma trận W cấp n × (n + 1) ta đã tìm các hàng trụ ở các bước lặp k = 1, , n–1 theo tiêu chí: hàng trụ i ≥ k gần hàng k nhất sao cho wik ≠ 0 Hệ này có nghiệm đúng là x1= 10, x2 = 1, nhưng phép khử Gauss đưa ra x2 =1.001; x1 = 4.780 do trụ a11 = 0.0003 “rất nhỏ” so với a12 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (13) 8. CHIẾN LƯỢC XOAY (tiếp) Xoay theo tỉ lệ riêng: Ở lần lặp thứ k của phép khử Gauss • Tính “cỡ” di của hàng i của ma trận A, với i = k, . . . , n. • Chọn ra phương trình xoay là một trong các phương trình trong n – k ứng cử viên còn lại mà hệ số của xk có giá trị tuyệt đối lớn nhất so với cỡ của phương trình đó Trong thuật toán Gauss, điều này có nghĩa là số nguyên j được lựa chọn là số nguyên nằm giữa k và n (thường là số nhỏ nhất) sao cho tỉ lệ riêng |a|maxd ij njk i   PHƯƠNG PHÁP SỐ-Bài 3 14 i ik j jk d |w| d |w|  với mọi i = k, , n HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (14) 9. PHÂN TÍCH MA TRẬN THÀNH TÍCH − GIẢI HỆ PTTT Sử dụng phép khử Gauss, không đổi chỗ các hàng, giải hệ Ax = b Bước 1: Phân tích A = LU, với U là ma trận tam giác trên, L là ma trận tam giác dưới, chứa các hệ số nhân của lần lặp thứ k Bước 2: - Giải hệ tam giác dưới tìm được y Vector p = (pi), i=1, ,n dùng để lưu thay đổi thứ tự của các hàng - Giải hệ tam giác trên Ux = y tìm được x  ,n,ib ip 1 ),( *bLy          ki k i kimw l ikik ik , , 0 1 , kink aam k kk k ikik    ;1,,1 ,/ )1()1(  PHƯƠNG PHÁP SỐ-Bài 3 15 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (15) 9. PHÂN TÍCH MA TRẬN THÀNH TÍCH − GIẢI HỆ PTTT THUẬT TOÁN THÉ XUÔI VÀ KHỬ NGƯỢC PHA 1: Biến đổi ma trận A ban đầu về tích các ma trận tam giác Với k = 1, , n–1 lặp các công việc: Tìm hàng trụ i ≥ k gần hàng k nhất sao cho nó có tỉ lệ riêng lớn nhất Nếu không có hàng i như vậy thì dừng (A không khả nghịch) Trái lại, thì đổi chỗ 2 hàng i và k, nếu i ≠ k, lưu thứ tự p mới của b Với i = k+1, ..., n, lặp các công việc (khử các phần tử dưới trụ): m = aik / akk (= Giá trị cần khử ở hàng i / Điểm trụ ở hàng trụ k) Với j = k+1, , n-1, lặp các công việc (giữ nguyên vế phải): aij = aij – m x akj Gán aik = m PHƯƠNG PHÁP SỐ-Bài 3 16 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (16) 9. PHÂN TÍCH MA TRẬN THÀNH TÍCH - GiẢI HỆ PTTT THUẬT TOÁN THÉ XUÔI VÀ KHỬ NGƯỢC PHA 2: Giải hệ tam giác dưới Ly = b* tìm y (L có đường chéo =1) Với k = 1, , n lặp lại việc gán     1k 1j jkjkpk yaby kk n 1kj jkjk k a xay x     PHƯƠNG PHÁP SỐ-Bài 3 17 Giải hệ tam giác trên Ux = y tìm được x Với k = n, n–1, , 1 lặp lại việc gán HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (17) 9. PHÂN TÍCH MA TRẬN THÀNH TÍCH − GIẢI HỆ PTTT PHƯƠNG PHÁP SỐ-Bài 3 18 Ví dụ 1: Giải hệ u + 2v = 5 4u + 9v = 21 với hệ số nhân là 4, lưu trong L. Vế phải được dùng để tìm y Giải hệ tam giác dưới Ly = b:                   1 5 21 5 14 01 yy Giải hệ tam giác trên Ux = y:                   1 3 1 5 10 21 xx        94 21 A HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (18) 9. PHÂN TÍCH MA TRẬN THÀNH TÍCH − GIẢI HỆ PTTT Ví dụ 2: Giải hệ Ax = b PHA 1 PHƯƠNG PHÁP SỐ-Bài 3 19                                                                                      13 01 001 L 5 3 7 5-3- - 154 3- l rt hàng 2k 5 3 7 - 154 l l 5312 3032 7154 *A 2i 5 4 d |a| d |a| max 3 5 3 d 5312 7154 3032 Axét1k 2 1 2 1 2 1 2 1 2 1 2 1 32 2 7- 2 3 2 1 2 1 2 1 2 1 2 1 31 2 1 21 2 21 i 1i 3i1 2 i u HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (19) 9. PHÂN TÍCH MA TRẬN THÀNH TÍCH − GIẢI HỆ PTTT Ví dụ 2: PHA 2 Giải hệ Ly = b* với PHƯƠNG PHÁP SỐ-Bài 3 20 Giải hệ Ux = y với                         0 7 1 01 001 2 1 y 5 3 7 3- *b|L 2 1 2 1                           0 1 3 500 0 154 x 0 7 -y|U 2 1 2 1 2 1 - HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (20) 9. PHÂN TÍCH MA TRẬN THÀNH TÍCH − GiẢI HỆ PTTT Nhận xét 1. Để nhận được nội dung cuối cùng của ma trận làm việc W, phép khử Gauss cần O(n3) phép cộng và bằng ấy phép nhân / chia. Thuật toán phân tích A = LU giải hệ PTTT chỉ cần n2 phép nhân / chia và n(n–1) phép cộng 2. Thuật toán phân tích A = LU giải hệ PTTT có thể dùng giải hệ phương trình Ax = b với các vế phải khác nhau 3. Để đổi chỗ hàng i cho hàng j của ma trận W ta sử dụng phép nhân trái ma trận hoán vị Pij (nhận được từ ma trận đơn vị I cùng cấp, bằng cách đổi chỗ các hàng i và j của I) với W. Ví dụ P32 và           010 100 001                      3 5 1 5 3 1                                300 560 142 560 300 142 010 100 001 PHƯƠNG PHÁP SỐ-Bài 3 21 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (21) 10. PHÉP LẶP ĐIỂM CỐ ĐỊNH GIẢI HỆ PTTT Nếu đặt thì hệ Ax = b còn có thể viết ở dạng phương trình ma trận b – Ax = 0 Giả sử tồn tại ma trận C xấp xỉ ma trận nghich đảo A-1 theo nghĩa ||I – CA|| < 1 Khi đó b – Ax = 0  C(b – Ax) = 0  x = x + C(b – Ax) Chọn g(x) = x + C(b – Ax) = (I – CA)x + Cb x x x x b b b b và n 2 1 n 2 1                            ,, a...aa ............ a...aa a...aa A mnm2m1 2n2221 1n1211              PHƯƠNG PHÁP SỐ-Bài 3 22 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (22) 10. PHÉP LẶP ĐIỂM CỐ ĐỊNH GIẢI HỆ PTTT Nhận xét: g(x) – g(y) = Cb + (I– CA)x – [Cb + (I – CA)y] = (I – CA)(x – y) Vậy ||g(x) – g(y)|| < ||I – CA|| ||x – y|| Chứng tỏ x-y được g co lại với tỉ lệ co K=||I – CA||< 1,  phép lặp x(m +1) = x(m) + C(b – Ax(m)) = Cb + (I – CA)x(m) , m = 0, 1, 2, . khởi đầu từ bất kì x(0) nào, cũng sẽ hội tụ về nghiệm duy nhất ξ của PT Ax = b với sai số ở mỗi bước lặp giảm xuống ít nhất theo tỉ lệ K = || I – CA||.  Tìm các ma trận C xấp xỉ ma trận nghịch đảo A-1? PHƯƠNG PHÁP SỐ-Bài 3 23 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (23) 10. PHÉP LẶP ĐIỂM CỐ ĐỊNH GIẢI HỆ PTTT Thuật toán Phép Lặp Điểm Cố Định đối với các hệ PTTT Cho hệ phương trình tuyến tính Ax = b bậc n. 1. Lấy một ma trận C bậc n thỏa mãn các điều kện (i) Đối với vectơ r cho trước, vectơ Cr tính được “một cách dễ dàng” (ii) Với một chuẩn ma trận nào đó, ||I – CA|| < 1 2. Lấy một vectơ n chiều x(0), (chẳng hạn, x(0) = 0) Với m = 0, 1, 2, , cho đến khi thỏa mãn điều kiện lặp, tính x(m+1) = x(m) + C(b – Ax(m)) PHƯƠNG PHÁP SỐ-Bài 3 24 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (23) 10. PHÉP LẶP ĐIỂM CỐ ĐỊNH GIẢI HỆ PTTT Sự hội tụ và sai số: ||x(m+1) - x(m|| = ||I – CA|| ||x(m) - x(m-1) || = K ||x(m) - x(m-1) || == Kn ||x(1) - x(0) | → 0 khi m → ∞. Bỏ qua sai số làm tròn, dãy x(0), x(1), x(2), hội tụ về nghiệm đúng ξ của hệ PTTT đã cho. Với K = ||I – CA|| < 1 ta suy ra được sai số ||x(m) – ξ|| = ||g(x(m-1))– g(ξ)|| ≤ K ||x(m-1) – ξ|| (*) Do ||x(m-1) – ξ|| ≤ ||x(m) – ξ|| + ||x(m-1) - x(m) || (**) → ||x(m) – ξ|| ≤ ||x(m) – x(m –1)|| * K / (1-K) PHƯƠNG PHÁP SỐ-Bài 3 25 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (24) 10. PHÉP LẶP ĐIỂM CỐ ĐỊNH GIẢI HỆ PTTT Tiêu chuẩn kết thúc lặp (a): ||x(m) – x(m –1)|| < ε hoặc (b): ||x(m) – x(m –1)|| / || x(m) || < ε với ε cho trước hoặc (c): m > M với M cho trước Vì với K = ||I – CA|| < 1 yếu tố ||x(m) – x(m –1)|| < ε bao hàm điều ||x(m) – ξ|| < ε Tiêu chuẩn cuối luôn có mặt trong bất cứ chương trình lặp nào để tránh việc lặp không kết thúc PHƯƠNG PHÁP SỐ-Bài 3 26 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (25) 10. PHÉP LẶP ĐIỂM CỐ ĐỊNH GIẢI HỆ PTTT Phép lặp Jacobi: Giả sử A là ma trận chéo trội thực sự, nghĩa là với i = 1, 2, , n Gọi D = diag(a11, a22, . . . ,ann) là ma trận đường chéo của A. Xét Chứng tỏ D-1 là một nghịch đảo xấp xỉ của A.  Sơ đồ lặp x(m +1) = x(m) + D–l (b – Ax(m)) ↔ x(m +1) = A’x(m) + b’ và công thức lặp tính các thành phần của nghiệm, j= 1, , n    ij ijii ||a| |a               ij iiji ij iiji i ||a/||a||a/||a 11 maxmax||ADI|| 1 (*)a/xabx ii ij (m) jjii )(m i            1 PHƯƠNG PHÁP SỐ-Bài 3 27 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (26) 10. PHÉP LẶP ĐIỂM CỐ ĐINH GIẢI HỆ PTTT Thuật toán lặp Jacobi: từ công thức lặp (*) ta suy ra thuật toán sau 1. Cấu tạo ma trận lặp A’ và b’ với 2. Lặp x(m+1) = A’ x(m) + b’ với nghiệm xấp xỉ ban đầu x(0) tùy ý 3. Điều kiện dừng: khi một trong 3 tiêu chuẩn (a), (b), (c) thỏa mãn .                                              nn n nn n nn n n n a b a b a b a a a a a a a a a a a a      22 2 11 1 21 22 2 22 21 11 1 11 12 0 0 0 b' A' PHƯƠNG PHÁP SỐ-Bài 3 28 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (27) 10x1 + x2 + x3 = 12 x1 + 10x2 + x3 = 12 x1 + x2 + 10x3 = 12 với x(0) = 0 và ma trận lặp Dãy này hội tụ về nghiệm đúng [1 1 l]T. Do x1 (0) x2 (0) x3 (0) 0 0 0 1.2 1.2 1.2 0.96 0.96 0.96 1.008 1.008 1.008 0.9984 0.9984 0.9984 1.00032 1.00032 1.00032 0.999936 0.999936 0.999936   0.000096ξxADI (6)1   ||||.max|||| 20101101i                          2.1 2.1 2.1 01.01.0 1.001.0 1.01.00 b' và A' PHƯƠNG PHÁP SỐ-Bài 3 29 10. PHÉP LẶP ĐIỂM CỐ ĐỊNH GIẢI HỆ PTTT Ví dụ, giải hệ sau bằng phép lặp Jacobi: NHẬN XÉT: Sai số này quá lớn vì thực ra ||ξ – x(6)|| = 0.000064 void lapJacobi (double eps, vector > a) { unsigned i, j, k , n = a.size(); vector xt(n,0), xs(n,0); // Khởi tạo vector nghiệm bool t; // Biến logic kiểm tra sai số k = 1; // Biến đếm số lần lặp do { for(i = 0; i <= n –1; i++) { xs[i] = 0; for (j = 0; j <= n – 1; j++) if (j != i) xs[i] = xs[i] – a[i][j] * xt[j]; xs[i] = (xs[i] + a[i][n]) / a[i][i]; } // Kiểm tra đ/k dừng t t = max {fabs(xs[i] – xt[i])} <= eps } while ( ( !t ) && (k <=100) ); } PHƯƠNG PHÁP SỐ-Bài 3 30 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (28) 10. PHÉP LẶP ĐIỂM CỐ ĐỊNH GIẢI HỆ PTTT Phép lặp Gauss–Seidel, hay phương pháp thay thế liên tiếp: Phân tích A = L + D + U với L là ma trận tam giác dưới thực sự, U là ma trận tam giác trên thực sự. Không mất TQ, giả sử ma trận đường chéo D khả nghịch, tức là aii≠ 0 với mọi i. Chọn C = (L + D) -1 làm nghịch đảo xấp xỉ của A (nới lỏng đ/k của C), nếu  Sơ đồ lặp x(m+1) = x(m) + ( L+ D)–1(b – Ax(m)) . và công thức lặp tính các thành phần của nghiệm 1||||||)(|| 11   ADI ADLI bULxD   )()()( mmm xx 11 n...,2,1,i, a bxaxa x ii ij ij i (m) jji 1)(m j 1)(m i          ji PHƯƠNG PHÁP SỐ-Bài 3 31 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (29) 10. PHÉP LẶP ĐIỂM CỐ ĐỊNH GIẢI HỆ PTTT Thuật toán lặp Gauss–Seidel : 1. Cấu tạo ma trận lặp A’ và vectơ b’ như trong phép lặp Jacobi, với giả thiết aii ≠ 0 với mọi i 2. Chọn một nghiệm xấp xỉ ban đầu x(0) 3. Với m = 1, 2, . . . , cho đến khi thỏa mãn, lặp các công việc sau: Với i = 1, , n, lặp các công việc sau: Xét ví dụ giải hệ trên, ta có bảng nghiệm x1 (0) x2 (0) x3 (0) 0 0 0 1.2 1.08 0.972 0.9948 1.0033 1.00019 0.99965 1.000016 1.0080033 i m j n i j ij m j i j ij m i b x a x a x ' ' ' ) 1 ( 1 ) ( 1 1 ) (           PHƯƠNG PHÁP SỐ-Bài 3 32 NHẬN XÉT: Phép lặp Gauss-Seidel hội tụ nếu ma trận hệ số A là chéo trội / tam giác dưới / xác định dương. Khi đó nó có tốc độ hội tụ cao hơn phép lặp Jacobi void lapSeidel (double eps, vector > a) { unsigned i, j, k, n = a.size(); vector xt(n,0), xs(n,0); // Khởi tạo vector nghiệm bool t; k = 1; // Đếm số lần lặp do { for (i = 0; i <= n –1; i++) { xs[i] = 0; for (j = 0; j <= n – 1; j++) { if ( j < I ) xs[i] = xs[i] – a[i][j] * xs[j]; else if ( j > I ) xs[i] = xs[i] – a[i][j] * xt[j]; } xs[i] = (xs[i] + a[i][n]) / b[i][i]; } // Kiểm tra đ/k dừng t = max {fabs(xs[i] – xt[i])} <= eps } while ( ( !t ) && (k <=100) ); } PHƯƠNG PHÁP SỐ-Bài 3 33 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH (30) 11. TỔNG KẾT CÁC PHƯƠNG PHÁP SỐ GiẢI HỆ PTTT Trực tiếp (khử Gauss và các cải thiện) : 1. Là các phương pháp giải đúng  không có sai số thường không đưa đến các nghiệm chính xác do mắc phải sai số làm tròn số khi tính toán. 2. Áp dụng được cho mọi hệ vuông Lặp: 1. Là phương pháp lặp điểm cố định, xuất phát từ một xấp xỉ ban đầu và dùng một thuật toán lặp phù hợp, đưa ra các nghiệm xấp xỉ liên tiếp ngày càng tốt hơn. 2. (i) Tốc độ hội tụ phụ thuộc vào đoán nhận ban đầu x(0) (ii) Phải xây dựng ma trận lặp thỏa mãn điều kiện hội tụ  Thường áp dụng cho các hệ PTTT có ma trận hệ số thưa hay có các tính chất đặc biệt như chéo trội, xác định dương, PHƯƠNG PHÁP SỐ-Bài 3 34 PHƯƠNG PHÁP SỐ-Bài 3 SỬ DỤNG SOLVER GIẢI HỆ PTTT (1) 35 PHƯƠNG PHÁP SỐ-Bài 3 SỬ DỤNG SOLVER GIẢI HỆ PTTT (2) 36 PHƯƠNG PHÁP SỐ-Bài 3 1. Nhập các giá trị đoán nhận ban đầu cho x1, x2, , xn 2. Nhập các vế trái của các phương trình f1, f2, , fn vào các ô riêng được liên kết với giá trị của các ẩn 3. Chọn công cụ Solver (Data Tab  Solver) a. Đặt “Target Cell” cho ô có công thức = 0 b. Chọn “Equal to Value of 0” c. Đối với “By Changing Cells,” trỏ tới các ô chứa giá trị ban đầu của x1, x2, , xn d. Thêm các ràng buộc (các phương trình) e. Nhấn nút “Options” , chọn “Assume linear model” f. Nhấn nút “Solve” SỬ DỤNG SOLVER GIẢI HỆ PTTT (3) 37 TÍNH MA TRẬN NGHỊCH ĐẢO (1) Phương pháp Gauss-Jordan: để giải AA-1= I tìm A-1, ta tìm từng cột của A-1: Chi phí tính toán chỉ tăng gấp 3 do tất cả các hệ PTTT đều có chung ma trận vế phải A Ví dụ: [A i1 i2 i3 ] = Đem 1/2 hàng 1 + hàng 2     IiiixxxAAA 321321 1               100210 010121 001012              100210 01 2 1 1 2 3 0 001012               100 0110 001012 3 2 3 1 3 4 2 1 3 2 PHƯƠNG PHÁP SỐ-Bài 3 38 Đem 2/3 hàng 2 + hàng 3 TÍNH MA TRẬN NGHỊCH ĐẢO (2) Phương pháp Gauss-Jordan: tiếp tục phép khử - các hàng được cộng với các hàng trên chúng, để tạo ra các giá trị 0 trên các giá trị trụ (4/3 x hàng 3 + hàng2) 2/3 x hàng 2 + hàng 1                  1 3 2 3 1 3 4 00 4 3 2 3 4 3 0 2 3 0 001012                 1 3 2 3 1 3 4 00 4 3 2 3 4 3 0 2 3 0 2 1 1 2 3 002                 4 3 2 1 4 1 100 2 1 1 2 1 010 4 1 2 1 4 3 001  321 xxxI PHƯƠNG PHÁP SỐ-Bài 3 39 chia mỗi hàng cho giá trị tại điểm trụ của nó PHƯƠNG PHÁP SỐ-Bài 3 SỬ DỤNG HÀM MINVERSE TÌM A-1 CHÚ Ý: Công thức phải được nhập như công thức mảng: 1. Chọn mảng các ô chứa A-1. 2. Nhập hàm MINVERSE tính A-1 với các tham số là các ô chứa ma trận A 3. Nhấn phím F2, rồi nhấn tiếp CTRL+SHIFT+ENTER 40 PHƯƠNG PHÁP SỐ-Bài 3 SỬ DỤNG HÀM MINVERSE TÌM A-1 41 TÍNH ĐỊNH THỨC Phương pháp khử Gauss: tiến hành phép khử xuôi ma trận A để đưa nó về dạng tam giác trên U Khi đó det(A) = (−1)nu11u22 unn Ví dụ 105 2 1 41A 00 0 154 0 -0 154 3-1-2 03-2 154 3-1-2 15-4 032 2 7- 2 3 2 1 2 1                                                             )( )()det( 5- - æi_hµng§ A æi_hµng§ 2 1 2 1 ** * 1 PHƯƠNG PHÁP SỐ-Bài 3 42 n là số lần đổi hàng PHƯƠNG PHÁP SỐ-Bài 3 SỬ DỤNG HÀM MDETERM TÌM det(A) 43 BÀI TẬP 1. Bài tập tính toán: 4.2-5, 4.2-6, 4.2-8, 4.3-4, 4.4-2, 4.4-4, 5.3-7 2. Bài tập lập trình: Viết các hàm trong C++ giải hệ Ax = b bằng phép khử Gauss có sử dụng chiến lược xoay theo tỉ lệ riêng và phân tích ma trận vuông thành tích các ma trận tam giác; giải hệ có ma trận chéo trội bằng các phép lặp Jacobi và Gauss-Seidel. PHƯƠNG PHÁP SỐ-Bài 3 44