Hướng dẫn giải hệ phương trình bằng matlab năm 2024
|
Trong chương này chúng ta sẽ xét các phương pháp số để giải các phương trình đại số tuyến tính dạng: 11 1 12 2 1n n 1 21 1 22 2 2n n 2 Show
n1 1 n2 2 nn n n a x a x a x b a x a x a x b a x a x a x b ####### ⎧ + + ⋅ ⋅ ⋅ + = ####### ⎪⎪ + + ⋅⋅ ⋅ + = ####### ⎨ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ####### ⎪ ####### ⎪⎩ + + ⋅ ⋅ ⋅ + = Các phương trình này có thể viết gọn dưới dạng: [A] [x] = [b] Trong đó: [ ]11 12 1n 21 22 2n n1 n2 nn a a a a a a A a a a ####### ⎡ ⋅⋅⋅ ⎤ ####### ⎢ ⋅⋅⋅ ⎥ ####### = ⎢ ⎥ ####### ⎢ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅ ⋅⋅⎥ ####### ⎢ ⋅⋅⋅ ⎥ ####### ⎣ ⎦ [ ]1 2 n b b b b ####### ⎡ ⎤ ####### ⎢ ⎥ ####### = ⎢ ⎥ ####### ⎢ ⋅ ⋅ ⋅⎥ ####### ⎢ ⎥ ####### ⎣ ⎦ [ ]1 2 n x x x x ####### ⎡ ⎤ ####### ⎢ ⎥ ####### = ⎢ ⎥ ####### ⎢ ⋅ ⋅ ⋅⎥ ####### ⎢ ⎥ ####### ⎣ ⎦ Ta sẽ xét 3 trường hợp: ) số phương trình bằng số ẩn số nên ma trận [A] là ma trận vuông ) số phương trình nhỏ hơn số ẩn số ) số phương trình lớn hơn số ẩn số §2. NGHIỆM CỦA HỆ PHƯƠNG TRÌNH ĐẠI SỐ TUYẾN TÍNH
[ x ] = [ A ] − 1 [ b] (1)nếu ma trận A không suy biến, nghĩa là định thức của ma trận khác không. Các lệnh MATLAB để giải hệ là (ctsys): clc A = [1 2;3 4]; b = [‐1;‐1]; x = A^‐1*b %x = inv(A)*b
136 phương trình m ít hơn số ẩn số n thì nghiệm không duy nhất. Giả sử m hàng của ma trận hệ số [A] là độc lập thì vec tơ n chiều có thể phân tích thành hai thành phần: [ x ] = [ x ] + + [ x]− (2)Trong đó một ma trận là ma trận không gian hàng của ma trận [A] và được viết dưới dạng tổ hợp của: [ x ] + = [ A] [T α ] (3)và ma trận kia là ma trận không gian không sao cho: [ A ][ x] − = 0 (4)Như vậy: [ A ] [( x ] + + [ x] − ) = [ A ][ A] [ T α +] [ A][ x ]− = [ A ][ A] [ T α =] [ b] (5)Do [A][A]T là ma trận không suy biến m × m có được bằng cách nhân ma trận m × n với ma trận n × m nên ta có thể giải phương trình đối với [α] để có: [ α] 0 = ⎣ ⎡ AA T ⎤⎦ − 1 [ b] (6)Thay (6) vào (3) ta có: [ α] 0 + = [ A ] [T α ] 0 = [ A] T ⎡⎣ AA T⎤⎦ − 1 [ b] (7)Điều này thoả mãn phương trình [A][x] = [b]. Tuy nhiên nó không là nghiệm duy nhất vì nếu thêm bất kì một vec tơ [x] thoả mãn (4) thì nó sẽ cũng là nghiệm. MATLAB dùng lệnh pinv để giải hệ (ctpinv) A = [1 2]; b = 3; x = pinv(A)*b
[ ]e = [ A ][ x ] − [ b] (8)Vậy thì bài tiám của ta là cực tiểu hoá hàm: J = 0 e 2 = 0 [ A][ x] − [ b] 2 = 0 ⎡⎣[ A ][ x] − [ b] ⎤⎦ T⎡⎣[ A ][ x ] −[ b]⎤⎦ (9)Ta tìm cực tiểu của J bằng cách cho đạo hàm theo x của (9) bằng không. [ ] [ ][ ] [ ] [ ] [ ] [ ] [ ] [ ]T 0 T 1 T J A A x b 0 x A A A b x ####### ∂ ⎡ ⎤ − ####### = ⎡⎣ − ⎤⎦ = =⎣ ⎦ ####### ∂ ####### (10) 138
a x a x a x b a x a x a x b a x a x a x b ####### ⎧ + + = ####### ⎪ + + = ####### ⎨ ####### ⎪⎩ + + = ####### (1) ) Trước hết ta khử x 1 ra khỏi các phương trình, ngoại trừ phương trình đầu tiên, bằng cách nhân phương trình đầu tiên với ai1/a 11 (i là chỉ số hàng) và trừ đi mỗi phương trình đó: (0) 11 1 (0) 12 2 (0) 13 3 (0) 1 (1) 22 2 (1) 23 3 (1) 2 (1) 32 2 (1) 33 3 (1) 3 a x a x a x b a x a x b a x a x b ####### ⎧ + + = ####### ⎪ + = ####### ⎨ ####### ⎪ + = ####### ⎩ ####### (2) Trong đó: a (0)ij = aij b(0)i = bi với i = 1, j = 1, 2, 3 (1) (0) (0)i1 (0) ij ij (0) 1j 11 a a a a = −a (1) (0) (0)i1 (0) i i (0) 1 11 b b a b = − a với i, j = 2, 3 Việc này gọi là lấy trụ tại a 11 và phần tử a 11 gọi là trụ. ) Tiếp theo ta khử x 2 trong phương trình thứ 3 của (2) bằng cách lấy phương trình thứ 2 nhân với a (1)i2 / a(1) 22 (i = 3) và trừ đi phương trình thứ 3: (0) 11 1 (0) 12 2 (0) 13 3 (0) 1 (1) 22 2 (1) 23 3 (1) 2 (2) 33 3 (2) 3 a x a x a x b a x a x b a x b ####### ⎧ + + = ####### ⎪ + = ####### ⎨ ####### ⎪ = ####### ⎩ ####### (3) Trong đó: (2) (1) (1)i2 (1) ij ij (1) 2 j 22 a a a a a ####### = − (2) (1) (1)i2 (1) i i (1) 2 22 b b a b a \= − với i, j = 3 (4) Quá trình này được gọi là thuật toán khử Gauss tiến và được tổng quát hoá thành: (k) (k 1) (k 1)ik (k 1) ij ij (k 1) kj kk (k) (k 1) (k 1)ik (k 1) i i (k 1) k kk a a a a i, j k 1,k 2,...,m a b b a b i k 1,k 2,...,m a − − − − − − − − ####### = − = + + ####### = − = + + ####### (5) Để thực hiện thuật toán khử Gauss ta dùng đoạn mã lệnh: 139 for k = 1:n‐ 1 for i= k+1:n if A(i, k) ̃= 0 lambda = A(i, k)/A(k, k); A(i, k+1:n) = A(i, k+1:n) ‐ lambdaA(k, k+1:n); b(i)= b(i) ‐ lambdab(k); end end end Sau khi có hệ phương trình dạng ta giác ta tìm nghiệm dễ dàng. Từ phương trình thứ 3 của (3) ta có: (2) 3 3 (2) 33 x b a \= (6a) Thay vào phương trình thứ 2 ta có: (1) 2 (1) 23 3 2 (1) 22 x b a x a \= − (6b) và cuối cùng từ phương trình thứ nhất ta có: (0) 3 (0) 1 (0) 11 1 j 2 1j j x 1 b a x a = ####### ⎛ ⎞ ####### = ⎜ − ⎟ ####### ⎝ ⎠ ∑ (6c)Ta cũng có thể tổng quát hoá quá trình tìm nghiệm bằng cách tính lùi và tìm nghiệm bằng: (i 1) m (i 1) i (i 1) i j i 1 ij j ii x 1 b a x i m,m 1,..., a − − − = + ####### ⎛ ⎞ ####### = ⎜ − ⎟ = − ####### ⎝ ⎠ ∑ (7)và tìm nghiệm bằng đoạn mã lệnh: for k = n:‐1: b(k) = (b(k) ‐ A(k, k+1:n)*b(k+1:n))/A(k, k); end Như vậy phương pháp Gauss gồm hai bước: ‐ khử theo thuật toán Gauss ‐ tìm nghiệm của phương trình dạng tam giác Đoạn mã lệnh để tráo hàng được viết trong hàm swaprows(): 141 AB(m, k) = 0; end end %Tim nghiem x(NA, :) = AB(NA, NA+1:N); for m = NA‐1: ‐1: x(m, :) = AB(m, NA + 1:N)‐AB(m, m + 1:NA)*x(m + 1:NA, :); %(2.2) end Để giải hệ phương trình ta dùng ctgauss clear all, clc A = [ 1 1 1;2 ‐ 1 ‐1; 1 1 ‐ 1 ]; b = [ 2 0 1 ]ʹ; x = gauss(A, b)
function x = gaussjordan(A, B) %Kich thuoc cua ma tran A, B la NA x va NA x NB. %Ham nay dung giai he Ax = B bang thuat toan loai tru Gauss‐Jordan NA = size(A, 2); 142 [NB1,NB] = size(B); if NB1 = NA error(ʹA va B phai co kich thuoc tuong ungʹ); end for i = 1:NA if A(i, i) == 0 % trao hang neu can swaprows(A, i, mx); end c = A(i, i); for j = i:NA A(i,j) = A(i, j)/c; end B(i) = B(i)/c; for k = 1:NA if k=i c = A(k, i); A(k, i:NA) = A(k, i:NA)‐A(i, i:NA)*c; B(k) = B(k) ‐ B(i)*c; end end end x = B; và dùng chương trình ctgaussjordan giải hệ: clear all, clc a = [ 5 3 1;2 ‐ 1 1; 1 ‐ 1 ‐ 1 ]; b = [9; 2; ‐ 1 ]; x = gaussjordan(a, b) §4. GIẢI HỆ PHƯƠNG TRÌNH BẰNG CÁCH PHÂN TÍCH MA TRẬN
144 Áp dụng hàm doolittlesol() giải hệ phương trình: 1 2 3 4 3 6 x 1 8 3 10 x 0 4 12 10 x 0 ####### ⎡ − ⎤ ⎡ ⎤ ⎡ ⎤ ####### ⎢ − ⎥ ⎢ ⎥ =⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢⎣ − − ⎥ ⎢⎦ ⎣ ⎥⎦ ⎢ ⎥⎣ ⎦ ta dùng chương trình ctdoolitle: a = [ 4 ‐ 3 6; 8 ‐ 3 10; ‐ 4 12 ‐ 10 ]; b = [1; 0; 0 ]; x = doolittlesol(a, b)
function x = croutsol(a, b) %Ham dung giai he pt AX = B bang thuat toan Crout % Cu phap: x = croutsol(a, b) n =size(a,1); [l,r] = crout(a); y(1,:) = b(1)/l(1, 1); for m = 2:n y(m,:) = (b(m) ‐ l(m, 1:m‐1)*y(1:m‐1,:))/l(m, m); end x(n, :) = y(n)/r(n, n); for m = n‐1: ‐1: x(m, :) = (y(m) ‐ r(m, m + 1:n)*x(m + 1:n, :))/r(m, m); end Khi giải phương trình ta chương trình ctcrout: clear all, clc a = [ 4 8 20; 6 13 16; 20 16 ‐ 91 ]; b = [24; 18; ‐ 110 ]; 145 x = croutsol(a, b)
function x = choleskisol(a, b) %Giai he pt bang thuat toan Choleski %Cu phap: x = choleskisol(a, b) n =size(a,1); l = choleski(a); r = lʹ; y(1,:) = b(1)/l(1, 1); for m = 2:n y(m,:) = (b(m) ‐ l(m, 1:m‐1)*y(1:m‐1, :))/l(m, m); end x(n, :) = y(n)/r(n, n); for m = n‐1: ‐1: x(m, :) = (y(m) ‐r(m, m + 1:n)*x(m + 1:n, :))/r(m, m); end Để giải hệ phương trình 1 1 1 4 2 2 x 5 2 2 4 x 10 2 4 11 x 27 ####### ⎡ − ⎤ ⎡ ⎤ ⎡ ⎤ ####### ⎢ − − ⎥ ⎢ ⎥ = ⎢ − ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢⎣ − ⎥ ⎢⎦ ⎣ ⎥⎦ ⎢⎣ ⎥⎦ ta dùng chương trình ctcholeski: clear all, clc a = [ 4 ‐ 2 2;‐ 2 2 ‐4;2 ‐ 4 11 ]; b = [6; ‐10; 27 ]; x = choleskisol(a, b) 147 để giảm bớt số lượng phần tử cần lưu trữ. Bây giờ ta phân tích ma trận theo thuật toán Doolittle: hàng k ‐ (ck‐ 1 /dk‐ 1 )×hàng k‐ 1 → hàng k k = 1, 2,..., n và: dk ‐ (ck‐ 1 /dk‐ 1 )×ek‐ 1 → dk Để hoàn tất thuật việc phân tích, ta lưu hệ số λ = ck‐ 1 /dk‐ 1 vào vị trí của ck‐ 1 trước đó ck‐ 1 /dk‐ 1 → ck‐ 1 Như vậy thuật toán phân tích ma trận là: for k = 2:n lambda = c(k‐1)/d(k‐1); d(k) = d(k) ‐ lambda*e(k‐1) c(k‐1) = lambda; end Sau đó ta tìm nghiệm của phương trình [L][R][X] = [B] bằng cách giải phương trình [L][Y] = [B] và sau đó là phương trình [R][X] = [Y]. Phương trình [L][Y] = [B] có dạng: 1 1 1 2 2 2 3 3 3 4 4 n 1 n n 1 0 0 0 0 y b c 1 0 0 0 y b 0 c 1 0 0 y b 0 0 c 0 0 y b 0 0 0 0 c − 1 y b ####### ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ =⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ####### L ####### L ####### L ####### L ####### M M M M L M L L để tìm nghiệm [Y] bằng cách thay thế tiến ta dùng đoạn lệnh: y(1) = b(1); for k = 2:n y(k) = b(k) ‐ c(k‐1)*y(k‐1); end Phương trình [R][X] = [Y] có dạng: 148 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 n n n d e 0 0 0 x y 0 d e 0 0 x y 0 0 d e 0 x y 0 0 0 d 0 x y 0 0 0 0 0 d x y ####### ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ =⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ####### L ####### L ####### L ####### L ####### M M M M L M L L để tìm nghiệm [X] bằng cách thay thế lùi ta dùng đoạn lệnh: x(n) = y(n); for k = n‐1:‐1: x(k) = (y(k) ‐ e(k)*x(k+1))/d(k); end Ta xây dựng hàm band3() để phân tích ma trận dạng đường chéo: function [c, d, e] = band3(c , d, e) % Phan tich ma tran A = [c\d\e]. % Cu phap: [c, d, e] = band3(c, d, e) n = length(d); for k = 2:n lambda = c(k‐1)/d(k‐1); d(k) = d(k) ‐ lambda*e(k‐1); c(k‐1) = lambda; end Ta viết hàm band3sol() dùng để giải hệ phương trình có ma trận [A] dạng đường chéo: function x = band3sol(c ,d, e, b) % Giai he A*x = b voi A = [c\d\e] la tich LU % Cu phap: x =band3sol(c, d, e, b) [c, d, e] = band3(c, d, e); n = length(d); for k = 2:n % thay the tien b(k) = b(k) ‐ c(k‐1)*b(k‐1); 150 [ ] [ ] [ ]1 1 1 1 22 n 2 n 1 n 2 n 2 n n 1 d e d e f f d d e f e d e f d − − − − − ####### ⎡ ⎤ ####### ⎢ ⎥ ⎡ ⎤ ⎡ ⎤ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### = ⎢ ⎥ = ⎢ ⎥ =⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ####### ⎢ ⎥ ⎢⎣ ⎥⎦ ####### ⎣ ⎦ ####### M M ####### M Ta thực hiện thuật toán biến đổi ma trận: hàng (k + 1) ‐ (ek/dk) × hàng k → hàng (k + 1) hàng (k + 2) ‐ (fk/dk) × hàng k → hàng (k + 2) Các số hạng bị thay đổi theo thuật toán này là: dk+1 ‐ (ek/dk) ek → dk+ ek+1 ‐ (ek/dk) fk → ek+ dk+2 ‐ (fk/dk) fk → dk+ và lưu trữ lại: ek/dk → ek fk/dk → fk sau khi đã biến đổi ma trận, ta giải hệ phương trình có ma trận tam giác. Hàm band5() dùng để phân tích ma trận: function [d, e, f] = band5(d, e, f) % A = [f\e\d\e\f]. % Cu phap: [d, e, f] = band5(d, e, f) n = length(d); for k = 1:n‐ 2 lambda = e(k)/d(k); d(k+1) = d(k+1) ‐ lambdae(k); e(k+1) = e(k+1) ‐ lambdaf(k); e(k) = lambda; lambda = f(k)/d(k); d(k+2) = d(k+2) ‐ lambdaf(k); f(k) = lambda; end lambda = e(n‐1)/d(n‐1); d(n) = d(n) ‐ lambdae(n‐1); e(n‐1) = lambda; 151 Ta viết hàm band5sol() để giải hệ phương trình: function x = band5sol(d, e, f, b) % Giai he A*x = b voi A = [f\e\d\e\f] % Cu phap: x = band5sol(d, e, f, b) [e,d,f ] = band5(e, d, f); n = length(d); b(2) = b(2) ‐ e(1)*b(1); for k = 3:n b(k) = b(k) ‐ e(k‐1)*b(k‐1) ‐ f(k‐2)*b(k‐2); end Để giải hệ phương trình 1 2 3 4 5 6 1 1 2 0 0 0 x 4 1 2 3 1 0 0 x 7 2 3 3 2 2 0 x 12 0 1 2 1 2 1 x 7 0 0 2 2 2 1 x 5 0 0 0 1 1 1 x 1 ####### ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ =⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ − ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎢ − ⎥ ⎢ ⎥ ⎢ ⎥ ####### ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ta dùng chương trình cban5eq clear all, clc d = [1 2 3 1 2 1]ʹ; e = [1 3 2 2 ‐1]ʹ; f = [2 1 2 1]ʹ; b = [4 7 12 7 5 1]; x = band5sol(d, e, f, b) §6. CÁC PHƯƠNG PHÁP LẶP ĐỂ GIẢI HỆ PHƯƠNG TRÌNH ĐẠI SỐ TUYẾN TÍNH Nói chung có hai phương pháp giải hệ phương trình đại số tuyến tính: phương pháp trực tiếp và phương pháp lặp. Các bài toán kĩ thuật thường đưa về hệ phương trình đại số tuyến tính có ma trận [A] thưa và lớn nên các phương pháp lặp rất thích hợp. Các phương pháp lặp được chia thành hai loại: phương pháp lặp tĩnh và phương pháp lặp động. 153 với sai số trong số liệu. Nó cho biết độ chính xác của kết quả từ phép nghịch đảo ma trận và nghiệm của hệ phương trình đại số tuyến tính).
154 BiCG. Nó được dùng cho hệ phương trình có ma trận hệ số không đối xứng. - Phương pháp Chebyshev: Phương pháp này tính lặp các đa thức với các hệ số được chọn để cực tiểu hoá chuẩn của số dư theo nghĩa min ‐ max. Ma trận hệ số phải xác định dương. Nó được dùng cho hệ phương trình có ma trận hệ số không đối xứng. Ta biết rằng tốc độ hội tụ của phép lặp phụ thuộc rất nhiều vào phổ của ma trận(các giá trị riêng của ma trận). Do vậy phép lặp thường đưa thêm một ma trận thứ hai để biến đổi ma trận hệ số thành ma trận có phổ thích hợp. Ma trận biến đổi như vậy gọi là ma trận điều kiện trước(preconditioner). Một preconditioner tốt sẽ cải thiện sự hội tụ của phương pháp lặp. Nhiều trường hợp, nếu không có preconditioner, phép lặp sẽ không hội tụ. Preconditioner đơn giản nhất chính là ma trận đường chéo mà Mi,j = Ai,j nếu i = j và các phần tử khác bằng zero. Ma trận như vậy gọi là ma trận điều kiện trước Jacobi. Trong tính toán, tồn tại hai loại ma trận điều kiện trước: ‐ ma trận [M] xấp xỉ ma trận [A] và làm cho việc giải hệ [M][X] = [B] dễ hơn giải hệ [A][X] = [B] ‐ ma trận [M] xấp xỉ [A]‐ 1 sao cho chỉ cần tính [M][B] là có [X] Phần lớn các ma trận [M] thuộc loại thứ nhất. §7. PHƯƠNG PHÁP LẶP JACOBI Xét hệ phương trình AX = F. Bằng cách nào đó ta đưa hệ phương trình về dạng X = BX + G trong đó: B = (bij)n,n G = (g 1 ,g 2 ,...,gn)T Chọn vectơ: X = ( x1(o),x2(o),....,xn(o) )T làm xấp xỉ thứ 0 của nghiệm đúng và xây dựng xấp xỉ X(m+1) = BX(m) + G ( m = 0,1,....) Người ta chứng minh rằng nếu phương trình ban đầu có nghiệm duy nhất và một trong ba chuẩn của ma trận B nhỏ hơn 1 thì dãy xấp xỉ hội tụ về nghiệm duy nhất đó. Cho một ma trận B, chuẩn của ma trận B, kí hiệu B là một trong 3 số : |
