Sử dụng phương pháp lặp đơn và lặp Jacobi để giải hệ phương trình đại số tuyến tính

TRƯỜNG ĐẠI HỌC BÁCH KHOA HÀ NỘI CHƯƠNG TRÌNH ĐÀO TẠO KỸ SƯ CHẤT LƯỢNG CAO BỘ MÔN: PHƯƠNG PHÁP TÍNH SỬ DỤNG PHƯƠNG PHÁP LẶP ĐƠN VÀ LẶP JACOBI ĐỂ GIẢI HỆ PHƯƠNG TRÌNH ĐẠI SỐ TUYẾN TÍNH Giảng viên hướng dẫn: TS. Hà Thị Ngọc Yến Nhóm 6 – Lớp Kỹ sư Chất lượng cao K60 - Nguyễn Văn Lâm - Nguyễn Khánh - Hạ Minh Đức - Trần Minh Dũng Hà Nội, 03/2017 MỞ ĐẦU Giải phương trình Đại số tuyến tính bằng phương pháp Gauss và Gauss – Jordan vẫn còn tồn tại một số nhược điểm sau: - Khối lượng tính toán lớn. - Mặc dù là phương pháp giải đúng, nhưng vẫn tồn tại sai số do quy tròn trong các bước tính,sai số tổng hợp đôi khi sẽ khá lớn. Mặt khác, phương pháp này không thể khống chế sai số theo ý muốn tính toán. Vì vậy, người ta đã tìm kiếm những phương pháp gần đúng. phương pháp lặp đơn và lặp Jacobi ra đời để giảm bớt khối lượng tính toán và khống chế sai số tốt hơn. PHẦN I: LẶP ĐƠN I. Ý TƯỞNG CHUNG Giống như phương pháp lặp để giải phương trình x g() x , ta cũng tiến hàn lặp trên phương trình x x   bằng cách: Chọn giá trị nghiệm ban đầu x x0 , thay lần lượt vào phương trình trên, ta được: x1 x 0   x2 x 1   ... Ta cần xây dựng các công thức để tính toán điều kiện hội tụ và sai số. II. PHƯƠNG PHÁP THỰC HIỆN Bước 1: Từ phương trình x x   , chọn giá trị nghiệm ban đầu x x0 , thay liên tiếp vào phương trình, ta được dãy các vector lặp: x1 x 0   x x   2 1 x3 x 2   ... Bước 2: Làm liên tiếp như vậy, ta có biểu thức truy hồi: xk x k1   Bước 3: Gọi xk là vector lặp thứ k, nếu dãy vector lặp có giới hạn: * lim xk  x k Thì giới hạn ấy là nghiệm đúng của phương trình Ax b hay x x   . Lưu ý: Để tạo điều kiện thuận lợi cho lặp đơn, ta có một số cách để biến đổi phương trình từ Ax b về x x   , 11  1n  1      trong đó:       ;          m1   mn  n  Chẳng hạn, từ Ax b, cộng x vào hai vế, ta được: Ax x  b  x Từ đây, ta có hai cách khác nhau để tìm được ma trận  và  : ()A I x  b  x x b () I  A x Cách 1:  ()AI  Cách 2:  ()IA        b   b III. ĐIỀU KIỆN HỘI TỤ 1. Định nghĩa chuẩn của vector: a) Định nghĩa Với vector x, ta định nghĩa ‖x‖ là một số không âm thỏa mãn: 1. ‖x‖ ≥ 0, dấu “ = ” ⇔=0 2. ‖x‖ = || . ‖x‖ , k là hằng số 3. ‖x+y‖ ≤ ‖x‖ + ‖y‖ b) Một số chuẩn thường dùng với vector n |x | + ‖‖(1)=  i i1 n |x | + ‖‖(2)=  i i1 + ‖‖(∞) = max {|1|,|2|,,||} c) Sự hội tụ của vector về một vector: Dãy vector {X(k)} được gọi là hội tụ tới vector X khi ()k lim xi x i  i k hoặc limXXk   0 k 2. Định nghĩa chuẩn của ma trận: a) Định nghĩa: Với là ma trận thực, ta định nghĩa ‖‖ là 1 số không âm thỏa mãn: 1. ‖‖ ≥ 0, dấu “ = ” ⇔=0 2. ‖‖ = || . ‖‖ , k là hằng số 3. ‖+‖ ≤ ‖‖ + ‖‖ b) Một số chuẩn thường dùng của ma trận: n A max aij : Chuẩn cột 1 j  i1 1 n 2  2 A a : Chuẩn Euclide 2  ij  i, j  n A max aij : Chuẩn hàng  i  j1 3. Giới thiệu về chuẩn của ánh xạ tuyến tính liên tục a, Ánh xạ tuyến tính liên tục Phép ứng A từ không gian n vào không gian m được gọi là một ánh xạ tuyến tính nếu nó có các tính chất sau đây i, A()()(),, x y  A x  A y  x y  n ii, A()(),, x  A x     x  n n ta gọi các vectơ e11,0,...,0 , e 2  0,1,...,0 ,..., en   0,0,...,1 trong  là các vectơ đơn vị. Ánh xạ tuyến tính sẽ được biểu diễn qua ảnh của các vecto trục đơn vị theo công thức sau. Ax( ) xAe1 1  xAe 2 2   ...  xAen n  , * m Mỗi A ei  là một phần tử trong  cho nên nó sẽ là một bộ gồm m số, ký hiệu là ai1, a i 2 ,..., a in  ta thiết lập một ma trận A gồm n hàng m cột, với các cột là bộ số A ei  tức là T  A e1, A e 2  ,..., A en   , hay a11 a 21... an 1    a a... a T  12 22n 2  ... ... ... ...    a1m a 2 m... a nm  Ma trận này được gọi là ma trận của ánh xạ tuyến tính A. nếu coi mỗi vecto là một mà trận cột thì ta có thế viết lại công thức * dưới dạng đơn giản là: x1  x  A x  Tx với x  2      xn  Khi n=m thì T là một ma trận vuông cấp n. ánh xạ tương ứng với nó là một không gian n vào chính nó. b, Chuẩn của ánh xạ tuyến tính liên tục Nếu A: XY,.,. là ánh xạ tuyến tính liên tục thì ta định nghĩa chuẩn của A bởi:  XY   Ax A  sup Y x X x x  0 X Từ định nghĩa này ta dễ thấy các tính chất sau i, AAAsupx  sup  x x1YY x  1 XX ii, nếu A tuyến tính liên tục thì Ax  A., x  x  X   Y X iii, nếu A tuyến tính và tồn tại số dương M sao cho Ax  M., x  x  X   Y X Thì A liên tục và AM c, Liên hệ với phương pháp lặp đơn Ta cần giả phương trình F x  0, trong đó F(x) là một hàm trên không gian định chuẩn nào đó và 0 được hiểu là phần tử 0 của không gian này. ở đây ta xét không gian n thì 0 là vecto 0,0,...,0T .ta biến đổi phương trình về dạng tương đương x=G(x) rồi thực hiện lặp đơn. Ta có thể phát biểu định lý sau : Định lý: giả sử y=G(x) là một hàm liên tục trên không gian định chuẩn nào đó và phép lặp * * xn G x n1  n=1,2 hội tụ tới x với xuất phát ban đầu x0 khi đó x là nghiệm của phương trình x=G(x), tức là ta có x** G x  . Chứng minh: từ xn G x n1  ,với lưu ý là hàm G(x) liên tục ,ta có : ** limxn lim G ( x n1 )  G lim x n  1  x  G x n n  n     4. Định lý về sự hội tụ của phương pháp: Nếu có 1 chuẩn của  thoả mãn  1 thì phương trình Ax b có nghiệm duy nhất x* và ()p * dãy x()k hội tụ tới x theo chuẩn tương ứng. Chứng minh: a) Phương trình Ax b có nghiệm duy nhất: Xét hệ thuần nhất : 0 hay x   x x  x    x ()()()()p p p p Theo giả thiết:  q 1 ()p (1 q ) x  0 ()p x  0 ()p Hay nói cách khác: Nếu ma trận  thoả mãn  1 thì: hệ x  x chỉ có nghiệm tầm ()p thường, hay ệ x x   có nghiệm duy nhất. * b) dãy x()k hội tụ tới x theo chuẩn tương ứng. Ta có: limxk1  lim (   x k ) x x  limxk1    lim x k x x  * Do: limxk1  lim x k  x x x  Như vậy: x** x   Mặt khác, từ x** x   , ta được: *** xk x  x k1   x    x k  1  x  x  x***   x  x  q x  x k()()()p()p k1 p k  1 p Truy hồi theo k , ta được: x x**  qk x  x k ()()p1 p Do 0q  1 nên: x x*  0 khi k   k ()p * Hay lim xk  x là điều phải chứng minh. k IV. CÔNG THỨC SAI SỐ Ta chứng minh hai công thức sai số sau: q x x*  x  x k()p 1 q k k1 ()p qk x x*  x  x k ()p 1 q 1 0 ()p Theo công thức truy hồi: x x   k1 k xk x k1   Trừ theo vế: xk1 x k () x k  x k  1 x  x ( x  x )   x  x k1 k k k  1 k k  1 q xk  x k1 Suy ra: xk1 x k  q x k  x k  1 * Với l N bất kỳ, ta được: xk l x k  l 1  q x k  l  1  x k  l  2 l Sử dụng phép truy hồi, ta được: xk l x k  l 1  q x k  x k  1 Mặt khác, xét: x x k l k ()()()xk l  x k  l 1  x k  l  1  x k  l  2   x k  1  x k Lấy chuẩn toàn bộ: xk l x k ()()()xk l  x k  l 1  x k  l  1  x k  l  2   x k  1  x k xk l  x k  l 1  x k  l  1  x k  l  2   x k  1  x k Hay: xklk x  x klkl   x  1  x kl   1  x kl   2   x k  1  x k l Sử dụng công thức: xk l x k  l 1  q x k  x k  1 đã chứng minh ở trên, ta được: l1 l  2 l  1 l  2 xxqqk l k 1 xxqq k 1  k  1 qxx k  k  1 1 ql Như vậy: x x  q x  x k l k1 q k k 1 Cố định k và cho l   , ta được: q x x*  x  x k1 q k k1 l Cho k 1 trong công thức xk l x k  l 1  q x k  x k  1 : l xl1  x l  q x 1  x 0 k1 Thay l k 1, ta được: xk x k1  q x 1  x 0 q Thay vào công thức x x*  x  x , được: k1 q k k1 q qk x x*  x  x  x  x k1q k k1 1  q 1 0 Hay: qk x x*  x  x k 1 q 1 0 PHẦN II: LẶP JACOBI I. ĐẶT VẤN ĐỀ Lặp Jacobi cung cấp một cách xử lý khi chuẩn của  không thoả mãn điều kiện lặp đơn, bằng cách dựa vào tính chất chéo trội của ma trận. II. CÁC LOẠI LẶP JACOBI Tuỳ vào đề bài, ta có thể sử dụng lặp Jacobi đối với hai loại ma trận, đó là: - Jacobi đối với ma trận chéo hàng. - Jacobi đối với ma trận chéo cột. 1. Jacobi chéo hàng: a) Nội dung phương pháp: a11 a 1n    Cho ma trận:    , khi đó ta định nghĩa chéo trội hàng:   an1  a nn  n aii  a ij với i 1, n j1, j  i Với một ma trận chéo trội hàng thì hiển nhiên aii 0  i  1, n . Xét phương trình Ax b, với aii  0 , chia cả hàng thứ i cho aii . Các ma trận trở thành: a   b  a12 a 1 11a 1n   a  11   11  a   b  211   2 **a   a  A22 , b   22              an1 a n 2   b n   1    ann a nn   a nn  Khi đó, phương trình trở thành: A** x b x  A** x  b  x x () I  A** x  b  ** IA  a a  0 12  1n a a  11 11  a a  210   2n * a a  Trong đó:   22 22         an1 a n 2     0  ann a nn  Các bước tiếp theo, ta tiến hành như đối với lặp đơn. b) Điều kiện hội tụ: n Do: aii  a ij j1, j  1 n a 1   ij j1, j  1 aii  * 1  Do đó, bằng cách lặp Jacobi, ta luôn thu được mà trận thoả mãn điều kiện lặp đơn. c) Công thức sai số: Ta áp dụng công thức sai số giống như lặp đơn: q x x*  x  x k()p 1 q k k1 ()p qk x x*  x  x k ()p 1 q 1 0 ()p 2. Jacobi chéo cột a) Nội dung phương pháp a11 a 1n    Cho ma trận:    , khi đó ta định nghĩa chéo trội cột:   an1  a nn  n ajj  aij , j  1, n i1, i  j Với một ma trận chéo trội hàng thì hiển nhiên aii 0  i  1, n . Từ hệ phương trình: Ax b, ta đặt ẩn phụ: ' x1 a 11 x 1 x'  a x 2 22 2 ... ' xn a nn x n Hệ phương trình trở thành: A** x'  b a a  1 12 1n a a  22 nn  b  a a  1 211  2n   **  b2 Trong đó: Aa a , b    11 nn             bn  an1 a n 2   1  a11 a 22  Vậy nên: x''** x   , với ****I  A,   b a a  0 12  1n a a  22 nn  a a  210   2n * a a  Khi đó:   11 nn         an1 a n 2     0  a11 a 22  Ta tiến hành lặp đơn với ma trận  * như trên, chú ý sau bước tìm x* , ta cần trả lại ẩn x cho phương trình. b) Điều kiện hội tụ n Ta đã có: ajj  aij , j  1, n i1, i  j n a Suy ra: 1ij  *  1 (luôn thoả mãn điều kiện lặp đơn).  1 i1, i  j a jj c) Công thức sai số Ta thực hiện lặp đơn với ****I  A,   b như trên, tuy nhiên ta có đáp số là nghiệm gần đúng của x' chứ không phải nghiệm của hệ ban đầu. Nên ta cần đánh giá để liên hệ 2 sai số: n n n ''k ***** k  k    k  k  xx'  xxi  i '  axx ii i  i  min axx ii i  i  min axx ii .      1 i1 i  11i  n  i  1 1  i  n x''k  x * Như vậy với đầu vào sai số  :   xk  x* 1 min aii 1i  n Lưu ý: sau khi giải ra nghiệm x' ta phải đưa về nghiệm x là nghiệm ban đầu đề bài yêu cầu. PHẦN III: THUẬT TOÁN, CHƯƠNG TRÌNH VÀ VÍ DỤ I. SƠ ĐỒ THUẬT TOÁN 1. Lặp đơn Nhập , , , , ≔ 0 Type equation here. Đ Det(A)=0 In ra cần tìm phương pháp khác S A#At Đ s s : = ‖‖() < 1 : = − , : = Đ : = ‖‖() < 1 Đ Cần tìm chu ẩ n khác đ ể thực hiện lặp đơn s Đ In ra và 2. Lặp Jacobi Chương trình chính Nhập , , , , S Đ Đ S Đ Đ S không trội hàng, đkhông trội cột Lặp jacobi hàng Lặp jacobi cột S Kết thúc 3. Lặp Jacobi chéo trội hàng Nhập , , , , For i1, n tg:,; A i i j,,:,/; A i j  A i j tg b i :/; b i tg end S Đ In ra và 4. Lặp Jacobi chéo trội cột Nhập , , , , For i1, n duongcheo j:,; A j j j, A i , j :  A i , j / duongcheo j ; end s In ra và II. CHƯƠNG TRÌNH 1. Lặp đơn if det(A)==0 disp('can tim phuong phap khac de giai bai toan nay'); return; end k = 0; if size(A) ~= size(A') disp('A khong la ma tran vuong'); return; end % Kiem tra cach bien doi va chuan su dung thoa man alpha = eye(length(A)) + A; beta = -b; if norm(alpha, inf) < 1 p = inf; else if norm(alpha, 1) < 1 p = 1; else if norm(alpha, 2) < 1 p = 2; else alpha = eye(length(A)) - A; beta = b; if norm(alpha, inf) < 1 p = inf; else if norm(alpha, 1) < 1 p = 1; else if norm(alpha, 2) < 1 p = 2; else disp('Can tim cach bien doi khac,chuan khac de lap don'); return; end end end end end end %Qua trinh lap q = norm(alpha, p); TRAPPIST1 = true; while TRAPPIST1 k = k + 1; x = alpha*x0 + beta; if q*norm((x - x0),p)/(1-q) < epsilon TRAPPIST1 = false; else x0 = x; end end disp('sai so'); disp(q*norm((x - x0),p)/(1-q) ); disp('ma tran alpha'); disp(alpha); disp('ma tran beta'); disp(beta); disp('loai chuan duoc dung la'); disp(p); disp('so buoc lap la'); disp(k); disp('nghiem gan dung cua he phuong trinh la'); disp(x); 2. Lặp Jacobi a) Chương trình chính for i=1:length(A) sum1=0; for j=1:length(A) sum1=sum1+abs(A(i,j)); end sumh=sum1-abs(A(i,i)); if abs(A(i,i))<=sumh u1=0; break; else u1=1; end end for j=1:length(A) sum2=0; for i=1:length(A) sum2=sum2+abs(A(i,j)); end sumc=sum2-abs(A(j,j)); if abs(A(j,j))<=sumc u2=0; break; else u2=1; end end if u1==1 && u2==1 disp('ma tran A vua cheo troi cot vua cheo troi hang'); jacobihang(A, b, x0, epsilon); else if u1==1 && u2==0 disp('ma tran A cheo troi hang'); jacobihang(A, b, x0, epsilon); else if u1==0&&u2==1 disp('ma tran A cheo troi cot'); jacobicot(A, b, x0, epsilon); else disp('ma tran A khong cheo troi cot va khong cheo troi hang,can tim phuong phap khac'); end end end b) Lặp Jacobi chéo trội hàng function [x, k] = jacobihang(A, b, x0, epsilon) k = 0; for i = 1:length(A) temp = A(i,i); A(i,:) = A(i,:)/temp; b(i) = b(i)/temp; end alpha = eye(length(A)) - A; beta = b; %Qua trinh lap q = norm(alpha, inf); %kiem tra sai so trong vong lap TRAPPIST1 = true; while TRAPPIST1 k = k + 1; x = alpha*x0 + beta; if (q*norm((x - x0),inf)/(1-q) < epsilon) TRAPPIST1 = false; else x0 = x; end end disp('so buoc lap la'); disp(k); disp('nghiem cua he la'); disp(x); c) Lặp Jacobi chéo trội cột function [x, k] = jacobicot(A, b, x0, epsilon) k = 0; for j = 1:length(A) duongcheo(j) = A(j,j); A(:,j) = A(:,j)/duongcheo(j); end alpha = eye(length(A)) - A; beta = b; %Qua trinh lap q = norm(alpha, 1); TRAPPIST1 = true; epsilon = epsilon*(min(duongcheo(:))); while TRAPPIST1 k = k + 1; x = alpha*x0 + beta; if (q*norm((x - x0),1)/(1-q) < epsilon) TRAPPIST1 = false; else x0 = x; end end for j=1:length(A) x(j) = x(j)/duongcheo(j); end disp('so buoc lap la'); disp(k); disp('nghiem cua he la'); disp(x); III. MỘT SỐ VÍ DỤ 1. Lặp đơn: 0.7x1 + 0.2x 2 + 0.1x 3 + 0.3x 4 =1  0.1x1 + 0.8x 2 + 0.3x 3 + 0.2x 4 + 0.1x 5  3  0.3x1 + 0.1x 2 + 0.9x 3 + 0.1x 4 + 0.2x 5  5 0.2x + 0.1x + 0.2x + 0.7x + 0.1x 7  1 2 3 4 5 0.1x + 0.1x + 0.2x + 0.3x + 0.8x 1  1 2 3 4 5 T Với vector ban đầu là x0  0 0 0 0 0 và   0,01. Sau khi sử dụng chương trình, ta được bảng sau: k 1 2 3 4 5 6 7 8 9 10 11 12 13 14 ()k x1 1 -1.9000 -2.3200 -3.1260 -3.3000 -3.4899 -3.5352 -3.5757 -3.5854 -3.5935 -3.5953 -3.5968 -3.5972 -3.5974 ()k x2 3 0.5000 0.8200 0.1420 0.0921 -0.0974 -0.1351 -0.1848 -0.1979 -0.2100 -0.2136 -0.2164 -0.2172 -0.2178 ()k x3 5 4.0000 5.6300 5.7500 6.2096 6.2964 6.4130 6.4413 6.4687 6.4760 6.4821 6.4838 6.4850 6.4854 ()k x4 7 7.5000 9.0100 9.1960 9.5351 9.5935 9.6695 9.6853 9.7019 9.7057 9.7092 9.7100 9.7107 9.7109 ()k x5 1 -2.3000 -2.3700 -3.1530 -3.2410 -3.4299 -3.4646 -3.5093 -3.5197 -3.5299 -3.5325 -3.5348 -3.5354 -3.5359 Sai Số 63 29.7 14.67 7.254 4.1364 1.7092 1.0488 0.4473 0.2459 0.1091 0.0542 0.0249 0.0114 0.0054 Như vậy, ta được hệ nghiệm gần đúng sau 14 vòng lặp: x1  3,5974  x2  0,2178  x3  6,4854 x  9,7109  4 x5  3,5359 Sai số 0,0054 2. Lặp Jacobi chéo trội hàng 9x1 + 2x 2 + x 3 + 3x 4 + x 5 =9  2x1 + 8x 2 + x 3 + 2x 4 + x 5  9   x1 + 3x 2 + 7x 3 + x 4 + x 5   3 2x + x + 3x + 10x + 2x  9  1 2 3 4 5 2x1 + x 2 + 3x 3 + 2x 4 + 9x 5   9 T Với   0,01 và chọn x0  0 0 0 0 0 Sau khi chạy chương trình, ta được kết quả sau: k 1 2 3 4 5 6 7 8 9 ()k x1 1.0000 1.2087 1.2090 1.2001 1.1935 1.1942 1.1926 1.1938 1.1929 ()k x2 1.1250 1.2786 1.2671 1.2604 1.2560 1.2575 1.2562 1.2573 1.2564 ()k x3 -0.4286 -0.7821 -0.8795 -0.8890 -0.8909 -0.8884 -0.8898 -0.8885 -0.8896 ()k x4 -0.9000 -0.8839 -0.8341 -0.8140 -0.8132 -0.8111 -0.8127 -0.8115 -0.8125 ()k x5 -1.0000 -1.0044 -0.9535 -0.9310 -0.9295 -0.9271 -0.9287 -0.9274 -0.9285 Sai Số 9 2.8286 0.7785 0.1806 0.0533 0.0203 0.0132 0.0105 0.0086 Như vậy ta được hệ nghiệm gần đúng sau 9 vòng lặp: x1 1,1929  x2 1,2564  x3  0,8896 x  0,8125  4 x5  0,9285 Sai số 0,0086 . 3. Lặp Jacobi chéo trội cột 7x1 + 3x 2 + 2x 3 + x 4 +2 x 5 =-5  x1 + 8x 2 + x 3 +2x 4 +3x 5  9  2x1 + x 2 + 9x 3 + x 4 +2x 5  1 2x + x + 2x + 8x + x 1  1 2 3 4 5  x1 + 2x 2 + 2x 3 +3x 4 +9x 5   6 T Chọn   0,01, vector đầu là x0  0 0 0 0 0 Ta được kết quả sau 9 vòng lặp: k 1 2 3 4 5 6 7 8 9 k x '1 -5 -7.3889 -8.2268 -8.5034 -8.6010 -8.6378 -8.6521 -8.6573 -8.6593 k x '2 9 11.3532 12.0504 12.3002 12.3975 12.4349 12.4485 12.4535 12.4551 k x '3 1 2.5119 3.2808 3.5887 3.7011 3.7415 3.7561 3.7615 3.7634 k x '4 1 1.7480 2.0374 2.1148 2.1312 2.1354 2.1368 2.1376 2.1378 k x '5 -6 -8.1329 -8.9964 -9.3304 -9.4508 -9.4923 -9.5070 -9.5121 -9.5141 Sai số 176 73.0794 27.6559 9.9654 3.5524 1.2821 0.4682 0.1725 0.0627 x’ -1.2370     1.5569  Ta thu được nghiệm x là x   0.4182     0.2672    -1.0571  0,067 Với sai số là  0,009 min(aii )