›
›

Phương pháp tính - Chương 2: Hệ phương trình tuyến tính ax = b

ppt 31 trang vanle 3360
Download
Bạn đang xem 20 trang mẫu của tài liệu "Phương pháp tính - Chương 2: Hệ phương trình tuyến tính ax = b", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

  • pptphuong_phap_tinh_chuong_2_he_phuong_trinh_tuyen_tinh_ax_b.ppt

Nội dung text: Phương pháp tính - Chương 2: Hệ phương trình tuyến tính ax = b

  1. BỘ MƠN TỐN ỨNG DỤNG - ĐHBK PHƯƠNG PHÁP TÍNH – BG SINH VIÊN CHƯƠNG 2 HỆ PHƯƠNG TRÌNH TUYẾN TÍNH Ax = b TS. NGUYỄN QUỐC LÂN (2/2006)
  2. NỘI DUNG A- CÁC PHƯƠNG PHÁP CHÍNH XÁC 1- PHƯƠNG PHÁP KHỬ GAUSS (PHẦN TỬ TRỤ) 2- PHÂN TÍCH NHÂN TỬ A = LU 3- PHÂN TÍCH CHOLESKY B- CÁC PHƯƠNG PHÁP LẶP 1- LẶP JACOBI 2- LẶP GAUSS - SEIDEL C- SỐ ĐIỀU KIỆN – HỆ ĐIỀU KIỆN XẤU
  3. TỔNG QUAN Hệ n phương trình bậc 1 (tuyến tính), n ẩn → Dạng Ax = b: a11 a12  a1n b1 x1 a a  a b x A = 21 22 2n , b = 2 , x = 2 = x x T    1 n an1 an2  ann bn xn a11 a12  a1n 0 a  a Đơn giản: Hệ tam giác A = 22 2n Giải lùi 0 0  ann T Hàng i: hi = [ai1 ai2 ain] . Biến đổi sơ cấp trên hàng hi → hi + khj: Nhân hj với k rồi cộng xuống hi (chỉ hi thay đổi)
  4. PHƯƠNG PHÁP KHỬ GAUSS Giải thuật: Biến đổi sơ cấp trên hàng → A: trên → Giải lùi 2x1 − 2x2 + 3x3 =1 VD: Giải hệ Xây dựng ma trận mở rộng 6x1 − 7x2 +14x3 = 5 4x1 −8x2 + 30x3 =14 Khử cột 1 với hệ số khử m 2 − 2 3 1 1j A = A| b = 6 − 7 14 5 a a   1 j ij m1 j = Tổng quát : mij = 4 −8 30 14 a11 aii 6 m = = 3 2 − 2 3 1 h → h −3h 2 − 2 3 1 21 2 2 2 1 6 − 7 14 5 → 0 −1 5 2 4 h → h − 2h m31 = = 2 4 −8 30 14 3 3 1 − 4 24 12 2
  5. GIẢI LÙI & PHẦN TỬ TRỤ − 4 Khử cột 2 với hệ số khử: m = = 4 32 −1 2 − 2 3 1 2 − 2 3 1 h3 → h3 − 4h2 0 −1 5 2 → 0 −1 5 2 0 − 4 24 12 0 0 4 4 Giải lùi với hệ tam giác trên thu được: 2x1 − 2x2 + 3x3 =1 x3 = 4 4 =1 2 − x + 5x = 2 x = (2 − 5x ) (−1) = 3 x = 3 2 3 2 3 4x3 = 4 x1 = (1+ 2x2 − 3x3 ) 2 = 2 1 (1) (2) Điều kiện: Khử cột 1: a11 0 & Khử cột 2: a22 0 & (3) Giải lùi: a33 0 Phần tử trụ (pivot) akk 0
  6. KHỬ GAUSS VỚI LỆNH MAPLE > with(linalg); # Khởi động gĩi lệnh Đại số tuyến tính > A := matrix(2,3,[2, 3, 4, 1, 2, 3]); # Nhập ma trận > m21 := A[2,1]/A[1,1]; # Tính hệ số khử > A := addrow(A,1,2,–m21) ; # Cộng hàng h2 → h2 – m21h1 > A := swaprow(A,1,2) ; # Nếu cần thiết, đổi hàng h2  h1 > x := backsup(A) ; # Hệ đã ở dạng tam giác trên: Giải lùi > AA := gausselim(A); # Lệnh gộp khử Gauss tồn ma trận VD: Giải hệ x1 − x2 + 2x3 − x4 = −8 2x1 − 2x2 + 3x3 − 3x4 = −20 x1 + x2 + x3 = −2 x1 − x2 + 4x3 + 3x4 = 4
  7. KHỬ GAUSS VỚI MA TRẬN “LẺ”: PIVOT ĐƠN VỊ VD: Giải hệ với phép khử 2.08x1 −1.3x2 = 0.608 Gauss, làm trịn 3 chữ số lẻ) − 0.7x1 + 2.08x2 −1.3x3 = −0.152 − 0.7x2 + 2.08x3 −1.3x3 = −0.168 − 0.7x3 + 2.08x4 =1.116 2.08 −1.3 0 0 1.264 − 0.7 2.08 −1.3 0 − 0.152 Ab= 0 − 0.7 2.08 −1.3 − 0.168 0 0 − 0.7 2.08 1.116 1.006 0.636 y = 0.593 0.736
  8. THỰC TẾ TÍNH TỐN: VẤN ĐỀ LÀM TRỊN SỐ VD: Giải hệ trên máy tính với phép làm trịn 4 chữ cĩ nghĩa 0.003x1 + 59.14x2 = 59.17 (E1) T Nghiệm chính xác: [10, 1] 5.291x1 − 6.130x2 = 46.78 (E2 ) Quy tắc làm trịn trên máy tính: Làm trịn chữ số cĩ nghĩa 12,34567 =1,234567101 1,235101 =12,35 a Trụ khử: a = 0.003 0 21 11 m21 = =1763,67 1764 a11 Biến đổi cột một: (E2) → (E2) – m21(E1) 0.003x1 + 59.14x2 = 59.17 x2 =1.001 : Tại sao??? −104300x2 = −104400 x1 = −10
  9. PHÂN TÍCH NHÂN TỬ (MATRIX FACTORIZATIONS) Ma trận vuơng A phân tích được thành dạng LU 1 0 0 0 * * * * Hệ Ax = b (LU)x = b * 1 0 0 0 * * * A =  Ux = y (1) * * 1 0 0 0 * * : 2 hệ tam giác Ly = b (2) * * * 1 0 0 0 *     L U Giải hệ đầu Giải 2 hệ : Ly = b (2) tìm y; Ux = y (1) tìm x Nhân A x b Nhân L Nhân U y
  10. VÍ DỤ Giả sử ma trận A phân tích được thành dạng LU như sau: 3 − 7 − 2 2 1 0 0 0 3 − 7 − 2 1 − 3 5 1 0 −1 1 0 0 0 − 2 −1 1 A = =  6 − 4 0 − 5 2 − 5 1 0 0 0 −1 1 − 9 5 − 5 12 − 3 8 3 1 0 0 0 −1 Sử dụng phân tích LU trên giải hệ Ax = b = [–9 5 7 11]T Giải Ly = b tìm y Giải Ux = y tìm x
  11. PHÂN TÍCH NHÂN TỬ A = LU Quan sát: Ma trận khử L và ma trận kết quả U. Xét tích L.U 1 0 0 2 − 2 3 3 1 0  0 −1 5 = 2 4 1 0 0 4 Kết quả: Nếu quá trình khử Gauss diễn ra bình thường (khơng đổi hàng), ma trận A của hệ Ax = b phân tích được thành tích LU: A = LU với ➢ L (lower): ma trận tam giác dưới, đường chéo chính bằng 1, chứa các hệ số khử ở vị trí khử ➢ U (upper): ma trận tam giác trên, cũng là ma trận kết quả nhận được sau quá trình khử Gauss
  12. GIẢI THUẬT TÌM LU (CROUT – DOOLITLE) Phân tích LU với đường chéo chính L bằng 1 Khử Gauss (khơng đổi hàng). Các hệ số khử tạo L, ma trận kết quả: U 2 − 2 3 m21 = 3 1 0 0 VD: A = 6 − 7 14 : m = 2 L = 3 1 0 31 4 −8 30 m32 = 4 2 4 1 L (hoặc U) cĩ đchéo chính = 1 G/thuật Doolitle (Crout) i−1 For j = 1 to n: i = 1 → j uij = aij − likukj k=1 Tự xem j−1 i = j +1 → n 1 SGK/ 35 lij = aij − likukj u jj k=1
  13. MINH HOẠ GIẢI THUẬT DOOLITLE (ĐCHÉO L = 1) 2 − 2 3 1 0 0 VD: A = 6 − 7 14 L = 1 0 ,U = 0 4 −8 30 1 0 0 j = 1: i = 1 u11 = a11 a a i = 2 21 i = 3 31  21 =  31 = u11 u11 j = 2: i = 1 i = 2 u12 = a12 u22 = a22 − 21u12 a −  u i = 3 32 31 12  32 = u22 j = 3: i = 1 i = 2 i = 3 u13 = a13 u23 = a23 −  21u13 i = 3 u33 = a33 −  31u13 −  32u23
  14. PHÂN TÍCH CHOLESKY Tương tự phân tích LU nhưng gọn hơn “phân nửa”! Ma trận vuơng A (n hàng, n cột) : A = [ aij ] xác định dương n n T n T  x = x1,, xn  R , x 0 : x Ax = aij xi x j 0 i=1 j=1 1 1 0 VD: A = 1 5 − 2 : ma trận (đối xứng) xác định dương vì 0 − 2 2 T 3 T 2 2 2  x = x1, x2 , x3  R x Ax = x1 +5x2 + 2x3 + 2x1x2 − 4x2 x3 0 A: XĐD n định thức con (hvuơng) trên đ/chéo chính đều > 0
  15. GIẢI THUẬT CHOLESKY Định lý: Ma trận A đối xứng xác định dương Tồn tại ma trận tam giác dưới B thoả mãn : A = BBT i−1 2 A đối xứng: For i = 1 to n do bii = aii − bik k=1 i−1 For j = i+1 to n do 1 b ji = a ji − b jkbik bii k=1 Ax = b (BBT)x = b BTx = y & By = b: 2 hệ (như LU) A k0 xác định dương (chỉ đối xứng): A = BBT cĩ thể chứa số phức 2 hệ BTx = y & By = b: phức. Nhưng nghiệm x: thực!
  16. MINH HOẠ GIẢI THUẬT CHOLESKY 1 1 0 0 0 VD: A = 1 5 − 2 B = 0 0 − 2 2 i = 1 b11 = a11 a a j = 2 21 j = 3 31 b21 = b31 = b11 b11 2 i = 2 b22 = a22 −b21 a − b b j = 3 32 31 21 b32 = b22 i = 3 2 2 b33 = a33 −b31 −b32
  17. TỔNG QUAN PHƯƠNG PHÁP LẶP Chương 1: Phương pháp lặp đơn với phương trình f(x) = 0 x = (x) f (x) = 0 xn+1 = (xn ) (x) − (y) q x − y : '(x) q 1 Hệ Ax = b x = Tx + c = (x), T: ma trận, c: vectơ. Đkiện:  (x) – (y) qx – y Dãy lặp: x(n+1) = Tx(n) + c T n Chuẩn vectơ, ma trận: x = [x1, x2 xn] R , A = [aij ] n  A = max a x = max xi   ij  1 i n 1 i n j=1  n n  x = x 1  i A 1 = max  aij  i=1 1 j n i=1 
  18. VÍ DỤ ➢ Tính các chuẩn vectơ và ma trận 1 −1 3 2 x = 3 A = 12 x = − 3 A = − 2 4 3 6 13 x 1 = A 1 = − 2 − 5 6 −1 ➢ Vectơ nào trong số hai vectơ sau xấp xỉ tốt nhất theo chuẩn , chuẩn một nghiệm hệ phương trình . − 0.1 0.3 x1 + 2x2 + 3x3 =1 (2) (1) (2) x x x = − 7.4 , x = − 6.8 2x1 + 3x2 + 4x3 = −1 → . 1 (1) 5.2 4.7 3x1 + 4x2 + 6x3 = 2 x x Tch. chuẩn vectơ, chuẩn ma trận: Chuẩn tích tích chuẩn Ax A  x (x)− (y) = T(x − y) T  x − y : T 1
  19. LẶP JACOBI (0) T Với vectơ x = [0, 0, 0] , tìm vectơ 10x1 + 3x2 + x3 = 7.5 (k) nghiệm xấp xỉ x của phép lặp − 3x1 +10x2 − x3 = 9 Jacobi với hệ sau. Dừng: x(k) “giống” − x1 + 2x2 + 8x3 = −2.5 x(k-1) (khoảng 0.3) . So với nghiệm = [0.5, 1, -0.5]T 1/ Rút x trên đường chéo chính Đưa về dạng x = Tx + c 3 1 7.5 x = − x − x + = −0.3x − 0.1x + 0.75 1 10 2 10 3 10 2 3 3 1 9 x2 = x1 + x3 + = 0.3x1 + 0.1x3 + 0.9 x = Tx + c 10 10 10 1 2 2.5 x = x − x − = 0.125x − 0.25x − 0.3125 3 8 1 8 2 8 1 2
  20. CƠNG THỨC LẶP JACOBI 0 − 3 − 1 10 10 0.75 4 3 4 T = 3 0 1 , T = max ,  = 1, c = 0.9 10 10 10 8 10 − 0.3125 1 − 2 0 8 8 (1) (0) (0) x1 = −0.3x2 − 0.1x3 + 0.75 = 0.75 (0) (1) (1) (0) (0) 2/ Từ x tính x : x2 = 0.3x1 + 0.1x3 + 0.9 = 0.9 (1) (0) (0) − 0.3125 x3 = 0.125x1 − 0.25x2 − 0.3125 = (k+1) (k) (k) x1 = −0.3x2 − 0.1x3 + 0.75 (k) (k+1) (k+1) (k) (k) Tổng quát: x x : x2 = 0.3x1 + 0.1x3 + 0.9 (k+1) (k) (k) x3 = 0.125x1 − 0.25x2 − 0.3125 (k ) q (1) (0) Sai số: Như lặp đơn với q = ||T|| : x − x* x − x 1− q
  21. LẶP JACOBI KHƠNG BIẾN ĐỔI MA TRẬN A 10x1 + 3x2 + x3 = 7.5 10x1 = − 3x2 − x3 + 7.5 Đ/ chéo Hệ Ax = b: − 3x1 +10x2 − x3 = 9 10x2 = 3x1 + x3 + 9 chính − x1 + 2x2 + 8x3 = −2.5 8x3 = x1 − 2x2 − 2.5 (k+1) (k ) (k ) : Giữ đ/chéo chính ở vế trái 10x1 = − 3x2 − x3 + 7.5 (k+1) (k ) (k ) (→ x(k+1)) ; Chuyển số hạng 10x2 = 3x1 + x3 + 9 (k+1) (k ) (k ) cịn lại sang vế phải (→ x(k)) 8x3 = x1 − 2x2 − 2.5 (k +1) (k ) (k ) (k) 10x1 + 3x2 + x3 = 7.5 : Thay x vào các số hạng (k ) (k +1) (k ) − 3x1 +10x2 − x3 = 9 ngồi đường chéo chính. (k ) (k ) (k +1) Xem x(k+1) là ẩn. Giải→x(k+1) − x1 + 2x2 + 8x3 = − 2.5
  22. TÍNH TỐN & KẾT QUẢ LẶP JACOBI (k +1) (k ) (k ) 10x1 + 3x2 + x3 = 7.5 10x1 = −3x2 − x3 + 7.5 (k +1) (k ) (k ) Hệ: − 3x1 +10x2 − x3 = 9 Lặp Jacobi : 10x2 = 3x1 + x3 + 9 − x + 2x + 8x = −2.5 (k+1) (k ) (k ) 1 2 3 8x3 = x1 − 2x2 − 2.5 (k ) (k ) (k +1) − 3x2 − x3 + 7.5 k 0 1 2 x1 = 10 (k) x1 0.0 0.75 (k ) (k ) (k +1) 3x1 + x3 + 9 x (k) 0.0 0.9 x2 = 2 10 (k) x3 0.0 –0.3125 (k ) (k ) (k +1) x1 − 2x2 − 2.5 x(k)-x(k-1) 0.9 x3 = 8 Ưu điểm Lặp Jacobi: Giải các hệ “thưa” (chứa rất nhiều số 0) n M/trận đ/c trội nghiêm ngặt: a a  i ĐK đủ : T 1 ii  ij j=1, j i
  23. LẶP GAUSS – SEIDEL Tương tự lặp Jacobi nhưng với thơng tin cập nhật hố − 3x(k ) − x(k ) + 7.5 x(k +1) = 2 3 Dùng x(k) 1 10 10x1 + 3x2 + x3 = 7.5 Lặp (k ) (k ) để tính (k +1) 3x1 + x3 + 9 Hệ: − 3x1 +10x2 − x3 = 9 x2 = Jacobi 10 giá trị − x + 2x + 8x = −2.5 1 2 3 (k ) (k ) (k+1) (k +1) x1 − 2x2 − 2.5 của x x3 = 8 (k ) (k ) (k +1) − 3x2 − x3 + 7.5 x1 = 10 ❖ x1 (mới): dùng x2 (cũ), x3 (cũ) Gauss 3x(k +1) + x(k ) + 9 (k +1) 1 3 ❖ x (mới): dùng x (mới), x (cũ) x2 = 2 1 3 Seidel 10 (k +1) (k +1) ❖ x (mới): dùng x (mới), x (mới) (k +1) x1 − 2x2 − 2.5 3 1 2 x3 = 8
  24. LẶP GAUSS – SEIDEL: SƠ ĐỒ TÁCH MA TRẬN Trình bày dạng khác: Xem x(k+1) là ẩn và chuyển sang vế trái (k +1) (k ) (k ) 10x1 + 3x2 + x3 = 7.5 10x1 = −3x2 − x3 + 7.5 (k+1) (k +1) (k ) − 3x1 +10x2 − x3 = 9 − 3x1 +10x2 = x3 + 9 − x + 2x + 8x = −2.5 (k+1) (k+1) (k+1) 1 2 3 − x1 + 2x2 + 8x3 = − 2.5 x(0) k = 0 → b(k ) k = k +1 Giải hệ → x(k+1) Gauss - Seidel: Biết x(k) → Tính vế phải b(k) → Giải hệ ra x(k+1)
  25. LẶP GAUSS – SEIDEL: VÍ DỤ TÁCH MA TRẬN Xét ví dụ lặp Gauss – Seidel, x(0) = [0, 0, 0]T. Cơng thức lặp: (k+1) (k ) (k ) 10x1 = −3x2 − x3 + 7.5  (k +1) (k +1) (k ) (k ) − 3x1 +10x2 = x3 + 9  → b (k+1) (k+1) (k+1) − x1 + 2x2 + 8x3 = − 2.5  k 0 1 2 x b x b x b (k) x1 0.0 (k) x2 0.0 (k) x3 0.0 (k) (k-1) x -x  Phép lặp Thay hệ Ax = b bằng giải liên tiếp nhiều hệ
  26. TỔNG KẾT LẶP JACOBI & GAUSS – SEIDEL 10x1 + 3x2 + x3 = 7.5 − 3x1 +10x2 − x3 = 9 − x1 + 2x2 + 8x3 = −2.5 x = Tx + c 10x1 = − 3x2 − x3 + 7.5 Lặp Jacobi Lặp Gauss– Seidel 10x2 = 3x1 + x3 + 9 8x3 = x1 − 2x2 − 2.5 (k+1) (k ) (k ) (k+1) (k ) (k ) 10x1 = − 3x2 − x3 + 7.5 10x1 = − 3x2 − x3 + 7.5 (k+1) (k ) (k ) (k+1) (k +1) (k ) 10x2 = 3x1 + x3 + 9 10x2 − 3x1 = x3 + 9 (k+1) (k ) (k ) (k +1) (k +1) (k +1) 8x = x − 2x − 2.5 3 1 2 8x3 − x1 + 2x2 = − 2.5
  27. HỆ PHƯƠNG TRÌNH BỊ NHIỄU Minh hoạ: Giải 2 hệ phương trình và nhận xét x + 2y = 2 Hệ “gần” nhau, nghiệm 2x + 3.9y = 2 “xa” nhau! Do detA 0: x + 2y = 2 det Ai xi = ? 2x + 4.1y = 2 det A x + 2y = 2 2x + 3.9y = 2 2x + 4.1y = 2
  28. VÍ DỤ WILSON: Ax = b, detA = 1 10 7 8 7 32 1 7 5 6 5 23 1 Ax = b : A = , b = x = 8 6 10 9 33 1 7 5 9 10 31 1 32.1 22.9 A(x +x) = b +b : b +b = x +x = 33.1 30.9 10 7 8.1 7.2 7.08 5.04 6 5 (A + A) x + x = b : A + A = x' 8 5.98 9.89 9 6.99 4.99 9 9.98
  29. SỐ ĐIỀU KIỆN CỦA HỆ Ax = b x b •“Nhiễu” vế phải A(x + x) = b + b x b x A •“Nhiễu” vế trái (A + A)(x + x) = b x + x A + A x b x A A . A−1  ; A . A−1  x b x + x A + A –1 Số điều kiện:  (A) = ||A|| . ||A || đặc trưng cho độ “nhạy cảm” của nghiệm hệ Ax = b đối với những thay đổi dù rất nhỏ trên b hoặc A Hệ điều kiện xấu (ill – conditionned):  (A) >> 1
  30. VÍ DỤ VD Wilson: 10 7 8 7 7 5 6 5 A = 8 6 10 9 7 5 9 10 25 − 41 10 − 6 − 41 68 −17 10 A−1 = 10 −17 5 − 3 − 6 10 − 3 2 VD: Tính số điều kiện theo chuẩn vơ cùng  (A) của ma trận 1.003 58.09 A = 5.550 321.8
  31. PHƯƠNG PHÁP TÌM MA TRẬN NGƯỢC x y Tính ma trận ngược −1 A = = c1 c2  z t T ➢ Giải hệ phương trình A.c1 = e1 = [1 0] (vectơ đơn vị thứ nhất) trên máy tính bỏ túi Cột 1 của ma trận ngược ➢ Vẫn trong chế độ giải hệ phương trình, giải tiếp hệ A.c2 T -1 = e2 = [0 1] (vectơ đơn vị thứ nhì) Cột 2 của A ➢ Trường hợp ma trận cấp 3: Giải 3 hệ Ac1 = e1, Ac2 = e2, Ac3 = e3 với e1, e2, e3 lần lượt là 3 vectơ đơn vị Tìm được 3 –1 vectơ nghiệm c1, c2, c3: 3 cột của ma trận ngược A cần tìm