Giáo trình môn Matlab

541 trang hoanguyen 6100

Download

Bạn đang xem 20 trang mẫu của tài liệu "Giáo trình môn Matlab", để 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:

giao_trinh_mon_matlab.pdf

Nội dung text: Giáo trình môn Matlab

CHƯƠNG 1: MATLAB CƠ BẢN §1. CÁC TOÁN TỬ CƠ BẢN CỦA MATLAB 1. Các toán tử cơ bản: Matlab là một phần mềm cao cấp dùng để giải các bài toán. Để khởi động MATLAB ta bấm đúp vào icon của nó. Các file MATLAB có dạng *.m và chỉ chạy trong môi trường MATLAB. MATLAB xử lí số liệu như là ma trận. Khi ta đánh lệnh vào cửa sổ lệnh, nó sẽ được thi hành ngay và kết quả hiện lên màn hình. Nếu ta không muốn cho kết quả hiện lên màn hình thì sau lệnh ta đặt thêm dấu “;”. Nếu lệnh quá dài, không vừa một dòng dòng có thể đánh lệnh trên nhiều dòng và cuối mỗi dòng đặt thêm dấu rồi xuống dòng. Khi soạn thảo lệnh ta có thể dùng các phím tắt : ↑ Ctrl‐P gọi lại lệnh trước đó ↓ Ctrl‐N gọi lệnh sau ← Ctrl‐B lùi lại một kí tự → Ctrl‐F tiến lên một kí tự Ctrl‐→ Ctrl‐R sang phải một từ Ctrl‐← Crtl‐L sang phải một từ home Ctrl‐A về đầu dòng end Ctrl‐E về cuối dòng esc Ctrl‐U xoá dòng del Ctrl‐D xoá kí tự tại chỗ con nháy đứng backspace Ctrl‐H xoá kí tự trước chỗ con nháy đứng ) Các phép toán cơ bản của MATLAB gồm: + cộng ‐ trừ * nhân / chia phải \ chia trái ^ luỹ thừa ‘ chuyển vị ma trận hay số phức liên hợp ) Các toán tử quan hệ : lớn hơn >= lớn hơn hoặc bằng == bằng 1
~= không bằng ) Các toán tử logic : & và | or ~ not ) Các hằng : pi 3.14159265 i số ảo j tương tự i eps sai số 2‐52 realmin số thực nhỏ nhất 2‐1022 realmax số thực lớn nhất 21023 inf vô cùng lớn NaN Not a number 2. Nhập xuất dữ liệu từ dòng lệnh: MATLAB không đòi hỏi phải khai báo biến trước khi dùng. MATLAB phân biệt chữ hoa và chữ thường. Các số liệu đưa vào môi trường làm việc của MATLAB được lưu lại suốt phiên làm việc cho đến khi gặp lệnh clear all. MATLAB cho phép ta nhập số liệu từ dòng lệnh. Khi nhập ma trận từ bàn phím ta phải tuân theo các quy định sau : • ngăn cách các phần tử của ma trận bằng dấu “,” hay dấu trống • dùng dấu “;” để kết thúc một hàng • bao các phần tử của ma trận bằng cặp dấu ngoặc vuông [ ] Để nhập các ma trận sau: ⎡⎤124 ⎡⎤ 1 ⎢⎥ ⎢⎥ A325=−⎢⎥ B1421 =⎣⎦⎡⎤ − C4 = ⎢⎥ ⎣⎦⎢⎥153 ⎣⎦⎢⎥ 7 ta dùng các lệnh: A = [ 1 2 3; 3 ‐2 4; 1 5 3] B = [ 1 4 2 1] C = [ 1; 4; 7] 3. Nhập xuất dữ liệu từ file: MATLAB có thể xử lí hai kiểu file dữ liệu: file 2
nhị phân *.mat và file ASCII *.dat. Để lưu các ma trận A, B, C dưới dạng file nhị phân ta dùng lệnh: save ABC A B C và nạp lại các ma trận A, B bằng lệnh: load ABC A B Nếu muốn lưu số liệu của ma trận B dưới dạng file ASCII ta viết: save b.dat B /ascii Ta viết chương trình ct1_1.m như sau: clear A = [1 2 3; 4 5 6] B = [3; ‐2; 1]; C(2) = 2; C(4) = 4 disp(’Nhan phim bat ky de xem nhap/xuat du lieu tu file’) save ABC A B C %luu A,B & C duoi dang MAT‐file co ten ’ABC.mat’ clear(’A’, ’C’) %xoa A va C khoi bo nho load ABC A C %doc MAT ‐ file de nhap A va C vao bo nho save b.dat B /ascii %luu B duoi dang file ASCII co ten ’b.dat’ clear B load b.dat %doc ASCII b x = input(’Nhap x:’) format short e x format rat, x format long, x format short, x 4. Nhập xuất dữ liệu từ bàn phím: Lệnh input cho phép ta nhập số liệu từ bàn phím. Ví dụ: 3
x = input(’Nhap x: ’) Lệnh format cho phép xác định dạng thức của dữ liệu. Ví dụ: format rat % so huu ti format long % so sẽ có 14 chu so sau dau phay format long e % so dang mu format hex % so dang hex format short e %so dang mu ngan format short %tro ve so dang ngan (default) Một cách khác để hiển thị giá trị của biến và chuỗi là đánh tên biến vào cửa số lệnh MATLAB. Ta cũng có thể dùng disp và fprintf để hiển thị các biến. Ví dụ: disp(ʹTri so cua x = ʹ), disp(x) Ta viết chương trình ct1_2.m như sau: clc f = input(ʹNhap nhiet do Fahrenheit[F]:ʹ); c = 5/9*(f ‐ 32); fprintf(ʹ%5.2f(do Fahrenheit) la %5.2f(do C).\nʹ, f, c) fid = fopen(ʹct1_2.datʹ, ʹwʹ); fprintf(fid, ʹ%5.2f(do Fahrenheit) la %5.2f(do C).\nʹ, f, c); fclose(fid); Trong trường hợp ta muốn nhập một chuỗi từ bàn phím, ta cần phải thêm kí tự s vào đối số. Ví dụ: ans = input(ʹBan tra loi hoac : ʹ,ʹsʹ) 5. Các hàm toán học: a. Các hàm toán học cơ bản: x exp(x) hàm e sqrt(x) căn bậc hai của x log(x) logarit tự nhiên 4
log10(x) logarit cơ số 10 abs(x) modun của số phức x angle(x) argument của số phức a conj(x) số phức liên hợp của x imag(x) phần ảo của x real(x) phần thực của x sign(x) dấu của x cos(x) sin(x) tan(x) acos(x) asin(x) atan(x) cosh(x) coth(x) sinh(x) tanh(x) acosh(x) acoth(x) asinh(x) atanh(x) b. Các hàm toán học tự tạo: MATLAB cho phép ta tạo hàm toán học và lưu nó vào một file để dùng như là hàm có sẵn của MATLAB. Ví dụ ta cần tạo hàm: 1 f(x)= 1 18x+ 2 và hàm: 22 ⎡⎤f(x,x)112 ⎡⎤x4x512+− f(x)2 ==⎢⎥⎢⎥2 ⎣⎦f(x,x)212 ⎣⎦2x112−−− 2x 3x 2.5 Muốn thế ta tạo ra file f1.m như sau: function y = f1(x) y = 1./(1+8*x.^2); và file f2.m: 5
function y = f2(x) y(1) = x(1)*x(1)+4*x(2)*x(2) ‐5; y(2) = 2*x(1)*x(1)-2*x(1)-3*x(2) -2.5; Khi nhập lệnh f1(2) ta có giá trị của hàm f1 tại x = 2. Khi nhập lệnh f2([2 4]) ta có giá trị của hàm f2 tại x1 = 2 và x2 = 4. Lệnh feval(‘f1’, 2) và feval(‘f2’, [2 4]) cũng cho kết quả tương tự. Cách thứ hai để biểu diễn một hàm toán học một biến trên dòng lệnh là tạo ra một đối tượng inline từ một biểu thức chuỗi. Ví dụ ta có thể nhập từ dòng lệnh hàm như sau: f1 = inline(’1./(1 + 8*x.^2)’,’x’); f1([0 1]), feval(f1, [0 1]) Ta cũng có thể viết: f1 = ʹ1./(1 + 8*x.^2)ʹ; x = [0 1]; eval(f1) Nếu hàm là đa thức ta chỉ cần nhập ma trận các hệ số từ số mũ cao nhất. Ví dụ với đa thức P4(x) = x4 + 4x3 + 2x + 1 ta viết: P = [1 4 0 2 1] Để nhân hai đa thức ta dùng lệnh conv; để chia 2 đa thức ta dùng lệnh deconv. Muốn tính trị số của đa thức ta dùng lệnh polyval và lệnh polyvalm dùng khi đa thức là ma trận. c. Các lệnh xử lí hàm: Lệnh fplot vẽ đồ thị hàm toán học giữa các giá trị đã cho. Ví dụ: fplot(‘f1’, [‐5 5 ]) grid on Cho một hàm toán học một biến, ta có thể dùng lệnh fminbnd của MATLAB để tìm cực tiểu địa phương của hàm trong khoảng đã cho. Ví dụ: 6
f = inline(ʹ1./((x ‐ 0.3).^2+0.01) + 1./((x ‐ 0.9).^2 + 0.04) ‐ 6 ʹ); x = fminbnd(f, 0.3, 1) Lệnh fminsearch tương tự hàm fminbnd dùng để tìm cực tiểu địa phương của hàm nhiều biến. Ta có hàm 3 biến lưu trong file three_var.m như sau: function b = three_var(v) x = v(1); y = v(2); z = v(3); b = x.^2 + 2.5*sin(y) ‐ z^2*x^2*y^2; Bây giờ tìm cực tiểu đối với hàm này bắt đầu từ x = ‐0.6 , y = ‐1.2 và z = 0.135 bằng các lệnh: v = [‐0.6 ‐1.2 0.135]; a = fminsearch(ʹthree_varʹ, v) Lệnh fzero dùng để tìm điểm zero của hàm một biến. Ví dụ để tìm giá trị không của hàm lân cận giá trị ‐0.2 ta viết: f = inline(ʹ1./((x ‐ 0.3).^2 + 0.01) + 1./((x ‐ 0.9).^2 + 0.04) ‐ 6ʹ); a = fzero(f, ‐0.2) Zero found in the interval: [‐0.10949, ‐0.264]. a = ‐0.1316 6. Các phép toán trên ma trận và vec tơ: a. Khái niệm chung: Giả sử ta tạo ra các ma trận a và b bằng các lệnh: a = [1 2 3; 4 5 6]; b = [3 ‐2 1]; Ta có thể sửa đổi chúng: 7
A = [a; 7 8 9] B = [b; [1 0 ‐1]]ʹ Toán tử ‘ dùng để chuyển vị một ma trận thực và chuyển vị liên hợp một ma trận phức. Nếu chỉ muốn chuyển vị ma trận phức, ta dùng thêm toán tử “.” nghĩa là phải viết “.’”. Ví dụ: C = [1 + 2*i 2 ‐ 4*i; 3 + i 2 ‐ 2*j]; X = Cʹ Y = C.’ b. Chỉ số: Phần tử ở hàng i cột j của ma trận m×n có kí hiệu là A(i, j). Tuy nhiên ta cũng có thể tham chiếu tới phần tử của mảng nhờ một chỉ số, ví dụ A(k) với k = i + (j ‐ 1)m. Cách này thường dùng để tham chiếu vec tơ hàng hay cột. Trong trường hợp ma trận đầy đủ thì nó được xem là ma trận một cột dài tạo từ các cột của ma trận ban đầu. Như vậy viết A(5) có nghĩa là tham chiếu phần tử A(2, 2). Để xác định kích thước của một ma trận ta dùng lệnh length(trả về kích thước lớn nhất) hay size(số hàng và cột). Ví dụ: c = [1 2 3 4; 5 6 7 8]; length(c) [m, n] = size(c) c. Toán tử “:” : Toán tử “:” là một toán tử quan trọng của MATLAB. Nó xuất hiện ở nhiều dạng khác nhau. Ví dụ: 1:10 tạo một vec tơ hàng chứa 10 số nguyên từ 1 đến 10. Lệnh: 100: ‐7: 50 tạo một dãy số từ 100 đến 51, giảm 7 mỗi lần. Lệnh: 0: pi/4: pi 8
tạo một dãy số từ 0 đến pi, cách đều nhau pi/4 Các biểu thức chỉ số tham chiếu tới một phần của ma trận. Viết A(1:k, j) là tham chiếu đến k phần tử đầu tiên của cột j. Ngoài ra toán tử “:” tham chiếu tới tất cả các phần tử của một hàng hay một cột. Ví dụ: B = A(:, [1 3 2 ]) tạo ra ma trận B từ ma trận A bằng cách đổi thứ tự các cột từ [1 2 3] thành [1 3 2] d. Tạo ma trận bằng hàm có sẵn: MATLAB cung cấp một số hàm để tạo các ma trận cơ bản: zeros tạo ra ma trận mà các phần tử đều là zeros z = zeros(2, 4) ones tạo ra ma trận mà các phần tử đều là 1 x = ones(2, 3) y = 5*ones(2, 2) rand tạo ra ma trận mà các phần tử ngẫu nhiên phân bố đều d = rand(4, 4) randn tạo ra ma trận mà các phần tử ngẫu nhiên phân bố trực giao e = randn(3, 3) magic(n) tạo ra ma trận cấp n gồm các số nguyên từ 1 đến n2 với tổng các hàng bằng tổng các cột n phải lớn hơn hay bằng 3. pascal(n) tạo ra ma trận xác định dương mà các phần tử lấy từ tam giác Pascal. pascal(4) eye(n) tạo ma trận đơn vị 9
eye(3) eye(m, n) tạo ma trận đơn vị mở rộng eye(3, 4) e. Lắp ghép: Ta có thể lắp ghép(concatenation) các ma trận có sẵn thành một ma trận mới. Ví dụ: a = ones(3, 3) b = 5*ones(3, 3) c = [a + 2; b] f. Xoá hàng và cột : Ta có thể xoá hàng và cột từ ma trận bằng dùng dấu []. Để xoá cột thứ 2 của ma trận b ta viết: b(:, 2) = [] Viết x(1: 2: 5) = [] nghĩa là ta xoá các phần tử bắt đầu từ đến phần tử thứ 5 và cách 2 rồi sắp xếp lại ma trận. g. Các lệnh xử lí ma trận: Cộng : X= A + B Trừ : X= A ‐ B Nhân : X= A * B : X.*A nhân các phần tử tương ứng với nhau Chia : X = A/B lúc đó X*B = A : X = A\B lúc đó A*X = B : X=A./B chia các phần tử tương ứng với nhau Luỹ thừa : X = A^2 : X = A.^2 Nghịch đảo : X = inv(A) Định thức : d = det(A) 7. Tạo số ngẫu nhiên: MATLAB có các lệnh tạo số ngẫu nhiên là rand và randn tạo ra các số ngẫu nhiên theo phân bố Gauss. rand(m, n) tạo ra ma trận các số ngẫu nhiên phân bố đồng nhất. randn(m, n) tạo ra ma trận các số ngẫu nhiên theo phân bố chuẩn Gauss. rand(3, 3) 10
randn(3, 3) 8. Các lệnh dùng lập trình: a. Các phát biểu điều kiện if, else, elseif: Cú pháp của if: if end Nếu cho kết quả đúng thì phần lệnh trong thân của if được thực hiện. Các phát biểu else và leseif cũng tương tự. Ví dụ: Ta xét chương trình) ct1_4. m để đoán tuổi như sau: clc disp(‘Xin chao! Han hanh duoc lam quen’); x = fix(30*rand); disp(‘Tuoi toi trong khoang 0 ‐ 30’); gu = input(‘Xin nhap tuoi cua ban: ‘); if gu x disp(‘Ban lon hon toi’); else disp(‘Ban bang tuoi toi’); end b. switch: Cú pháp của switch như sau : switch case n1 : case n2 : . . . . . . . . . . . . . . . case nn : otherwise : end c. while: vòng lặp while dùng khi không biết trước số lần lặp. Cú pháp của nó như sau: 11
while end Xét chương trình in ra chuoi “Xin chao” lên mà hình với số lần nhập từ bàn phím ct1_5.m như sau: clc disp(ʹxin chaoʹ); gu = input(ʹNhap so lan in: ʹ); i = 0; while i ~= gu disp([ʹXin chaoʹ i]); i = i + 1 end d. for: vòng lặp for dùng khi biết trước số lần lặp. Cú pháp như sau: for = : : Ta xây dựng chương trình đoán số ct1_6.m: clc x = fix(100*rand); n = 7; t = 1; for k = 1:7 num = int2str(n); disp([ʹBan co quyen du doan ʹ, num, ʹ lanʹ]); disp(ʹSo can doan nam trong khoang 0 ‐ 100ʹ); gu = input(ʹNhap so ma ban doan: ʹ); if gu x disp(ʹSo ban doan lon honʹ); else disp(ʹBan da doan dung. Xin chuc mungʹ); t = 0; break; end 12
n = n ‐ 1; end if t > 0 disp(ʹBan khong doan ra roiʹ); numx = int2str(x); disp([ʹDo la so: ʹ, numx]); end e. break: phát biểu break để kết thúc vòng lặp for hay while mà không quan tâm đến điều kiện kết thúc vòng lặp đã thoả mãn hay chưa. §2. ĐỒ HOẠ TRONG MATLAB 1. Các lệnh vẽ: MATLAB cung cấp một loạt hàm để vẽ biểu diễn các vec tơ số liệu cũng như giải thích và in các đường cong này. plot đồ họa 2‐D với số liệu 2 trục vô hướng và tuyến tính plot3 đồ họa 3‐D với số liệu 2 trục vô hướng và tuyến tính polar đồ hoạ trong hệ toạ độ cực loglog đồ hoạ với các trục logarit semilogx đồ hoạ với trục x logarit và trục y tuyến tính semilogy đồ hoạ với trục y logarit và trục x tuyến tính plotyy đồ hoạ với trục y có nhãn ở bên trái và bên phải 2. Tạo hình vẽ: Hàm plot có các dạng khác nhau phụ thuộc vào các đối số đưa vào. Ví dụ nếu y là một vec tơ thì plot(y) tạo ra một đường thẳng quan hệ giữa các giá trị của y và chỉ số của nó. Nếu ta có 2 vec tơ x và y thì plot(x, y) tạo ra đồ thị quan hệ giữa x và y. t = [0: pi/100: 2*pi] y = sin(t); plot(t, y) grid on polar(t, y) 3. Đặc tả kiểu đường vẽ: Ta có thể dùng các kiểu đường vẽ khác nhau khi vẽ hình. Muốn thế ta chuyển kiểu đường vẽ cho hàm plot. Ta viết chương trình ct1_7.m tạo ra đồ thị hàm hình sin: 13
t = [0: pi/100: 2*pi]; y = sin(t); plot(t, y, ’. ‘) % vẽ bằng đường chấm chấm grid on 4. Đặc tả màu và kích thước đường vẽ: Để đặc tả màu và kích thước đường vẽ ta dùng các tham số sau: LineWidth độ rộng đường thẳng,tính bằng số điểm MarkerEdgeColor màu của các cạnh của khối đánh dấu MarkerFaceColor màu của khối đánh dấu MarkerSize kích thước của khối đánh dấu Màu được xác định bằng các tham số: Mã Màu Mã Màu r red m magenta g green y yellow b blue k black c cyan w white Các dạng điểm đánh dấu xác định bằng: Mã Kiểu đánh dấu Mã Kiểu đánh dấu + dấu cộng . điểm o vòng tròn x chữ thập * dấu sao s hình vuông d hạt kim cương v điểm tam giác hướng xuống ^ điểm tam giác hướng lên tam giác sang phải h lục giác p ngũ giác Các dạng đường thẳng xác định bằng: Mã Kiểu đường Mã Kiểu đường ‐ đường liền : đường chấm chấm ‐‐ đường đứt nét ‐. đường chấm gạch 14
Ta xét chương trình ct1_8.m như sau: x = ‐pi : pi/10 : pi; y = tan(sin(x)) ‐ sin(tan(x)); plot(x, y, ʹ‐‐rs’, ʹLineWidthʹ, 2, ʹMarkerEdgeColorʹ, ʹkʹ, ʹMarkerFaceColorʹ, ʹgʹ, ʹMarkerSizeʹ, 10) Chương trình này sẽ vẽ đường cong y = f(x) có các đặc tả sau : ‐ đường vẽ là đường đứt nét(‐‐) ‐ khối đánh dấu hình vuông (s), đường vẽ màu đỏ(r) ‐ đường vẽ rộng 2 point ‐ các cạnh của khối đánh màu đen ‐ khối đánh dấu màu green ‐ kích thước khối đánh dấu 10 point 5. Thêm đường vẽ vào đồ thị đã có: Để làm điều này ta dùng lệnh hold. Khi ta đánh lệnh hold on thì MATLAB không xoá đồ thị đang có. Nó thêm số liệu vào đồ thị mới này. Nếu phạm vi giá trị của đồ thị mới vượt quá các giá trị của trục toạ độ cũ thì nó sẽ định lại tỉ lệ xích. 6. Chỉ vẽ các điểm số liệu: Để vẽ các điểm đánh dấu mà không nối chúng lại với nhau ta dùng đặc tả nói rằng không có các đường nối giữa các điểm, nghĩa là ta gọi hàm plot chỉ với đặc tả màu và điểm đánh dấu. Ta xét chương trình ct1_9.m như sau: x = ‐pi : pi/10 : pi; y = tan(sin(x)) ‐ sin(tan(x)); plot(x, y, ʹsʹ, ʹMarkerEdgeColorʹ, ʹkʹ) 7. Vẽ các điểm và đường: Để vẽ cả các điểm đánh dấu và đường nối giữa chúng ta cần mô tả kiểu đường và kiểu điểm. Ta xét chương trình ct1_10.m: x = 0:pi/15:4*pi; y = exp(2*sin(x)); plot(x, y, ʹ‐rʹ, x, y, ʹokʹ) dùng vẽ đường cong y = f(x) có đường nối liền, màu đỏ. Điểm đánh dấu là 15
chữ o có màu đen. 8. Vẽ với hai trục y: Lệnh plotyy cho phép tạo một đồ thị có hai trục y. Ta cũng có thể dùng plotyy để cho giá trị trên hai trục y có kiểu khác nhau nhằm tiện so sánh. Ta xét chương trình ct1_11.m: t = 0:900; A = 1000; b = 0.005; a = 0.005; z2 = sin(b*t); z1 = A*exp(‐a*t); [haxes, hline1, hline2] = plotyy(t, z1, t, z2,ʹsemilogyʹ, ʹplotʹ); 9. Vẽ đường cong với số liệu 3 ‐ D: Nếu x, y, z là 3 vec tơ có cùng độ dài thì plot3 sẽ vẽ đường cong 3D. Ta viết chương trình ct1_12.m: t = 0:pi/50:10*pi; plot3(sin(t),cos(t),t) axis square; grid on 10. Đặt các thông số cho trục: Khi ta tạo một hình vẽ, MATLAB tự động chọn các giới hạn trên trục toạ độ và khoảng cách đánh dấu dựa trên số liệu dùng để vẽ. Tuy nhiên ta có thể mô tả lại phạm vi giá trị trên trục và khoảng cách đánh dấu theo ý riêng. Ta có thể dùng các lệnh sau: axis đặt lại các giá trị trên trục toạ độ axes tạo một trục toạ độ mới với các đặc tính được mô tả get và set cho phép xác định và đặt các thuộc tính của trục toạ độ đang có gca trở về trục toạ độ cũ MATLAB chọn các giới hạn trên trục toạ độ và khoảng cách đánh dấu dựa trên số liệu dùng để vẽ. Dùng lệnh axis có thể đặt lại giới hạn này. Cú pháp của lệnh: axis[ xmin , xmax , ymin , ymax] Ta xét chương trình ct1_13.m như sau: 16
x = 0:0.025:pi/2; plot(x, tan(x), ʹ‐roʹ) axis([0 pi/2 0 5]) MATLAB chia vạch trên trục dựa trên phạm vi dữ liệu và chia đều. Ta có thể mô tả cách chia nhờ thông số xtick và ytick bằng một vec tơ tăng dần. Ví dụ xét chương trình ct1_14.m: x = ‐pi: .1: pi; y = sin(x); plot(x, y) set(gca, ʹxtickʹ, ‐pi :pi/2:p); set(gca, ʹxticklabelʹ, {ʹ‐piʹ, ʹ‐pi/2ʹ, ʹ0ʹ, ʹpi/2ʹ, ʹpiʹ}) 11. Ghi nhãn lên các trục toạ độ: MATLAB cung cấp các lệnh ghi nhãn lên đồ hoạ gồm : title thêm nhãn vào đồ hoạ xlabel thêm nhãn vào trục x ylabel thêm nhãn vào trục y zlabel thêm nhãn vào trục z legend thêm chú giải vào đồ thị text hiển thị chuỗi văn bản ở vị trí nhất định gtext đặt văn bản lên đồ hoạ nhờ chuột \bf bold font \it italics font \sl oblique font (chữ nghiêng) \rm normal font Các kí tự đặc biệt xem trong String properties của Help. Ta dùng các lệnh xlabel , ylabel , zlabel để thêm nhãn vào các trục toạ độ. Ta có thể thêm văn bản vào bất kì chỗ nào trên hình vẽ nhờ hàm text. Ta có chương trình ct1_15.m: x = ‐pi: .1: pi; y = sin(x); plot(x, y) xlabel(ʹt = 0 to 2\piʹ, ʹFontsizeʹ, 16) ylabel(ʹsin(t)ʹ, ʹFontsizeʹ, 16) 17
title(ʹ\it{Gia tri cua sin tu zero đến 2 pi}ʹ, ʹFontsizeʹ, 16) text(3*pi/4, sin(3*pi/4),ʹ\leftarrowsin(t ) = 0.707ʹ, ʹFontSizeʹ, 12) 12. Định vị văn bản trên hình vẽ: Ta có thể sử dụng đối tượng văn bản để ghi chú các trục ở vị trí bất kì. MATLAB định vị văn bản theo đơn vị dữ liệu trên trục. Ví dụ để vẽ hàm y = Aeαt với A = 0.25 , t = 0 đến 900 và α = 0.005 ta viết chương trình ct1_16.m: t = 0: 900; plot(t, 0.25*exp(‐0.005*t)) plot(t, y) text(300, .25*exp(‐.005*300), ʹ\bullet\leftarrow\fontname{times}0.25{\ite}^{‐ 0.005{\itt}} tai, {\itt} = 300ʹ, ʹFontSizeʹ, 14)%ghi chu tai t = 300 Tham số HorizontalAlignment và VerticalAlignment định vị văn bản so với các toạ độ x, y, z đã cho. 13. Đồ hoạ đặc biệt: a. Khối và vùng: Đồ hoạ khối và vùng biểu diễn số liệu là vec tơ hay ma trận. MATLAB cung cấp các hàm đồ hoạ khối và vùng : bar hiển thị các cột của ma trận m*n như là m nhóm, mỗi nhóm có n bar barh hiển thị các cột của ma trận m*n như là m nhóm, mỗi nhóm có n bar nằm ngang bar3 hiển thị các cột của ma trận m*n như là m nhóm, mỗi nhóm có n bar dạng 3D bar3h hiển thị các cột của ma trận m*n như là m nhóm, mỗi nhóm có n bar dạng 3D nằm ngang Mặc định, mỗi phần tử của ma trận được biểu diễn bằng một bar. Ta xét chương trình ct1_17.m: y = [5 2 1 6 7 3 8 6 3 5 5 5 1 5 8]; 18
bar(y) b. Mô tả dữ liệu trên trục: Ta dùng các hàm xlabel và ylabel để mô tả các dữ liệu trên trục. Ta xét chương trình ct1_18.m: nhdo = [29 23 27 25 20 23 23 27]; ngay = 0: 5: 35; bar(ngay, nhdo) xlabel(ʹNgayʹ) ylabel(ʹNhiet do (^{o}C)ʹ) set(gca,ʹYLimʹ,[15 30],ʹLayerʹ,ʹtopʹ) grid on set(gca,ʹYLimʹ,[15 30]) Mặc định,phạm vi giá trị của trục y là từ 0 đến 30. Để xem nhiệt độ trong khoảng từ 15 đến 30 ta thay đổi phạm vi giá trị của trục y: set(gca,ʹYLimʹ,[15 30],ʹLayerʹ,ʹtopʹ) và trên đồ thị, phạm vi giá trị của trục y đã thay đổi. c. Xếp chồng đồ thị: Ta có thể xếp chồng số liệu trên đồ thị thanh bằng cách tạo ra một trục khác trên cùng một vị trí và như vậy ta có một trục y độc lập với bộ số liệu khác. TCE = [515 420 370 250 135 120 60 20]; nhdo = [29 23 27 25 20 23 23 27]; ngay = 0:5:35; bar(ngay, nhdo) xlabel(ʹNgayʹ) ylabel(ʹNhiet do (^{o}C)ʹ) Để xếp chồng một số liệu lên một đồ thị thanh ở trên, có trục thứ 2 ở cùng vị trí như trục thứ nhất ta viết: h1 = gca; và tạo trục thứ 2 ở vị trí trục thứ nhất trước nhất vẽ bộ số liệu thứ 2: h2 = axes(ʹPositionʹ,get(h1,ʹPositionʹ)); 19
plot(days,TCE,ʹLineWidthʹ,3) Để trục thứ 2 không gây trở ngại cho trục thứ nhất ta viết: set(h2,ʹYAxisLocationʹ,ʹrightʹ,ʹColorʹ,ʹnoneʹ,ʹXTickLabelʹ,[]) set(h2,ʹXLimʹ,get(h1,ʹXLimʹ),ʹLayerʹ,ʹtopʹ) Để ghi chú lên đồ thị ta viết: text(11,380,ʹMat doʹ,ʹRotationʹ,‐‐55,ʹFontSizeʹ,16) ylabel(ʹTCE Mat do (PPM)ʹ) title(ʹXep chong do thiʹ,ʹFontSizeʹ,16) (lưu trong ct1_19.m) d. Đồ hoạ vùng: Hàm area hiển thị đường cong tạo từ một vec tơ hay từ một cột của ma trận. Nó vẽ các giá trị của một cột của ma trận thành một đường cong riêng và tô đầy vùng không gian giữa các đường cong và trục x. ta xét chương trình ct1_20.m: Y = [5 1 2 8 3 7 9 6 8 5 5 5 4 2 3]; area(Y) hiển thị đồ thị có 3 vùng, mỗi vùng một cột. Độ cao của mỗi đồ thị vùng là tổng các phần tử trong một hàng. Mỗi đường cong sau sử dụng đường cong trước làm cơ sở. Để hiển thị đường chia lưới ta dùng lệnh: set(gca,ʹLayerʹ,ʹtopʹ) set(gca,ʹXTickʹ,1:5) grid on f. Đồ thị pie: Đồ thị pie hiển thị theo tỉ lệ phần trăm của một phần tử của một vec tơ hay một ma trận so với tổng các phần tử. Các lệnh pie và pie3 tạo ra đồ thị 2D và 3D. ta xét chương trình ct1_21.m: X = [19.3 22.1 51.6; 34.2 70.3 82.4; 61.4 82.9 90.8; 20
50.5 54.9 59.1; 29.4 36.3 47.0]; x = sum(X); explode = zeros(size(x)); [c,offset] = max(x); explode(offset) = 1; h = pie(x,explode) %A = [ 1 3 6]; %pie3(A) Khi tổng các phần tử trong đối số thứ nhất bằng hay lớn hơn 1, pie và pie3 chuẩn hoá các giá trị. Như vậy cho vec tơ x, mỗi phần có diện tích xi / sum(xi ) với xi là một phần tử của x. Giá trị được chuẩn hoá mô tả phần nguyên của mỗi vùng. Khi tổng các phần tử trong đối số thứ nhất nhỏ hơn 1, pie và pie3 không chuẩn hoá các phần tử của vec tơ x. Chúng vẽ một phần pie. x = [.19 .22 .41]; pie(x) g. Làm hình chuyển động: Ta có thể tạo ra hình chuyển động bằng 2 cách • tạo và lưu nhiều hình khác nhau và lần lượt hiển thị chúng • vẽ và xoá liên tục một đối tượng trên màn hình,mỗi lần vẽ lại có sự thay đổi. Với cách thứ nhất ta thực hiện hình chuyển động qua 3 bước: • dùng hàm moviein để dành bộ nhớ cho một ma trận đủ lớn nhằm lưu các khung hình. • dùng hàm getframes để tạo các khung hình. • dùng hàm movie để hiển thị các khung hình. Sau đây là ví dụ sử dụng movie để quan sát hàm fft(eye(n)).Ta tạo chương trình ct1_22.m như sau : axis equal M = moviein(16, gcf); set(gca, ʹNextPlotʹ, ʹreplacechildrenʹ) h = uicontrol(ʹstyleʹ, ʹsliderʹ, ʹpositionʹ,[100 10 500 20], ʹMinʹ, 1, ʹMaxʹ, 16) for j = 1:16 plot(fft(eye(j + 16))) 21
set(h, ʹValueʹ, j) M(:, j) = getframe(gcf); end clf; axes(ʹPositionʹ, [0 0 1 1]); movie(M, 30) Bước đầu tiên để tạo hình ảnh chuyển động là khởi gán ma trận. Tuy nhiên trước khi gọi hàm moviein, ta cần tạo ra các trục toạ độ có cùng kích thước với kích thước mà ta muốn hiển thị hình. Do trong ví dụ này ta hiển thị các số liệu cách đều trên vòng tròn đơn vị nên ta dùng lệnh axis equal để xác định tỉ lệ các trục. Hàm moviein tạo ra ma trận đủ lớn để chứa 16 khung hình. Phát biểu: set(gca, ʹNextPlotʹ, ʹreplacechildrenʹ) ngăn hàm plot đưa tỉ lệ các trục về axis normal mỗi khi nó được gọi. Hàm getframe không đối số trả lại các điểm ảnh của trục hiện hành ở hình hiện có. Mỗi khung hình gồm các số liệu trong một vec tơ cột. Hàm getframe(gcf) chụp toàn bộ phần trong của một cửa sổ hiện hành. Sau khi tạo ra hình ảnh ta có thể chạy chúng một số lần nhất định ví dụ 30 lần nhờ hàm movie(M, 30) . Một phương pháp nữa để tạo hình chuyển động là vẽ và xoá, nghĩa là vẽ một đối tượng đồ hoạ rồi thay đổi vị trí của nó bằng cách thay đổi toạ độ x, y và z một lượng nhỏ nhờ một vòng lặp. Ta có thể tạo ra các hiệu ứng khác nhau nhờ các cách xoá hình khác nhau. Chúng gồm: • none MATLAB không xoá đối tượng khi nó di chuyển • background MATLAB xoá đối tượng bằng cách vẽ nó có màu nền • xor MATLAB chỉ xoá đối tượng Ta tạo ra M‐file có tên là ct1_23.m như sau: A = [ ‐8/3 0 0; 0 ‐10 10; 0 28 ‐1 ]; y = [35 ‐10 ‐7]ʹ; h = 0.01; p = plot3(y(1), y(2), y(3),ʹ.ʹ, ʹEraseModeʹ, ʹnoneʹ, ʹMarkerSizeʹ, 5); axis([0 50 ‐25 25 ‐25 25]) 22
hold on for i = 1:4000 A(1,3) = y(2); A(3,1) = ‐y(2); ydot = A*y; y = y + h*ydot; set(p, ʹXDataʹ, y(1), ʹYDataʹ, y(2), ʹZDataʹ, y(3)) % thay doi toa do drawnow i = i + 1; end 13. Đồ hoạ 3D: a.Các lệnh cơ bản: Lệnh mesh và surf tạo ra lưới và mặt 3D từ ma trận số liệu. Gọi ma trận số liệu là z mà mỗi phần tử của nó z(i, j) xác định tung độ của mặt thì mesh(z) tạo ra một lưới có màu thể hiện mặt z còn surf(z) tạo ra một mặt có màu z. b. Đồ thị các hàm hai biến: Bước thứ nhất để thể hiện hàm 2 biến z=f(x,y) là tạo ma trận x và y chứa các toạ độ trong miền xác định của hàm. Hàm meshgrid sẽ biến đổi vùng xác định bởi 2 vec tơ x và y thành ma trận x và y. Sau đó ta dùng ma trận này để đánh giá hàm. Ta khảo sát hàm sin(r)/r. Để tính hàm trong khoảng ‐8 và 8 theo x và y ta chỉ cần chuyển một vec tơ đối số cho meshgrid: [x,y] = meshgrid(‐8:.5:8); r = sqrt(x.^2 + y.^2) + 0.005; ma trận r chứa khoảng cách từ tâm của ma trận. Tiếp theo ta dùng hàm mesh để vẽ hàm. z = sin(r)./r; mesh(z) c. Đồ thị đường đẳng mức: Các hàm contour tạo, hiển thị và ghi chú các đường đẳng mức của một hay nhiều ma trận. Chúng gồm: clabel tạo các nhãn sử dụng ma trận contour và hiển thị nhãn contour hiển thị các đường đẳng mức tạo bởi một giá trị cho trước của ma trận Z. 23
contour3 hiển thị các mặt đẳng mức tạo bởi một giá trị cho trước của ma trận Z. contourf hiển thị đồ thị contour 2D và tô màu vùng giữa 2 các đường contourc hàm cấp thấp để tính ma trận contour Hàm meshc hiển thị contour và lưới và surfc hiển thị mặt contour. [X,Y,Z] = peaks; contour(X,Y,Z,20) Mỗi contour có một giá trị gắn với nó. Hàm clabel dùng giá trị này để hiển thị nhãn đường đồng mức 2D. Ma trận contour chứa giá trị clabel dùng cho các đường contour 2D. Ma trận này được xác định bởi contour, contour3 và contourf. Để hiển thị 10 đường đẳng mức của hàm peak ta viết: Z = peaks; [C,h] = contour(Z,10); clabel(C,h) title({ʹCac contour co nhanʹ,ʹclabel(C,h)ʹ}) Hàm contourf hiển thị đồ thị đường đẳng mức trên một mặt phẳng và tô màu vùng còn lại giữa các đường đẳng mức. Để kiểm soát màu tô ta dùng hàm caxis và colormap. Ta viết chương trình ct1_26.m: Z = peaks; [C, h] = contourf(Z, 10); caxis([‐20 20]) colormap autumn; title({ʹContour co to mauʹ, ʹcontourf(Z, 10)ʹ}) Các hàm contour(z, n) và contour(z, v) cho phép ta chỉ rõ số lượng mức contour hay một mức contour cần vẽ nào đó với z là ma trận số liệu, n là số đường contour và v là vec tơ các mức contour. MATLAB không phân biệt giữa vec tơ một phần tử hay đại lượng vô hướng. Như vậy nếu v là vec tơ một phần tử mô tả một contour đơn ở một mức hàm contour sẽ coi nó là số lượng đường contour chứ không phải là mức contour. Nghĩa là, contour(z, v) cũng như contour(z, n). Để hiển thị một đường đẳng mức ta cần cho v là một 24
vec tơ có 2 phần tử với cả hai phần tử bằng mức mong muốn. Ví dụ để tạo ra một đường đẳng mức 3D của hàm peaks ta viết chương trình ct1_27.m: xrange = ‐3: .125: 3; yrange = xrange; [X,Y] = meshgrid(xrange, yrange); Z = peaks(X, Y); contour3(X, Y, Z) Để hiển thị một mức ở Z = 1, ta cho v là [1 1] v = [1 1] contour3(X, Y, Z, v) Hàm ginput cho phép ta dùng chuột hay các phím mũi tên để chọn các điểm vẽ. Nó trả về toạ độ của vị trí con trỏ. Ví dụ sau sẽ minh hoạ các dùng hàm ginput và hàm spline để tạo ra đường cong nội suy hai biến. Ta tạo một M‐file có tên ct1_28.m như sau: disp(ʹChuot phai tro cac diem tren duong veʹ) disp(ʹChuot trai tro diem cuoi cua duong veʹ) axis([0 10 0 10]) hold on x = []; y = []; n = 0; but = 1; while but = =1 [xi,yi,but] = ginput(1); plot(xi, yi, ʹgoʹ) n = n + 1; x(n, 1) = xi; y(n,1) = yi; end t = 1:n; ts = 1: 0.1: n; xs = spline(t, x, ts); 25
ys = spline(t, y, ts); plot(xs, ys, ʹc‐ʹ); hold off 14. Vẽ các vectơ: Có nhiều hàm MATLAB dùng hiển thị các vec tơ có hướng và vec tơ vận tốc. Ta định nghĩa một vec tơ bàng cách dùng một hay 2 đối số. Các đối số mô tả thành phần x và thành phần y của vec tơ. Nếu ta dùng 2 đối số thì đối số thứ nhất sẽ mô tả thành phần x và đối số thứ ha mô tả thành phần y. Nếu ta chỉ dùng một đối số thì MATLAB xử lí nó như một số phức, phần thực là thành phần x và phần ảo là thành phần y. Các hàm vẽ vec tơ gồm: compass vẽ các véc tơ bắt đầu từ gốc toạ độ của hệ toạ độ cực feather vẽ các vec tơ bắt đầu từ một đường thẳng quiver vẽ các vec tơ 2D có các thành phần (u, v) quiver3 vẽ các vec tơ 3D có các thành phần (u, v, w) a. Hàm compass: Ta xét ví dụ vẽ hướng và tốc độ gió. Các vec tơ xác định hướng (góc tính bằng độ) và tốc độ gió (km/h) là: hg = [45 90 90 45 360 335 360 270 335 270 335 335]; td = [6 6 8 6 3 9 6 8 9 10 14 12]; Ta biến đổi hướng gió thành radian trước khi biến đổi nó thành toạ độ vuông góc. hg1 = hg * pi/180; [x, y] = pol2cart(hg1, td); compass(x, y) và tạo ra ghi chú trên đồ thị: gc = {ʹHuong gio và suc gio tai san bay Da Nangʹ) text(–28, 15, gc) b. Hàm feather: Hàm feather hiển thị các vec từ bắt đầu từ một đường thẳng song song với trục x. Ví dụ để tạo ra các vec tơ có góc từ 900 đến 00 và cùng độ dài ta viết chương trình ct1_30.m: theta = 90: –10: 0; 26
r = ones(size(theta)); trước khi vẽ, chuyển các số liệu sang toạ độ vuông góc và tăng độ lớn thành r để dễ nhìn: [u, v] = pol2cart(theta*pi/180, r*10); feather(u, v) axis equal Nếu đối số là số phức z thì feather coi phần thực là x và phần ảo là y. Ta xét chương trình ct1_31.m: t = 0: 0.3: 10; s = 0.05 + i; Z = exp(–s*t); feather(Z) c. Hàm quiver: Hàm quiver hiển thị các vec tơ ở các điểm đã cho trong mặt phẳng. Các vec tơ này được xác định bằng các thành phần x và y. Ví dụ để tạo ra 10 contour của hàm peaks ta dùng chương trình ct1_32.m: n = –2.0: .2: 2.0; [X,Y,Z] = peaks(n); contour(X, Y, Z, 10) Bây giờ dùng hàm gradient để tạo các thành phần của vec tơ dùng làm đối số cho quiver: [U, V] = gradient(Z, .2); Đặt hold on để thêm đường contour: hold on quiver(X,Y,U,V) hold off 27
d. Hàm quiver3: Hàm quiver3 hiển thị các vec tơ có các thành phần (u,v,w) tại điểm (x, y, z). Ví dụ ta biểu diễn quỹ đạo của một vật được ném đi theo t. Phương trình của chuyển động là: at2 z(t) = v t + 0 2 Ta viết chương trình ct1_33.m. Trước hết ta gán vận tốc ban đầu và gia tốc a: v0 = 10; % Van toc ban dau a = –32; % gia toc Tiếp theo tính z tại các thời điểm: t = 0:.1:1; z = vz*t + 1/2*a*t.^2; Tính vị trí theo hướng x và y: vx = 2; x = vx*t; vy = 3; y = vy*t; Tính các thành phần của vec tơ vận tốc và hiển thị bằng các dùng quiver3: u = gradient(x); v = gradient(y); w = gradient(z); scale = 0; quiver3(x, y, z, u, v, w, scale) axis square §3. GIAO DIỆN ĐỒ HOẠ 1. Khái niệm chung: Để tiện dụng ta có thể tạo nên giao diện đồ hoạ(GUI ‐ Graphic User Interface) giữa người dùng và MATLAB. Trong giao diện này ta có thể xuất dữ liệu dưới 2 dạng: văn bản và đồ hoạ. Mỗi một GUI có một hay nhiều layout(diện mạo). Việc tạo GUI tạo nên một công cụ đồ hoạ phục vụ 28
nhập xuất dữ liệu một cách trực giác, rất thuận tiện. Ngoài ra có thể dùng GUI để giám sát các quá trình, hiển thị các đối tượng. 2. Nhập xuất kí tự, số liệu ra GUI: a. Tạo khung hình: Ta xét các lệnh sau(ct1_35.m): f = input(ʹNhap nhiet do(do K): ʹ); c = (f ‐ 32)*5/9; fprintf(1,ʹnhiet do(do C) la: %g\nʹ, c) Ba dòng lệnh trên thực hiện các công việc sau: ‐ nhập giá trị đầu vào ‐ thực hiện phép tính quy đổi nhiệt độ ‐ xuất kết quả ra màn hình Bây giờ ta tìm cách cài các dòng lệnh trên sao cho chúng thực hiện trên khuôn khổ một khung đồ hoạ có dạng như trên Các lệnh sau(ct1_36.m) thực hiện công việc trên: set(gcf,ʹDefaultUicontrolUnitʹ, ʹNormalizedʹ) frame_1 = uicontrol(gcf, ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.1 0.8 0.3]); frame_2 = uicontrol(gcf, ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.6 0.8 0.3]); set(frame_1, ʹBackgroundColorʹ, [0.5 0.5 0.5]); set(frame_2, ʹBackgroundColorʹ, [0.5 0.5 0.5]); 29
text_f = uicontrol(gcf, ʹStyleʹ, ʹTextʹ, ʹStringʹ, ʹFahrenheit: ʹ, ʹPositionʹ, [0.3 0.7 0.2 0.05],ʹHorizontalAlignmentʹ,ʹLeftʹ); edit_f = uicontrol(gcf, ʹStyleʹ, ʹEditʹ, ʹStringʹ, ʹ168.0ʹ, ʹPositionʹ, [0.6 0.7 0.1 0.05 ], ʹHorizontalAlignmentʹ, ʹRightʹ, ʹCallbackʹ, ʹct1_38ʹ); text_c1 = uicontrol(gcf,ʹStyleʹ, ʹTextʹ, ʹStringʹ, ʹCelcius: ʹ, ʹPositionʹ, [0.3 0.3 0.2 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); text_c2 = uicontrol(gcf,ʹStyleʹ, ʹTextʹ, ʹStringʹ, ʹ100.0ʹ, ʹPositionʹ, [0.6 0.3 0.1 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); Bây giờ ta sẽ xem các lệnh trên hoạt động như thế nào. Các lệnh sau: set(gcf,ʹDefaultUicontrolUnitʹ, ʹNormalizedʹ) frame1 = uicontrol(gcf,ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.1 0.8 0.3]); frame2 = uicontrol(gcf,ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.6 0.8 0.3]); set(frame1,ʹBackgroundColorʹ, [0.5 0.5 0.5]); set(frame2,ʹBackgroundColorʹ, [0.5 0.5 0.5]); tạo hai khung hình chữ nhật trong cửa sổ Figure hiện hành với nền màu xám. Hai khung (Frames) có toạ độ các góc dưới trái là (0.1, 0.1) và (0.1, 0.6), cùng chiều cao 0.3 đơn vị và bề rộng 0.8 đơn vị. Đơn vị được tính bằng % của kích cỡ ngoài của Figure. Vậy ta có thể diễn giải như sau: ‐ Khung thứ nhất có góc trái dưới tại điểm có toạ độ 10% chiều ngang và 10% chiều cao của khung ngoài Figure. ‐ Khung thứ 2 có góc trái phía dưới tại điểm có toạ độ ứng với 10% chiều ngang và 60% chiều cao của khung ngoài Figure. ‐ Cả hai khung có chiều cao bằng 30% chiều cao và bề ngang bằng 80% bề ngang của khung ngoài Figure. 30
b. Dùng lệnh edit và text để nhập xuất kí tự và số liệu: Trên đây ta đã dùng lệnh uicontrol để tạo và xác định vị trí hai khung hình. Đoạn lệnh sau sử dụng uicontrol để viết chuỗi kí tự “Fahrenheit” lên khung bên trên: text_ f = uicontrol(gcf,ʹStyleʹ,ʹTextʹ,ʹStringʹ,ʹFahrenheit: ʹ, ʹPositionʹ,[0.3 0.7 0.2 0.05],ʹHorizontalAlignmentʹ,ʹLeftʹ); Chuỗi kí tự “Fahrenhaeit” được đặt vào đúng vị trí dồn trái của ô có Position ghi trong đoạn chương trình trên. Đoạn lệnh sau dùng Edit để viết chuỗi kí tự “68.0” vào vị trí bên cạnh của “Fahrenheit”. Chuỗi kí tự có vị trí dồn phải trong ô (Position Box). edit_f = uicontrol(gcf,ʹStyleʹ, ʹEditʹ, ʹStringʹ, ʹ168.0ʹ, ʹPositionʹ, [0.6 0.7 0.1 0.05 ], ʹHorizontalAlignmentʹ, ʹRightʹ, ʹCallbackʹ, ʹct1_38ʹ); Do sử dụng edit, chuỗi kí tự “68.0” là chuỗi có thể viết lại được trực tiếp trên GUI. Sau khi nhấn nút trên, giá trị mới viết lại được tiếp nhận và MATLAB sẽ gọi lệnh viết trong phần callback ct1_38.m. Cuối cùng ta còn phải dùng uicontrol để tạo ta chuỗi text, hiển thị chuỗi “Celcius” và “20.0” trong khung bên dưới. text_c1 = uicontrol(gcf,ʹStyleʹ,ʹTextʹ,ʹStringʹ,ʹCelcius: ʹ, ʹPositionʹ,[0.3 0.3 0.2 0.05],ʹHorizontalAlignmentʹ,ʹLeftʹ); text_c2 = uicontrol(gcf,ʹStyleʹ,ʹTextʹ,ʹStringʹ,ʹ20.0ʹ,ʹPositionʹ, [0.6 0.3 0.1 0.05],ʹHorizontalAlignmentʹ,ʹLeftʹ); c. Tự động cập nhật giá trị lên GUI: Để hoàn thiện ví dụ GUI ta thực hiện chương trình với nhiệm vụ tính quy đổi từ độ K sang độ C và tự động điền kết quả vào các ô bên cạnh chuỗi Celcius. Đoạn mã sau phục vụ mục đích callback (hoàn trả giá trị) được lưu vào file ct1_38.m và có nội dung như sau: f = get(edit_f, ʹStringʹ); f = str2num(f); c = ( f ‐ 32)*5/9; 31
c = num2str(c); set(text_c2, ʹStringʹ,c); Đoạn mã trên nhận giá trị do lệnh uicontrol “edit” đọc vào dưới dạng chuỗi (string) và sau đó: ‐ biến đổi từ dạng string sang dạng số ‐ tính quy đổi từ nhiệt độ fahrenheit sang nhiệt độ celcius ‐ biến đổi từ số sang string ‐ xuất kết quả dưới dạng string ra GUI nhờ text_c2 3. Nhập số liệu từ thanh trượt: Ngoài cách nhập số liệu từ bàn phím, ta có thể nhập số liệu từ thanh trượt. Ta muốn tạo ra một giao diện như sau: Trong giao diện này, con trượt sẽ làm thay đổi giá trị nhiệt độ đua vào và nhiệt độ quy đổi tính theo độ C cũng sẽ thay đổi tương ứng. Các lệnh tạo ra giao diện này (ct1_37.m) là: set(gcf, ʹDefaultUicontrolUnitʹ, ʹNormalizedʹ) frame_1 = uicontrol(gcf, ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.1 0.8 0.3]); frame_2 = uicontrol(gcf, ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.6 0.8 0.3]); set(frame_1, ʹBackgroundColorʹ ,[0.5 0.5 0.5]); set(frame_2, ʹBackgroundColorʹ, [0.5 0.5 0.5]); text_ f = uicontrol(gcf, ʹStyleʹ, ʹTextʹ, ʹStringʹ, ʹFahrenheit: ʹ,ʹPositionʹ, [0.3 0.7 0.2 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); edit_f = uicontrol(gcf, ʹStyleʹ, ʹEditʹ, ʹStringʹ, ʹ168.0ʹ.,,, 32
ʹPositionʹ, [0.6 0.7 0.1 0.05 ], ʹHorizontalAlignmentʹ, ʹRightʹ, ʹCallbackʹ, ʹct1_38ʹ); text_c1 = uicontrol(gcf,ʹStyleʹ, ʹTextʹ, ʹStringʹ, ʹCelcius: ʹ, ʹPositionʹ, [0.3 0.3 0.2 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); text_c2 = uicontrol(gcf,ʹStyleʹ, ʹTextʹ, ʹStringʹ, ʹ100.0ʹ, ʹPositionʹ, [0.6 0.3 0.1 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); slider_f = uicontrol(gcf,ʹStyleʹ, ʹSliderʹ, ʹMinʹ, 32.0, ʹMaxʹ, 212.0, ʹValueʹ, 68.0, ʹPositionʹ, [0.6 0.8 0.2 0.05], ʹCallbackʹ, ʹct1_39; ct1_38ʹ); Để tạo thanh trượt ta dùng lệnh: slider_f = uicontrol(gcf,ʹStyleʹ,ʹSliderʹ,ʹMinʹ,32.0,ʹMaxʹ,212.0, ʹValueʹ,68.0,ʹPositionʹ,[0.6 0.8 0.2 0.05], ʹCallbackʹ,ʹct1_39; ct1_38ʹ); Như vậy Callback có thể gọi một chuỗi các lệnh MATLAB, phân cách nhau bằng dấu chấm than hay dấu phẩy. Chuỗi callback gọi ct1_39.m: f = get(slider_f,ʹValueʹ); f = num2str(f); set(edit_f,ʹStringʹ,f,ʹCallBackʹ,ʹct1_40; ct1_38ʹ); với tác dụng nhập nhiệt độ giữ tại ‘Value’ của slider_f vào vị trí bên cạnh ô chứa chuỗi “Fahrenheit”. Sau đó Callback gọi tiếp ct1_38.m để tính quy đổi giá trị nhiệt độ và gán vào ô cạnh chuỗi “Celcius”. File ct1_40.m như sau: f = get(edit_f,ʹStringʹ); f = str2num(f); set(slider_f,ʹValueʹ,f); 33
có nhiệm vụ cập nhật giá trị giữ tại ‘Value’ của slider_f để rồi sau đó ct1_38.m làm nốt phần việc còn lại: tính đổi nhiệt độ và gán vào vị trí cạnh ô chứa chuỗi “Celcius”. 4. Chọn lựa khi xuất số liệu: a. Khái niệm chung: Ngoài khả năng xuất dữ liệu cố định theo kiểu string hay kiểu số, ta có thể xuất dữ liệu theo một danh mục nào đó. Để minh hoạ, ta tạo file ct1_41.m như sau: f = input(ʹNhap nhiet do: ʹ); r = f + 459.7; c = (f ‐ 32)*5/9; k = c + 273.15; choice = input([ʹNhap 1 cho Rankieʹ, ʹ2 cho Celciusʹ, ʹ3 cho Kelvin: ʹ]); if choice = = 1 fprintf(1, ʹNhiet do (do R) la: %g\nʹ, r); elseif choice = = 2 fprintf(2, ʹNhiet do (do C) la: %g\nʹ, c); elseif choice = = 3 fprintf(2, ʹNhiet do (do C) la: %g\nʹ, c); end Từ cửa sổ lệnh, nhập lệnh ct1_41 thì MATLAB sẽ hỏi nhiệt độ và đích quy đổi rồi hiển thị kết quả. Tuy nhiên công cụ GUI của MATLAB cho phép ta thực hiện việc lựa chọn thuận lợi hơn. Ta có thể chọn một trong 4 phương xuất dữ liệu sau đây: ‐ dùng popupmenu ‐ dùng list box ‐ dùng radio button ‐ dùng check box b. Dùng popupmenu: Ta tạo ra giao diện như sau: 34
Các lệnh thực hiện công việc trên (ct1_42.m) là: set(gcf, ʹDefaultUicontrolUnitʹ, ʹNormalizedʹ) frame_1 = uicontrol(gcf, ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.1 0.8 0.3]); frame_2 = uicontrol(gcf, ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.6 0.8 0.3]); set(frame_1, ʹBackgroundColorʹ, [0.5 0.5 0.5]); set(frame_2, ʹBackgroundColorʹ ,[0.5 0.5 0.5]); text_f = uicontrol(gcf,ʹStyleʹ, ʹTextʹ, ʹStringʹ, ʹFahrenheit: ʹ, ʹPositionʹ, [0.3 0.7 0.2 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); edit_f = uicontrol(gcf,ʹStyleʹ, ʹEditʹ, ʹStringʹ, ʹ168.0ʹ, ʹPositionʹ, [0.6 0.7 0.1 0.05 ], ʹHorizontalAlignmentʹ, ʹRightʹ, ʹCallbackʹ, ʹct1_38ʹ); popup_c = uicontrol(gcf, ʹStyleʹ,ʹPopupmenuʹ, ʹStringʹ,ʹRankine|Celcius|Kelvinʹ, ʹValueʹ,2, ʹPositionʹ,[0.3 0.3 0.2 0.05], ʹCallbackʹ,ʹct1_43; ct1_45ʹ); text_c2 = uicontrol(gcf, ʹStyleʹ, ʹTextʹ, 35
ʹStringʹ, ʹ100.0ʹ, ʹPositionʹ, [0.6 0.3 0.1 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); slider_f = uicontrol(gcf, ʹStyleʹ, ʹSliderʹ, ʹMinʹ, 32.0, ʹMaxʹ, 212.0, ʹValueʹ, 68.0, ʹPositionʹ, [0.6 0.8 0.2 0.05], ʹCallbackʹ, ʹct1_39; ct1_45ʹ); Khi kích chuột vào Popupmenu , có ba khả năng chọn lựa sẽ xuất hiện. Tiếp tục nháy chuột vào một trong 3 khả năng đó , Popupmenu biến mất chỉ còn lại đơn vị được chọn. Khi dùng chuột kéo thanh trượt ở frame phía trên, ta có được giá trị quy đổi sang đơn vị được chọn hiển thị ở phía dưới. Trong đoạn mã trên, giá trị ‘Value’ đặt sẵn là 2. Khi Callback gọi ct1_43.m: choice = get(popup_c,’Value’); thì giá trị của biến choice được đưa tới ‘Value’. Sau đó Callback gọi tiếp ct1_45.m để xem kết quả giữ trong choice. File ct1_45.m như sau: f = get(edit_f, ʹStringʹ); f = str2num(f); r = f + 459.7; c = (f ‐ 32)*5/9; k = c + 273.15; choice = input([ʹNhap 1 cho Rankieʹ, ʹ2 cho Celciusʹ, ʹ3 cho Kelvin: ʹ]); if choice = = 1 t = r; elseif choice = = 2 t = c; elseif choice = = 3 t = k end t = num2str(t); set(text_c2, ʹStringʹ,t); 36
Bằng cách thay ‘Popupmenu’ bằng ‘Radiobutton’ uicontrol ta có phương án Radiobutton. Giao diện sẽ có dạng: Các lệnh thực hiện công việc này (ct1_46.m) là: set(gcf, ʹDefaultUicontrolUnitʹ, ʹNormalizedʹ) frame_1 = uicontrol(gcf, ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.1 0.8 0.3]); frame_2 = uicontrol(gcf, ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.6 0.8 0.3]); set(frame_1,ʹBackgroundColorʹ, [0.5 0.5 0.5]); set(frame_2,ʹBackgroundColorʹ, [0.5 0.5 0.5]); text_f = uicontrol(gcf, ʹStyleʹ, ʹTextʹ, ʹStringʹ, ʹFahrenheit: ʹ,ʹPositionʹ, [0.3 0.7 0.2 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); edit_f = uicontrol(gcf, ʹStyleʹ, ʹEditʹ, ʹStringʹ,ʹ168.0ʹ, ʹPositionʹ, [0.6 0.7 0.1 0.05 ], ʹHorizontalAlignmentʹ, ʹRightʹ, ʹCallbackʹ,ʹct1_41ʹ); strings = [ʹRankineʹ; ʹCelciusʹ; ʹKelvineʹ]; show = [ 0; 1; 0]; ys = [ 3; 2; 1]*0.075 + 0.075; for i = 1:3 radio_c(i) = uicontrol(gcf, ʹStyleʹ, ʹRadiobuttonʹ, 37
ʹStringʹ, strings(i), ʹValueʹ, show(i), ʹPositionʹ, [0.3 ys(i) 0.2 0.05], ʹCallbackʹ, ʹ ct1_47; ct1_45ʹ); end text_c2= uicontrol(gcf, ʹStyleʹ, ʹTextʹ, ʹStringʹ,ʹ100.0ʹ, ʹPositionʹ, [0.6 0.3 0.1 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); slider_f = uicontrol(gcf, ʹStyleʹ, ʹSliderʹ, ʹMinʹ,32.0, ʹMaxʹ, 212.0, ʹValueʹ, 68.0, ʹPositionʹ, [0.6 0.8 0.2 0.05], ʹCallbackʹ, ʹct1_39; ct1_45ʹ); File ct1_47.m: for i = 1:3 if gcbo = = radio_c(i) choice = i; set(radio_c(i), ʹValueʹ, 1); elseif set(radio_c(i), ʹValueʹ, 0); end; end; Đoạn lệnh trên là một vòng lặp, so sánh số (handle) Callback thu được (giá trị do hàm gcbo trả về) với handle của mỗi nút. Nút nào có số trùng sẽ được đóng (turn on, ‘Value’ = 1) và nút nào khác số sẽ bị ngắt (turn off,’Value’ = 0). Cuối cùng Callback gọi ct1_45.m để thực hiện việc tính quy đổi được chọn và hiển thị kết quả. Điểm khác duy nhất là khi chọn, Popupmenu chỉ chứa một phần tử thì radiobutton có thể đồng thời chứa nhiều phần tử. Cuối cùng ta xét phương án dùng listbox. Giao diện cần tạo như sau: 38
Các mã tạo ra giao diện trên (ct1_48.m) là: set(gcf, ʹDefaultUicontrolUnitʹ, ʹNormalizedʹ) frame_1 = uicontrol(gcf, ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.1 0.8 0.3]); frame_2 = uicontrol(gcf, ʹStyleʹ, ʹFrameʹ, ʹPositionʹ, [0.1 0.6 0.8 0.3]); set(frame_1, ʹBackgroundColorʹ, [0.5 0.5 0.5]); set(frame_2, ʹBackgroundColorʹ, [0.5 0.5 0.5]); text_f = uicontrol(gcf, ʹStyleʹ, ʹTextʹ, ʹStringʹ, ʹFahrenheit: ʹ, ʹPositionʹ, [0.3 0.7 0.2 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); edit_f = uicontrol(gcf, ʹStyleʹ, ʹEditʹ, ʹStringʹ, ʹ168.0ʹ, ʹPositionʹ, [0.6 0.7 0.1 0.05 ], ʹHorizontalAlignmentʹ, ʹRightʹ, ʹCallbackʹ, ʹct1_38ʹ); listbox_c = uicontrol(gcf, ʹStyleʹ, ʹListboxʹ, ʹStringʹ, ʹRankine|Celcius|Kelvinʹ, ʹValueʹ, 2, ʹPositionʹ, [0.3 0.3 0.2 0.05], ʹCallbackʹ, ʹct1_49;ct1_45ʹ); text_c2 = uicontrol(gcf, ʹStyleʹ, ʹTextʹ, ʹStringʹ, ʹ100.0ʹ, ʹPositionʹ, [0.6 0.3 0.1 0.05], ʹHorizontalAlignmentʹ, ʹLeftʹ); slider_f = uicontrol(gcf, ʹStyleʹ, ʹSliderʹ, ʹMinʹ,32.0, ʹMaxʹ, 212.0, ʹValueʹ, 68.0, ʹPositionʹ, [0.6 0.8 0.2 0.05], ʹCallbackʹ, ʹct1_39; ct1_45ʹ); 5. Công cụ đồ hoạ tạo GUI 39
a. Tạo GUI bằng công cụ đồ hoạ: Trên đây ta đã xem xét cách tạo GUI bằng phương pháp thủ công. Ta có thể tạo GUI bằng công cụ đồ hoạ. Khi nhập lệnh guide ta gọi trình đồ hoạ (Graphics User Interface Development Environment) để soạn thảo layout. Kết quả đầu tiên là ta có một layout rỗng như sau: Soạn thảo Alignment thuộc tính Chạy thử Soạn menu Vùng thiết kế Các phần tử Việc đầu tiên là ta thiết kế giao diện mong muốn. Ta sẽ dùng chuột kéo các phần tử cần dùng từ bên trái và thả vào layout rỗng bên phải. Ta có thể dịch chuyển các phần tử này đế các vị trí mong muốn và cân chỉnh bằng công cụ Alignment. Với mỗi phần tử ta cấn xác định thuộc tính cho nó bằng cách bấm đúp vào phần tử hay bấm vào công cụ soạn thảo thộc tính Sau khi thiết kế xong ta lưu nó lại. Lúc này MATLAB tự động tạo ra file *.fig dùng lưu giao diện vừa tạo và file *.m chưa các mã lệnh cần thực hiện. Việc cuối cùng là viết các mã lệnh vào file *.m. Trong quá trình thiết kế ta có thể chạy thử xem sau mỗi bước thiết kế đã đạt yêu cầu chưa bằng cách bấm vào ô chạy thử b. Một số ví dụ tạo GUI: ) Đếm số lần bấm chuột: Ta thiết kế một giao diện như sau: 40
Ta muốn là khi bấm chuột, số lần bấm sẽ được đếm và ghi lại. Trước hết ta gọi guide và có được một layout rỗng. Vào Property Inspector (ô soạn thảo thuộc tính) và ghi vào Name chuỗi ʺct1_52ʺ và chấp nhận thuộc tích Tag mặc định của nó là figure1; dùng Font chữ mặc định, cỡ chữ 12, bold. Ta dùng ô Edit Text để ghi lại số lần bấm. Ta vào Property Inspector rồi chọn String. Ta nhập vào ô này chuỗi ʺSo lan bam chuot: 0ʺ. Ta ghi vào ô Tag chuỗi ʺeditmotʺ và cũng dùng Font chữ mắc định, cỡ chữ 12 và bold. Tiếp theo kéo Pushbutton vào layout và soạn thảo thuộc tính cho nó với Font chữ mặc định, cỡ chứ 12, bold. Trong thuôc tính String ghi chuỗi ʺ Bam chuotʺ; ghi và Tag chuỗi ʺpushbuttonmotʺ. Như vậy là ta đã thiết kế xong. Bây giờ ta lưu lại với tên là ct1_52.fig và ct1_52.m. Nhiệm vụ tiếp theo là ghi các lệnh cần thiết vào file ct1_52.m. File này đã được MATLAB tự động tạo ra. Ta phải thêm vào đó các mã lệnh để khi bấm chuột thì số lần bấm được thể hiện trên ô Edit Text. Ta sẽ ghi các mã lệnh này vào phần: function varargout = pushbuttonmot_Callback(h, eventdata, handles, varargin) do lệnh cần được thực hiện khi gọi pushbutton. Nội dung của ct1_52.m là: function varargout = Ct1_52(varargin) if nargin = = 0 fig = openfig(mfilename,ʹreuseʹ); set(fig, ʹColorʹ, get(0, ʹdefaultUicontrolBackgroundColorʹ)); 41
handles = guihandles(fig); guidata(fig, handles); if nargout > 0 varargout{1} = fig; end elseif ischar(varargin{1}) try [varargout{1:nargout}] = feval(varargin{:}); catch disp(lasterr); end end function varargout = pushbuttonmot_Callback(h, eventdata, handles, varargin) persistent dem;%bien dem la persistent de no ton tai giua lan goi ham if isempty(dem) dem = 0; end dem = dem + 1; str = sprintf(ʹSo lan bam chuot: %dʹ,dem); set(handles.editmot,ʹStringʹ,str); ) Chuyển đổi từ độ Fahrenheit sang độ Celcius: Ta thiết kế một GUI để chuyển đổi nhiệt độ. Giao diện có dạng như sau: Thuộc tính của Layout được ghi Name: ct1_53 còn các thuộc tính khác là mặc định. 42
Ta dùng hai Frame với các Tag là frmmot và frame2. Các thuộc tính khác chấp nhận giá trị mặc định. Edit Text thứ nhất có các thuộc tính FontName: Arial, FontSize: demi, FơntWeight: demi, String: Fahrenheit, Tag: editmot còn các thuộc tính khác là mặc định. Edit Text thứ hai có các thuộc tính FontName: Arial, FontSize: demi, FơntWeight: demi, String: để trống, Tag: edithai còn các thuộc tính khác là mặc định. Edit Text thứ ba có các thuộc tính FontName: Arial, FontSize: demi, FơntWeight: demi, String: Celcius, Tag: editba còn các thuộc tính khác là mặc định. Edit Text thứ tư có các thuộc tính FontName: Arial, FontSize: demi, FơntWeight: demi, String: để trống, Tag: editbon còn các thuộc tính khác là mặc định. Sau khi thiết kế xong, lưu nó với tên ct3_18.fig. MATLAB tạo thêm ct1_53.m. Bây giờ ta cần viết mã cho nó. Nhiệm vụ của đoạn mã là khi ta nhập nhiệt độ Fahrenheit vào ô Edit text thứ hai thì trong ô Edit Text thứ 4 phải xuất hiện giá trị nhiệt độ Celcius tương ứng. Do vậy nội dung của ct1_53.m là: function varargout = Ct1_53(varargin) if nargin == 0 % LAUNCH GUI fig = openfig(mfilename,ʹreuseʹ); set(fig,ʹColorʹ,get(0,ʹdefaultUicontrolBackgroundColorʹ)); handles = guihandles(fig); guidata(fig, handles); if nargout > 0 varargout{1} = fig; end elseif ischar(varargin{1}) try [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard catch disp(lasterr); end end function varargout = edithai_Callback(h, eventdata, handles, varargin) f = get(handles.edithai,ʹStringʹ); 43
f = str2num(f); c = (f ‐ 32)*5/9; c = num2str(c); set(handles.editbon,ʹStringʹ,c); Trong đó đoạn mã cần viết nằm trong đoạn: function varargout = edithai_Callback(h, evendata, handles, varargin) Các lệnh khác là do MATLAB tự động tạo ra. ) Dùng slider để nhập số liệu: Ta dùng ví dụ chuyển đổi nhiệt độ trên nhưng bây giờ sẽ thêm slider để thay đổi nhiệt độ đầu vào. Giao diện sẽ có dạng: Như vậy ta cần 5 phần tử, trong đó có một phần tử là slider và 4 phần tử Edit Text. Layout có thuộc tính Name: ct1_54, còn các thuộc tính khác ta chấp nhận giá trị mặc định. Slider có thuộc tính Max: 1.0 và Min: 0.0. Edit Text thứ nhất có thuộc tính FontSize: 12, FơntWeight: bold, String: Fahrenheit còn các thuộc tính khác chấp nhận giá trị mặc định. Edit Text thứ 2 có thuộc tính FontSize: 12, FơntWeight: bold, String: để trống. Edit Text thứ 3 có thuộc tính FontSize: 12, FơntWeight: bold, String: Celcius. 44
Edit Text thứ 4 có thuộc tính FontSize: 12, FơntWeight: bold, String: để trống. (Các thuộc tính mà ta không nhắc đến có nghĩa là chấp nhận giá trị mặc định). Layout được lưu với tên ct1_54.fig. Bây giờ ta viết mã cho phần ct1_54.m mà MATLAB đã tự động tạo ra. Nhiệm vụ của nó là nhận giá trị thay đổi từ con trượt, cập nhật cho Edit Text 2 và Edit Text 4. Ta có nội dung của ct1_54.m: function varargout = ct1_54(varargin) if nargin = = 0 fig = openfig(mfilename,ʹreuseʹ); handles = guihandles(fig); guidata(fig, handles); if nargout > 0 varargout{1} = fig; end elseif ischar(varargin{1}) try [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard catch disp(lasterr); end end function varargout = slider1_Callback(h, eventdata, handles, varargin) f = get(handles.slider1,ʹValueʹ);%nhan gia tri tu con truot f = f*180 + 32;%tinh ra do Fahrenheit a = num2str(f);%bien lai thanh chuoi set(handles.edit2,ʹStringʹ,a);%ghi vao o do Fahrenheit b = (f‐32)*5/9;%doi thanh do Celcius b = num2str(b);%doi lai thanh chuoi set(handles.edit4,ʹStringʹ,b);%ghi vao o do Celcius ) Xuất số liệu có lựa chọn: Ta vẫn dùng ví dụ trên nhưng bây giờ nhiệt độ quy đổi có thể được tính theo thang nhiệt độ Kenvine, Celcius hay 45
Rankine. Để có thể chọn lựa ta dùng một trong các phương án: Popupmenu, Rdiobutton, Listbox hay Checkbox. Giao diện khi dùng Popupmenu như sau: Như vậy là ta cần một Slider, ba Edit Text và một Popupmenu. Layout có thuộc tính Name: ct13_55. Slider có thuộc tính Max: 1 và Min: 0 Edit Text thứ nhất có thuộc tính FontSize: 12, FơntWeight: bold và String: Fahrenheit. Edit Text thứ hai có thuộc tính FontSize: 12, FơntWeight: bold và String để trống. Edit Text thứ 3 có thuộc tính FontSize: 12, FơntWeight: bold và String để trống. Popupmenu có thuộc tính FontSize: 12, FontWeight: bold. Để ghi vào thuộc tính String ta bấm đúp chuột vào icon của nó và viết 3 dòng: Kelvine, Celcius và Rankine. File được lưu với tên ct1_55.fig. Vấn đề còn lại là viết mã trong file ct1_55.m. Mã cần thực hiện nhận giá trị từ Slider, xem Popupmenu nào được chọn để hiển thị nhiệt độ tương ứng. File ct1_55.m như sau: function varargout = ct1_55(varargin) if nargin == 0 % LAUNCH GUI fig = openfig(mfilename,ʹreuseʹ); set(fig,ʹColorʹ,get(0,ʹdefaultUicontrolBackgroundColorʹ)); handles = guihandles(fig); guidata(fig, handles); 46
if nargout > 0 varargout{1} = fig; end elseif ischar(varargin{1}) try [varargout{1:nargout}] = feval(varargin{:}); catch disp(lasterr); end end function varargout = slider1_Callback(h, eventdata, handles, varargin) f = get(handles.slider1,ʹValueʹ); f = f*180 + 32; a = num2str(f); set(handles.edit2,ʹStringʹ,a); r = f + 495.7; c = (f ‐ 32)*5/9; k = c + 273.15; chon = get(handles.popupmenu1,ʹValueʹ); if chon = = 1 t = k; elseif chon = = 2 t = c; elseif chon = = 3 t = r; end t = num2str(t); set(handles.edit3,ʹStringʹ,t); Tiếp theo ta xét trường hợp dùng listbox. Thay vì dùng Popupmenu ta dùng Listbox. Các phần tử khác và thuộc tính của nó không thay đổi. Thuộc tính Name của Layout là ct1_56. Ta vào ô String của Listbox và ghi vào đó 3 dòng Kelvine, Celcius và Rankine. Giao diện như sau: 47
File được lưu với tên ct1_56.fig. Tiếp theo viết lệnh cho ct1_56.m. Ta có file này như sau: function varargout = ct1_56(varargin) if nargin = = fig = openfig(mfilename,ʹreuseʹ); set(fig,ʹColorʹ,get(0,ʹdefaultUicontrolBackgroundColorʹ)); handles = guihandles(fig); guidata(fig, handles); if nargout > 0 varargout{1} = fig; end elseif ischar(varargin{1}) try [varargout{1:nargout}] = feval(varargin{:}); catch disp(lasterr); end end function varargout = slider1_Callback(h, eventdata, handles, varargin) f = get(handles.slider1,ʹValueʹ); f = f*180 + 32; a = num2str(f); set(handles.edit2,ʹStringʹ,a); r = f + 495.7; 48
c = (f ‐ 32)*5/9; k = c + 273.15; chon = get(handles.listbox1,ʹValueʹ); if chon = = 1 t = k; elseif chon = = 2 t = c; elseif chon = = 3 t = r; end t = num2str(t); set(handles.edit3,ʹStringʹ,t); Ta tiếp tục xét phương án dùng Radiobutton. Giao diện có dạng: Ta dùng ba Radiobutton thay cho Listbox. Radiobutton thứ nhất có thuộc tính FontSize: 12, FơntWeight: bold và String: Rankine. Radiobutton thứ 2 có thuộc tính FontSize: 12, FơntWeight: bold và String: Celcius. Radibutton thứ 3 có thuộc tính FontSize: 12, FơntWeight: bold và String: Kelvine. Các phần tử khác và thuộc tính của chúng vẫn như cũ. Layout có thuộc tính Name: ct1_57. Lưu GUI với tên ct1_57.fig. Tiếp theo ta viết các mã lệnh trong ct1_57.m: function varargout = ct1_57(varargin) if nargin = = 0 fig = openfig(mfilename,ʹreuseʹ); set(fig,ʹColorʹ,get(0,ʹdefaultUicontrolBackgroundColorʹ)); handles = guihandles(fig); 49
guidata(fig, handles); if nargout > 0 varargout{1} = fig; end elseif ischar(varargin{1}) try [varargout{1:nargout}] = feval(varargin{:}); catch disp(lasterr); end end function mutual_exclude(off) set(off,ʹValueʹ,0); function varargout = slider1_Callback(h, eventdata, handles, varargin) global chon f = get(handles.slider1,ʹValueʹ); f = f*180 + 32; a = num2str(f); set(handles.edit2,ʹStringʹ,a); r = f + 495.7; c = (f ‐ 32)*5/9; k = c + 273.15; if chon = = 1 t = r; elseif chon = = 2 t = c; elseif chon == 3 t = k; end t = num2str(t); set(handles.edit3,ʹStringʹ,t); function varargout = radiobutton1_Callback(h, eventdata, handles, varargin) global chon; off = [handles.radiobutton2, handles.radiobutton3]; mutual_exclude(off); chon = 1; function varargout = radiobutton2_Callback(h, eventdata, handles, varargin) global chon; off = [handles.radiobutton1, handles.radiobutton3]; 50
mutual_exclude(off); chon = 2; function varargout = radiobutton3_Callback(h, eventdata, handles, varargin) global chon; off = [handles.radiobutton1, handles.radiobutton2]; mutual_exclude(off); chon = 3; on lnh: function mutual_exclude(off) set(off,'Value',0); lm cho 3 nút lnh tr thnh mt nhóm. Các câu lnh: off = [handles.radiobutton1, handles.radiobutton2]; mutual_exclude(off); lm cho khi chn mt nút Radiobutton ny thì không chn c nút khác na. Cui cùng ta xét phng án dùng Checkbox. Giao din nh sau: Cng nh phng án dùng Radiobutton, ta dùng ba Checkbox. Checkbox th nht có thuc tính FontSize: 12, FntWeight: bold v String: Rankine. Checkbox th hai có thuc tính FontSize: 12, FntWeight: bold v String: Celcius. Checkbox th ba có thuc tính FontSize: 12, FntWeight: bold v String: Kelvine. Các phn t khác không co gì thay i. Ta lu file vi tên ct1_58.fig. Tip theo ta vit mã lnh cho ct1_58.m: function varargout = ct1_58(varargin) if nargin = = 0 fig = openfig(mfilename,'reuse'); set(fig,'Color',get(0,'defaultUicontrolBackgroundColor')); handles = guihandles(fig); guidata(fig, handles); 51
if nargout > 0 varargout{1} = fig; end elseif ischar(varargin{1}) try [varargout{1:nargout}] = feval(varargin{:}); catch disp(lasterr); end end function mutual_exclude(off) set(off,'Value',0); function varargout = slider1_Callback(h, eventdata, handles, varargin) global chon f = get(handles.slider1,'Value'); f = f*180 + 32; a = num2str(f); set(handles.edit2,'String',a); r = f + 495.7; c = (f - 32)*5/9; k = c + 273.15; if chon = = 1 t = r; elseif chon = = 2 t = c; elseif chon = = 3 t = k; end t = num2str(t); set(handles.edit3,'String',t); function varargout = checkbox1_Callback(h, eventdata, handles, varargin) global chon; off = [handles.checkbox2, handles.checkbox3]; mutual_exclude(off); chon = 1; function varargout = checkbox2_Callback(h, eventdata, handles, varargin) global chon; off = [handles.checkbox1, handles.checkbox3]; mutual_exclude(off); chon = 2; function varargout = checkbox3_Callback(h, eventdata, handles, varargin) global chon; off = [handles.checkbox2, handles.checkbox1]; mutual_exclude(off); chon = 3; ) GUI có dùng ho: Ta xây dng mt GUI dùng v th hm y=tsin(t). Giao din nh sau: 52
Ta dùng mt Axes, bn Pushbutton to nên giao din ny. Khi nhn Plot, th ca hm y = tsin(t) c v. Khi nhn Grid on, th c chia li. Khi nhn Grod off, li b xoá. Nhn Close óng th. Layout có thuc tính Name: ct1_59, HandleVisibility: callback. Các Pushbutton đều có thuộc tính FontSize: 12, FơntWeight: bold và các String là các tên lệnh. GUI được lưu với tên file là ct1_59.fig. Tiếp theo ta soạn thảo lệnh cho ct1_59.m: function varargout = ct1_59(varargin) if nargin = = 0 fig = openfig(mfilename,ʹreuseʹ); handles = guihandles(fig); guidata(fig, handles); if nargout > 0 varargout{1} = fig; end elseif ischar(varargin{1}) try [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard catch disp(lasterr); end end function varargout = pushbutton1_Callback(h, eventdata, handles, varargin) grid on function varargout = pushbutton2_Callback(h, eventdata, handles, varargin) grid off 53
function varargout = pushbutton3_Callback(h, eventdata, handles, varargin) close function varargout = pushbutton4_Callback(h, eventdata, handles, varargin) t = 0:0.01:20; y = t.*sin(t); plot(t,y); Tip theo ta xét mt GUI có giao din nh sau: Nhim v ca GUI l v th ca hm peaks theo các dng khác nhau( mesh, surf v contour) vi các Colormap khác nhau(hsv, hot, gray, prism, cool, winter v summer). Vic v các dng th thc hin nh các Pushbutton. Vic chn Colormap thc hin nh Listbox. Layout có thuc tính Name: ct1_60 v thuc tính HandleVisbility: on. Các Pushbutton u có thuc tính FontSize: 12 v FntWeight: bold. Ta lu GUI vi tên ct1_60.fig. Mã trong ct1_60.m gm: function varargout = ct1_60(varargin) if nargin = = 0 fig = openfig(mfilename,'reuse'); set(fig,'Color',get(0,'defaultUicontrolBackgroundColor')); handles = guihandles(fig); guidata(fig, handles); if nargout > 0 varargout{1} = fig; end elseif ischar(varargin{1}) try [varargout{1:nargout}] = feval(varargin{:}); catch 54
disp(lasterr); end end function varargout = pushbutton1_Callback(h, eventdata, handles, varargin) z = peaks(40); chon = get(handles.listbox1,'Value'); if chon = =1 colormap(hsv(256)); elseif chon = =2 colormap(hot(256)); elseif chon = =3 colormap(gray(256)); elseif chon = =4 colormap(prism(256)); elseif chon = =5 colormap(cool(256)); elseif chon = =6 colormap(winter(256)); elseif chon = =7 colormap(summer(256)); end mesh(z); function varargout = pushbutton2_Callback(h, eventdata, handles, varargin) z = peaks(40); chon = get(handles.listbox1,'Value'); if chon = =1 colormap(hsv(256)); elseif chon = =2 colormap(hot(256)); elseif chon = =3 colormap(gray(256)); elseif chon = =4 colormap(prism(256)); elseif chon = =5 colormap(cool(256)); elseif chon = =6 colormap(winter(256)); elseif chon = =7 colormap(summer(256)); end surf(z); function varargout = pushbutton3_Callback(h, eventdata, handles, varargin) z = peaks(40); chon = get(handles.listbox1,'Value'); if chon = =1 colormap(hsv(256)); 55
elseif chon = =2 colormap(hot(256)); elseif chon = =3 colormap(gray(256)); elseif chon = = 4 colormap(prism(256)); elseif chon = = 5 colormap(cool(256)); elseif chon = = 6 colormap(winter(256)); elseif chon = = 7 colormap(summer(256)); end contour(z); ) GUI có dùng đồ hoạ: Ta xây dựng một GUI dùng menu. Giao diện của GUI như sau: Menu Draw gồm các menu con Mesh, Contour và Close. GUI được lưu trong file ct1_61.fig và chương trình được lưu trong file ct1_61.m: function varargout = ct1_61(varargin) gui_Singleton = 1; gui_State = struct(ʹgui_Nameʹ, mfilename, ʹgui_Singletonʹ, gui_Singleton, ʹgui_OpeningFcnʹ, @ct1_61_OpeningFcn, ʹgui_OutputFcnʹ, @ct1_61_OutputFcn, ʹgui_LayoutFcnʹ, [] , 56
ʹgui_Callbackʹ, []); if nargin && ischar(varargin{1}) gui_State.gui_Callback = str2func(varargin{1}); end if nargout [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else gui_mainfcn(gui_State, varargin{:}); end handles.output = hObject; function varargout = ct1_61_OutputFcn(hObject, eventdata, handles) varargout{1} = handles.output; function mnumesh_Callback(hObject, eventdata, handles) z = peaks(40); mesh(z); function Untitled_3_Callback(hObject, eventdata, handles) z = peaks(40); contour(z); function mnuclose_Callback(hObject, eventdata, handles) clf close function mnudraw_Callback(hObject, eventdata, handles) 57
CHƯƠNG 2: MA TRẬN §1. MỘT SỐ KHÁI NIỆM ( Ma trận [A] gọi là đối xứng nếu [A]T = [A] ( Cho một ma trận vuông [A], cấp n. Ta nói ma trận [A] không suy biến (non singular) nếu ma trận có thể nghịch đảo được hay nói cách khác, định thức của ma trận khác không. ( Ma trận Hermite là một ma trận vuông có các phần tử là số phức bằng chuyển vị liên hợp của nó, nghĩa là phần tử ở hàng i cột j bằng số phức ∗ T liên hợp của phân tử ở hàng j cột i ⎣⎦⎡⎤AA= ⎣⎡ ⎦⎤ . Ví dụ ma trận ⎡⎤32j+ []A = ⎢⎥ là ma trận Hermite. ⎣⎦2j− 1 ( Ma trận Householder là một ma trận vuông dạng: 2 T [][]HE=− [][] UU [][]UUT Trong đó v là vec tơ cột khác zero ( Ma trận [A] gọi là trực giao nếu [A]T[A] = [E] ∗ T ( Ma trận phức [U] gọi là ma trận unita nếu ⎣⎡UU⎦⎣⎦⎤⎡⎤⎣⎦⎡⎤= E. Ví dụ ma ⎡⎤1j+−+ 1j ⎢⎥22 trận []U = ⎢⎥ là ma trận unita ⎢⎥1j+− 1j ⎣⎦⎢⎥22 ( Một ma trận chỉ có một cột gọi là một vec tơ ( Chuẩn của một vec tơ X, kí hiệu là X , là một số thực thoả mãn: ‐ X > 0 ‐ cX= c X ‐ XY+≤ X + Y Giả thiết X = [x1, x2, ,xn]T, ta thường dùng một trong 3 chuẩn sau đây: ‐ Xmaxx= j 1 j n = ‐ Xx2 ∑ j j1= 58
n 2 = ‐ Xx3 ∑ j j1= ( Chuẩn của một ma trận [A], kí hiệu là A, là một số thực thoả mãn: ‐ A > 0 ‐ cA= c A ‐ AB+≤ A + B ‐ AB≤ A B Ta thường dùng một trong 3 chuẩn sau đây: n ‐ Amaxa= i,j 1 i ∑ j1= n ‐ Amaxa= i,j 1 j ∑ i1= n 2 = ‐ Aa3 ∑ i,j i,j= 1 ( Ma trận [A] gọi là xác định dương nếu với vec tơ [x] bất kì ta có: [][][]xAx0T > ( Ma trận [A] gọi là nửa xác định dương nếu với vec tơ [x] bất kì ta có: [][][]xAx0T ≥ Ta định nghĩa ma trận xác định âm và nửa xác định âm một cách tương tự. ( Hạng của ma trận là cấp của ma trận con của ma trận ấy có định thức khác không còn mọi ma trận con cấp cao hơn đều có định thưc bằng không(ma trận con là ma trận có được bằng cách xoá một số hàng và cột của ma trận ban đầu). §2. BIẾN ĐỔI HOUSEHOLDER 1. Ma trận Householder: Ta biến đổi ma trận [A] về dạng có các phần tử thuộc đường chéo chính, các phần tử phía trên và phía dưới đường chéo chính khác zero, còn các phần tử còn lại bằng zero(ma trận ba đường chéo) bằng cách dùng phép biến đổi Householder. Phép biến đổi Householder dùng ma trận Householder. [][]UUT [][]HE=− (1) Q 59
Trong đó: 11T 2 QUUU==[][] [] (2) 22 Do [H] đối xứng nên: TT T ⎛⎞⎛⎞[][]UU[][] UU [][][][][]HH==− HH⎜⎟⎜⎟ E[] E − ⎝⎠⎝⎠QQ TT [][]UUT [][][]UUU( )[] U =−[]E2 + QQ2 [][]UUT [U2QU]()[ ]T =−[]E2 + =[] E QQ2 Từ đây ta thấy [H] cũng là ma trận trực giao. Cho [X] là vec tơ bất kỳ và khảo sát phép biến đổi [H][X]. Chọn: [U] = [X] + k[I1] (3) Trong đó: T kX=±[] []I1001 = ⎣⎡ L ⎦⎤ Ta có: T ⎧⎫T ⎛⎞[][]UU ⎪⎪[][]UX( + kI[]1 ) [][]HX=−⎜⎟ [] E[] X =−⎨⎬ [] E[] X ⎝⎠QQ⎩⎭⎪⎪ UXXT + kIXT 2 [][][]( []1 []) []Uk()+ kX[]1 =−[]XX =−[] QQ Nhưng: T 2 T TT2 2Q=+([] X k[] I11) ([] X + k[] I) =[] X + k([] X[][] I 1111 + I[] X) + k[][] I I 222 =+k 2kx11 += k 2(k + kx ) Như vậy: T [][][][]HX=−=−=− X U kI[]1 ⎣⎡ k 0 0L 0⎦⎤ (4) nghĩa là phép biến đổi loại trừ tất cả các phần tử của [X] trừ phần tử đầu tiên. 2. Biến đổi Householder một ma trận đối xứng: Bây giờ ta áp dụng phép biến đổi cho ma trận [A] đối xứng: T TT ⎡⎤10[] ⎡⎤⎡⎤aX11[] a 11 [] X ⎣⎦⎡⎤PA1 []==⎢⎥⎢⎥⎢⎥ (5) ⎣⎦[]0H [ ]⎣⎦⎣⎦[]X[] A′ [][][] HX HA[]′ 60
Trong đó [X] là cột đầu tiên của [A] với phần tử đầu tiên bị bỏ đi. [A’] có được từ [A] bằng cách bỏ đi cột và hàng đầu tiên. Ma trận [H] cấp (n ‐1) được xây dựng theo các công thức (1) ÷ (3). Do (4) ta thấy phép biến đổi này làm cột đầu tiên của [A] trở thành: ⎡⎤a11 ⎢⎥−k ⎢⎥ ⎡⎤a11 ⎢⎥= ⎢⎥0 ⎣⎦[][]HH ⎢⎥ ⎢⎥M ⎣⎦⎢⎥0 Phép biến đổi: ⎡⎤T aHX11 ()[][] ⎣⎦⎡⎤PAP11[] ⎣⎦⎡⎤=→⎢⎥[] A (6) ⎣⎦⎢⎥[][][]HX HA[]′ [] H sẽ đường chéo hoá hàng đầu tiên và cột đầu tiên của ma trận [A]. Sơ đồ biến đổi của ma trận 4×4 là: 1 0 0 0 a11 a12 a13 a14 1 0 0 0 a11 ‐k 0 0 0 a21 0 ‐k [Q][A’] × × = 0 [Q] a31 [A’] 0 [Q] 0 [Q] 0 a41 0 0 Hàng và cột thứ 2 của ma trận [A] được biến đổi tiếp bằng cách dùng phép biến đổi đối với phần bên phải, phía dưới của ma trận. Phép biến đổi này có thể biểu diễn bằng []PAP22[][]→ [ A,] trong đó: T ⎡⎤[]E02 [] []P2 = ⎢⎥ (7) ⎣⎦[]0H [ ] với [E2] là ma trận đơn vị 2×2 và [H] là ma trận (n ‐ 2)×(n ‐ 2) có được bằng cách chọn [X] từ (n ‐ 2) phần tử phía dưới của cột thứ 2 của ma trận [A]. Thực hiện (n ‐ 2) phép biến đổi: T ⎡⎤[]E0i [] []Pi = ⎢⎥ i = 1, 2, , n ‐ 2 ⎣⎦[]0H [ ] để có được ma trận ba đường chéo(tridiagonal). Ta có: 61
T ⎛⎞[][]UU [AU′][] TT []AH′′[]=[] A⎜⎟[] E − =−[] A ′[] U =−[] A ′[][] VU ⎝⎠QQ Trong đó: []AU′ [] []V = (8) Q Do vậy: T ⎛⎞[][]UU T []HA[]′′[] H=−⎜⎟ [] E()[] A −[][] VU ⎝⎠Q T TT[UU][ ] =−[]AVU′′[][] −()[] AVU −[][] Q TTT′ T [][]UUA( [ ]) [][][] UUVU( )[] =−[]AVU′ [][] − + QQ =−−+[]AVUUV2gUU′ [][ ]TT[ ][ ] [ ][ ] T Trong đó: [][]UVT g = (9) 2Q Đặt: [W] = [V] ‐ g[U] (10) Ta thấy ngay phép biến đổi có dạng: []HA[]′′[] H=−[] A[] WU[ ]TT −[ UW][ ] (11) Thuật toán có thể tóm lại như sau: ‐ Cho [A’] là ma trận vuông cấp (n ‐ i) có được từ phần dưới bên phải của ma trận [A] T ‐ Đặt ⎣⎦⎡⎤Xa= ⎣⎦⎡⎤i++ 1,i a i 2,iL a n,i ‐ Tính []X. Cho k = []X nếu x1 > 0 và k = ‐ [X] nếu x1 < 0 T ‐ Cho ⎣⎦⎡⎤Ukxx=+⎣⎦⎡⎤12L x ni− 2 []U ‐ Tính Q = 2 []AU′ [] ‐ Tính []V = Q [][]UVT ‐ Tính g = 2Q 62
‐ Tính [W] = [V] ‐ g[U] ‐ Tính []AA=−[]′ [][] WUUWTT −[ ][ ] ‐ Đặt aai,i++ 1==− i 1,i k Ta xây dựng hàm housetrans() để thực hiện thuật toán trên: function A = housetrans(A) % Bien doi Householder ma tran A thanh ma tran % ba đường chéo dang[c\d\c]. % De co c va d dung d = diag(A), c = diag(A,1). n = size(A, 1); for k = 1:n‐2 u = A(k+1:n, k); uMag = sqrt(dot(u, u)); if u(1) < 0; uMag = ‐uMag; end u(1) = u(1) + uMag; A(k+1:n, k) = u; % Luu u vao phan duoi cua A. H = dot(u, u)/2; v = A(k+1:n,k+1:n)*u/H; g = dot(u, v)/(2*H); v = v ‐ g*u; A(k+1:n, k+1:n) = A(k+1:n, k+1:n) ‐ v*uʹ ‐ u*vʹ; A(k, k+1) = ‐uMag; end k = zeros(n); for i = 1:n k(i, i) = A(i, i); end for i = 1:n‐1 k(i, i+1) = A(i, i+1); k(i+1, i) = A(i, i+1); end A = k; 63
Để tính ma trận ba đường chéo theo phép biến đổi Householder ta dùng chương trình cthousetrans.m: clear all, clc a = [ 1 2 3 4; 2 9 3 5; 3 3 3 7; 4 5 7 6]; b = householder(a) d = diag(b) c = diag(b, 1) §3. BIẾN ĐỔI THÀNH MA TRẬN HESSENBERG Nếu ma trận [A] là ma trận đối xứng, phương pháp Householder có thể được sử dụng để biến đổi nó thành ma trận đồng dạng đối xứng ba đường chéo. Nếu ma trận [A] không đối xứng, phương pháp Householder biến đổi ma trận [A] thành ma trận đồng dạng Hessenberg. Ma trận Hessenberg là ma trận có dạng: ⎡⎤aaa11 12 13L a 1,n ⎢⎥aaaL a ⎢⎥21 22 23 2n []H = ⎢⎥0a32 a 33L a 2n ⎢⎥ ⎢⎥MMML M ⎢⎥ ⎣⎦000L ann Ta thực hiện phép biến đổi Householder trên ma trận [A] và có được: [Q][H][Q’] = [A] trong đó [Q] là ma trận trực giao (ta gọi đây là phân tích Hessenberg ma trận [A]) . Thuật toán có thể tóm lại như sau: ‐ Cho [Q] là ma trận đơn vị cấp n T ‐ Đặt ⎣⎦⎡⎤X0a= ⎣⎦⎡⎤i2,i+ L a n,i ‐ Tính []X. Cho α= []X nếu ai+2,i > 0 và α = ‐ [X] nếu ai+2,i < 0 T ‐ Cho ⎣⎦⎡⎤U0=α+⎣⎦⎡⎤ x2niL x− 2 []U ‐ Tính β= 2 []UU[]′ ‐ Tính []PE=− [] β 64
‐ Tính []QQP= [][] ‐ Tính [APAP] = [][][] Ta xây dựng hàm hessenberg() để thực hiện phép phân tích trên: function [H, Q] = hessenberg(a) [n, n] = size(a); q = eye(n); for k = 1:n ‐ 2 alfa = 0; for j = k+1:n alfa = alfa + a(j, k)^2; end alfa = sign(a(k+1, k))*sqrt(alfa); u = zeros(1, n); u(k+1:n) = a(k+1:n, k); u(k+1) = u(k+1) + alfa; beta = .5*u*uʹ; p = eye(n); for i = 1:n p(i, 1:n) = p(i, 1:n) ‐ (u(i)*u(1:n))/beta; end q = q*p; a = p*a*p; end H = a; Q = q; Để phân tích ma trận ta dùng chương trình cthessenberg.m: clear all, clc a = [ 1 2 3 4; 5 6 7 4; 6 4 8 9; 3 5 7 9]; [H, Q] = hessenberg(a) §4. PHÂN TÍCH MA TRẬN THEO PHƯƠNG PHÁP DOOLITTLE 65
Một ma trận không suy biến [A] gọi là phân tích được thành tích hai ma trận [L] và [R] nếu: [A] = [L] [R] Việc phân tích này, nếu tồn tại, là không duy nhất. Nếu ma trận [L] có các phần tử nằm trên đường chéo chính bằng 1, ta có phép phân tích Doolittle. Nếu ma trận [R] có các phần tử nằm trên đường chéo chính bằng 1, ta có phép phân tích Crout. Nếu [R] = [L]T (hay [L] = [R]T) ta có phép phân tích Choleski. Với ma trận bậc 3, [L] và [R] có dạng: ⎡⎤100 ⎡ rrr11 12 13 ⎤ ==⎢⎥ ⎢ ⎥ []Ll⎢⎥21 10 [] R ⎢ 0rr 22 23 ⎥ ⎢⎥ ⎢ ⎥ ⎣⎦ll131 32 ⎣ 00r 33 ⎦ Để tìm lij và rij ta thực hiện phép nhân. Sau khi nhân ta có: ⎡⎤rr11 12 r 13 =++⎢⎥ []Arlrlr⎢⎥11 21 12 21 22 rlr 13 21 23 ⎢⎥ ⎣⎦rl11 31 rl 12 31+++ rl 22 32 rl 13 31 rl 23 32 r 33 Bây giờ ta thực hiện phép khử Gauss đối với phương trình trên. Đầu tiên ta chọn hàng thứ nhất làm trụ và thực hiên phép biến đổi: hàng 2 ‐ l21 × hàng 1 (khử a21) → hàng 2 hàng 3 ‐ l31 × hàng 1 (khử a31) → hàng 3 kết quả ta có: ⎡⎤rr11 12 r 13 = ⎢⎥ []A0rr12223⎢⎥ ⎢⎥ ⎣⎦0rlrlr22 32 23 32+ 33 Sau đó ta lấy hàng thứ hai làm trụ và thực hiện biến đổi: hàng 3 ‐ l32 × hàng 2 (khử a32) → hàng 3 và có: ⎡⎤rrr11 12 13 = ⎢⎥ []A0rr22223⎢⎥ ⎢⎥ ⎣⎦00r33 Như vậy ta thấy ngay rằng ma trận [R] là ma trận có được khi thực hiện loại trừ Gauss tiến ma trận [A] và các phần tử của [L] là các nhân tử dùng khi 66
loại trừ aij. Điều đó có nghĩa là để tìm ma trận [L] và [R] ta dùng phép khử Gauss tiến. Ta xây dựng hàm doolittle() để thực hiện loại phân tích Doolittle. function [l,r] = doolittle(A) %Phan tich ma tran A thanh A = L*U n = size(A, 1); u = zeros(n); for k = 1:n‐1 for i = k+1:n if A(i, k)~= 0.0 lambda = A(i, k)/A(k, k); A(i, k+1:n) = A(i, k+1:n) ‐ lambda*A(k, k+1:n); A(i, k) = lambda; end end end l = tril(A); for i = 1:n l(i, i) = 1; end l = triu(A); for i = 1:n l(i,i) = A(i, i); end §5. PHÂN TÍCH MA TRẬN THEO PHƯƠNG PHÁP CROUT Tương tự như thuật toán Doolittle, ta có thể phân tích ma trận [A] theo thuật toán Crout thành tích của ma trận [L] và [R]. Các ma trận bậc 3 theo Crout có dạng: ⎡⎤l0011 ⎡⎤ 1rr 12 13 ==⎢⎥ ⎢⎥ []Lll0⎢⎥21 22 [] R01r ⎢⎥ 23 ⎢⎥ ⎢⎥ ⎣⎦lll31 32 33 ⎣⎦ 001 Để tìm lij và rij ta thực hiện phép nhân. Sau khi nhân ta có: 67
⎡⎤llr11 11 12 lr 11 13 =++⎢⎥ []Allrllrlr⎢⎥21 21 12 22 21 13 22 23 ⎢⎥ ⎣⎦llrllrlrl31 31 12+++ 32 31 13 32 23 33 Như vậy: a11 = 1. r11 + 0.0 + 0.0 = r11 ; a12 = r12 ; a13 = r13 a21 = l21r11 ; a22 = l21r12 + r22 ; a23 = l31r11 a31 = l31r11 ; a32 = l31r12 ; a33 = l31r13 + l32r23 + r33 Một cách tổng quát ta có : với j > i : lij = rji = 0 với i = 1 : r1j = a1j (j = 1 tới n) lj1 = aj1/r11 (j = 1 tới n) với i = 2 tới n i−1 rij = aij − ∑likrkj ( j = i tới n) k=1 i−1 a ji − ∑l jkrki k=1 l ji = (j = i tới n) rii Ta xây dựng hàm crout() để phân tích ma trận theo thuật toán Crout: function [l, r] = crout(a) n = size(a, 1); l = zeros(n); r = zeros(n); for i = 1:n r(1, i) = a(1, i); l(i, i) = 1.; l(i, 1) = a(i, 1)/a(1, 1); end for k = 2:n r(k, k:n) = a(k, k:n) ‐ l(k, 1:k)*r(1:k, k:n); if k~= n 68
for i = 1:n l(i, k) = (a(i, k)‐ l(i, 1:k‐1)*r(1:k‐1, k))/r(k, k); end end end §6. PHÂN TÍCH MA TRẬN THEO PHƯƠNG PHÁP CHOLESKI Thuật toán Choleski cho phép phân tích ma trận [A] thành tích hai ma trận: [A] = [L][L]T. Thuật toán này đòi hỏi: ‐ [A] là ma trận thực, đối xứng ‐ [A] là ma trận xác định dương Ta vuông [A] cấp 3 theo thuật toán Choleski: ⎡⎤⎡⎤⎡⎤aaa11 12 13 l 11 00lll 11 21 31 ⎢⎥⎢⎥⎢⎥= ⎢⎥⎢⎥⎢⎥aaa21 22 23 ll 21 22 00ll 22 32 ⎢⎥⎢⎥⎢⎥ ⎣⎦⎣⎦⎣⎦aaa31 32 33 lll 31 32 33 00l 33 Sau khi thực hiện phép nhân ta có: 2 ⎡⎤a11 a 12 a 13⎡⎤ l 11 ll 11 21 ll 11 31 ⎢⎥ ⎢⎥=+22 + ⎢⎥aaa21 22 23⎢⎥ llll 11 21 21 22 llll 21 31 22 32 ⎢⎥⎢⎥222 ⎣⎦aaa31 32 33⎣⎦ lllllllll 11 31 21 31+++ 22 32 31 32 33 Vế phải là ma trận đối xứng. Cân bằng các phần tử của hai ma trận ta có: l11== a 11 la/l 21 21 11 la/l 31 = 31 11 222 lall(all)/llall22=− 22 21 32 =− 32 21 31 22 33 =−− 33 31 32 Tổng quát, với ma trận cấp n, ta có: j T [][]LL=++⋅⋅⋅+= lli1 j1 ll i2 j2 ll ik jk i ≥ j ()ij ∑ k1= Cân bằng với phần tử của ma trận [A] ta có: j aij==+=∑ l ik l jk i j,j 1, ,n j 1,2, ,n k1= Do ma trận [L] là ma trận tam giác trái nên đối với cột thứ nhất ta có: la11== 11 la/l i1 i1 11 Đối với cột khác, rút lij ra khỏi tổng ta có: 69
j−1 allllij=+∑ ik jk ij jj k1= Nếu i = j (phần tử trên đường chéo) thì: j1− 2 ljj=− a jj∑ l jk j = 2,3, ,n k1= và phần tử nằm ngoài đường chéo: ⎛⎞j1− 1 lij=−⎜⎟ a ij∑ l ik l jk j = 2, 3, , n i =++ j 2, j 3, ,n ⎝⎠k1= ljj Dựa vào thuật toán trên ta xây dựng hàm choleski() function L = choleski(A) % Phan tich ma tran a thanh A = LL’. % Cu phap: L = choleski(A) f = posdef(A); if f == 0 error(ʹMa tran khong xac dinh duong!ʹ); return end n = size(A, 1); for j = 1:n temp = A(j, j) ‐ dot(A(j, 1:j‐1),A(j, 1:j‐1)); if temp < 0.0 error(ʹMa tran khong xac dinh duongʹ) end A(j, j) = sqrt(temp); for i = j+1:n A(i, j)=(A(i, j) ‐ dot(A(i, 1:j‐1),A(j, 1:j‐1)))/A(j, j); end end L = tril(A); function f = posdef(M) %Kiem tra lieu ma tran M co xac dinh duong hay kong isposdef = true; 70
for i=1:length(M) if ( det( M(1:i, 1:i) ) <= 0 ) isposdef = false; break; end end f = isposdef;% 0 neu sai, 1 neu dung §7. PHÂN TÍCH QR BẰNG THUẬT TOÁN HOUSEHOLDER Cho ma trận [A], phân tích QR của nó cho ta: [A] = [Q]*[R] Trong đó [Q] là ma trận trực giao và [R] là ma trận tam giác phải. Ta dùng biến đổi Householder để tìm các ma trận [Q] và [R]. [][][]HHn1−− n2⋅⋅⋅ HAR 1 [] =[ ] (1) Như vậy: −1 −−11 []AHH=([][][]−− ⋅⋅⋅ H) [] RH =[] ⋅⋅⋅ [ H − ][ HR − ][] n1 n2 1 1 n2 n1 (2) =[][HHHRQR1n2n1 ⋅⋅⋅−− ][][][][] = Tích của tất cả các ma trận Householder: []QH= [][1n2n1L H−− ][] H (3) không những đối xứng mà còn trực giao như mỗi ma trận [Hk]: T T [][]QQ=⋅⋅⋅([][ H H−− ][] H) [][ H ⋅⋅⋅ H −− ][] H 1n2n11n2n1 TT T =[HHn1−− ][] n2 ⋅⋅⋅ [][][ HHHH 1 1 ⋅⋅⋅ n2 −− ][] n1 =[] E Ta xây dựng hàm qrdecom() để phân tích ma trận: function [Q, R] = qrdecom(A) %Phan tich QR n = size(A, 1); R = A; Q = eye(n); for k = 1:n ‐ 1 H = householder(R(:, k), k); R = H*R; %Pt.(1) Q = Q*H; %Pt.(3) 71
end Hàm householder() dùng để tạo ra ma trận Householder: function H = householder(x, k) % Tao ma tran Householder n = length(x); tmp = sum(x(k+1:n).^2); g = sqrt(x(k)^2 + tmp); c = sqrt((x(k) + g)^2 + tmp); u = zeros(n, 1); u(k) = (x(k) + g)/c; u(k + 1:n) = x(k + 1:n)/c; H = eye(n) ‐ 2*u*uʹ; %ma tran Householder Để phân tích ma trận ta dùng chương trình ctqrdecom.m: clear all, clc a = [4 1 3 ‐2; 1 ‐2 4 1; 3 4 1 2; ‐2 1 2 3]; [q, r] = qrdecom(a) §8. PHÂN TÍCH QR BẰNG THUẬT TOÁN QUAY GIVENS Kỹ thuật quay Givens là một phương pháp để phân tích ma trận [A] thành tích của ma trận [Q] và ma trận [R] bằng cách làm cho các phần tử lần lượt bằng zero cho đến khi có được ma trận tam giác phải. Ý tưởng là dùng một ma trận quay đơn giản 2 × 2 đặt dọc theo đường chéo chính của một ma trận đơn vị và làm cho một phần tử của ma trận bằng zero. Các phần tử của ma trận quay để quay một vec tơ ngược chiều kim đồng hồ một góc θ là: ⎡⎤cosθ− sin θ []Qθ = ⎢⎥ ⎣⎦sinθθ cos Nếu ta muốn quay vec tơ [x1 x2]T và muốn làm cho x2 bằng zero rồi quaytheo chiều kim đồng hồ một góc θ(hay ngược chiều kim đồng hồ một góc ‐θ) trong đó: 72
x θ=arctg 2 x1 thì ma trận quay để thực hiện phép quay này theo chiều kim đồng hồ một góc θ là: ⎡⎤cosθθ sin []Qθ = ⎢⎥ ⎣⎦−θsin cos θ Trong đó: x x cosθ= c = 1 sinθ= s = 2 22 22 xx12+ xx12+ Do đó: 1 ⎡⎤⎡⎤xx12 cs []Q == θ 22⎢⎥⎢⎥−−xx sc xx12+ ⎣⎦⎣⎦21 Chú ý là như mong muốn: ⎡⎤xx22+ 12 22 ⎡⎤⎡xcxsx112+ ⎤⎢⎥⎡ xx+ ⎤ []Q ===xx22+ ⎢ 12⎥ θ ⎢⎥⎢xsxcx−+ ⎥⎢⎥12 ⎣⎦⎣212 ⎦⎢⎥⎣⎢ 0 ⎦⎥ ⎣⎦0 Nếu A là ma trận m × n, ta sẽ xem điều gì xảy ra khi ta thay các phần tử của [Q] vào ma trận con xác định bằng các cột và hàng thứ i, các cột và hàng thứ j. Nói cách khác ta thay ma trận 2 × 2 này dọc theo đường chéo chính tại một số điểm: ⎡1LLL 000⎤ ⎢ ⎥ ⎢ MO M M MO M⎥ ⎧δkl ki,lj≠≠ ⎪ ⎢0cs0LLL⎥ ⎪ c k, l== i; k,l j ⎢ ⎥ == []Gkl ⎨ ⎢ MM MOMMM⎥ ⎪ s ki;lj== ⎢0sc0LLL− ⎥ ⎩⎪−s kj;li== ⎢ ⎥ ⎢ MM0 MMOM⎥ ⎣⎢0001LLL⎦⎥ Như vậy [G] là ma trận đơn vị m × m ngoại trừ các giá trị đã bị thay thế: gii = gjj = c gij = ‐gij = s Điều này sẽ tạo ra ma trận unita: [G]T[G] = [E] nghĩa là: 73
∑gglk lp=δ kp l và đòi hỏi: c2 + s2 = 1 Điều này đúng vì cos2θ + sin2θ = 1 ∀θ. Khi ma trận này được áp dụng cho ma trận m × n ta có: ⎧ ⎪∑δ=klaaki,j lp kp ≠ ⎪ l bkp==∑∑ga kl lp⎨ ga il lp =+= ca ip sa jp k i ll⎪ ⎪ ∑gajl lp= −+ sa ip ca jp k = j ⎩⎪ l Như vậy ma trận mới chỉ bị thay đổi ở hàng i và cột j. Ta chọn s và c sao cho các phần tử ở cột r và hàng j bằng zero: a jr a s = c = ir 22 22 aajrir+ aajrir+ Như vậy ta sẽ có: −+aa ab b ==jr ir ir jr 0 jr 22 aajr+ ir Ta xây dựng hàm givens() để thực hiện thuật toán trên: function [Q, R] = givens(A); % Phan tich QR bang thuat toan quay Givens n = size(A, 1); Q = eye(n); for j = 1:n‐1 for i = n:‐1:j+1 z = 1/sqrt(A(i‐1, j)^2 + A(i, j)^2); c = A(i‐1, j)*z; s = A(i, j)*z; A(i‐1:i,:) = [c s; ‐s c]*A(i‐1:i,:); Q(i‐1:i,:) = [c s; ‐s c]*Q(i‐1: i,:); end end R = A; 74
Q = Qʹ; Để phân tích một ma trận ta dùng chương trình ctgivens.m: clear all, clc A = [17 24 30 17; 8 13 20 7; 2 10 8 6; ‐23 ‐43 ‐54 ‐26]; [Q, R] = givens(A) §9. PHÂN TÍCH QR BẰNG THUẬT TOÁN GRAM ‐ SCHMIDT Ta có thể thực hiện việc phân tích ma trận [A] thành tích các ma trận [Q] và [R] bằng cách trực giao hoá các cột của ma trận [A]. Ta gọi các cột của ma trận [A] là a1, ,an. Từ các vec tơ này ta muốn có n vec tơ trực giao v1, ,vn. Vec tơ trực giao đầu tiên được chọn là: va11= Để có vec tơ thứ hai, ta dùng y2 nhưng trừ bớt đi phần y2 cùng chiều với v2. Như vậy ta có: vyba21=− 1 với b được chọn sao cho v1 trực giao với v2: vv12=−=−= v(a 1 2 bv) 1 va 12 bvv 11 0 hay: va b = 12 vv11 Tiếp tục quá trình đến bước thứ k ta có: k1− va va=− ik v kk∑ i i1= vvii Như vậy thuật toán gồm các bước: a1 ‐ ra,q11== 1 1 r11 - lặp từ k = 2 đến n ⎛⎞k1− qakk=−⎜⎟∑ rqr ikikk ⎝⎠i1= với T rqaik= i k và rkk được chọn sao cho q1k = , nghĩa là: 75
za=−kkik qr rzkk = Ta xây dựng hàm qrgramschmidt() để thực hiện thuật toán trên: function [Q, R] = qrgramschmidt(A); % Phan tich mt bang thuat toan Gram ‐ Schmidt [m,n] = size(A); R(1,1) = norm(A(:, 1)); Q(:,1) =A(:, 1)/R(1, 1); for k = 2:n R(1:k‐1, k) = Q(1:m, 1:k‐1)ʹ*A(1:m, k); z = A(1:m, k) ‐ Q(1:m, 1:k‐1)*R(1:k‐1, k); R(k,k) = norm(z); Q(1:m,k) = z/R(k, k); end Để phân tích một ma trận ta dùng chương trình chương trình ctqrgamschmidt.m: clear all, clc a = [ 1 2 3 4 5; 6 7 8 9 0; 3 4 5 6 7; 8 9 0 1 2; 2 4 6 8 1]; [q, r] = qrgramschmidt(a) §10. PHÂN TÍCH MA TRẬN THEO GIÁ TRỊ RIÊNG Cho ma trận [A], ta có: [A][X] = λ[X] Nếu ta đặt [U] là một ma trận mà các cột của nó là các vec tơ riêng của ma trận [A] và ma trận [Λ] là ma trận đường chéo có các phần tử trên đường chéo chính là λi thì: [A][U] = [Λ][U] hay: [A] = [U][Λ][U]‐1 Dạng này của ma trận được gọi là dạng phân tích theo giá trị riêng và vec tơ riêng. Ta dùng chương trình cteigdecom.m để phân tích ma trận: clear all, clc 76
a = [ 1 3 5; 3 4 9; 5 9 6]; [L, U] = eigjacobi(a) §11. PHÂN TÍCH LQ Cho ma trận [A]T, ta có thể phân tích QR ma trận này thành: [A]T = [Q1][R1] Do ([Q][R])T = [R1]T[Q1]T nên: ([A]T)T = [A] = [L][Q] và ta nhận được phân tích LQ của ma trận [A]. Ta xây dựng hàm lqdecom() để thực hiện thuật toán này: function [Q, L] = lqdecom(A) A = Aʹ; [Q, L] = qrdecom(A); L = Lʹ; Q = Qʹ; Để phân tích một ma trận ta dùng chương trình ctlqdecom.m: clear all, clc a = [ 1 3 5; 2 4 6; 7 8 9]; [Q, L] = lqdecom(a) §12. PHÂN TÍCH JORDAN 1. Ma trận có thể đường chéo hoá: Ma trận [A] gọi là có thể đường chéo hoá nếu và chỉ nếu tồn tại phép biến đổi đồng dạng [V] sao cho [A] = [V][Λ][V]‐1 trong đó [Λ] là ma trận đường chéo [Λ] = diag(λ1, λ2, , λn). Điều kiện cần để [A] có thể đường chéo hoá là [A] có n vec tơ riêng độc lập tuyến tính. Điều kiện đủ để [A] có thể đường chéo hoá là [A] có n giá trị riêng phân biệt vì khi [A] có n giá trị riêng phân biệt thì các vec tơ riêng tương ứng là độc lập tuyến tính. Số lần lặp lại mi của giá trị riêng λi gọi là vô số đại số (algebraic multiplicity) của λi, kí hiệu là AM(λi ). Số vec tơ riêng độc lập tuyến tính tương ứng với giá trị riêng λi gọi là vô số hình học (geometric multiplicity) của λi, kí hiệu là GM(λi ). 77
2. Dạng Jordan: Khi không thể tìm được n giá trị riêng phân biệt, nghĩa là ma trận [A] không có n vec tơ riêng độc lập tuyến tính thì ma trận [A] không thể đường chéo hoá. Tuy nhiên, nếu có phép biến đổi đồng dạng [M] biến đổi [A] thành [J]: [A] = [M][J][M]‐1 Trong đó [J] là ma trận gần đường chéo: [J] = diag(J1, , Jn) ⎡⎤λi 10L 0 ⎢⎥01λ ⎢⎥i OM []Ji = ⎢⎥MOλi OM ⎢⎥ λ ⎢⎥MOO i 1 ⎢⎥ ⎣⎦MLL0 λi là khối Jordan và ma trận [J] được gọi là dạng Jordan kinh điển của ma trận [A]. Số khối Jordan bằng số vec tơ riêng độc lập tuyến tính của ma trận [A], nghĩa là bằng GM(λi). Cụ thể, mỗi vec tơ riêng độc lập tuyến tính tương ứng với mỗi khối. Do vậy nếu ma trận [A] có các vec tơ riêng độc lập tuyến tính thì dạng Jordan trùng với dạng đường chéo của ma trận [S]‐1[A][S] = [Λ] trong đó [Λ] = diag(λ1, , λn) và [S] có các cột là các vec tơ riêng của [A] 3. Xây dựng dạng Jordan của ma trận [A]: Khi [A] không có n vec tơ riêng độc lập tuyến tính để tạo ra các cột của ma trận [M] thì ta có thể thêm các vec tơ độc lập tuyến tính vào các vec tơ riêng để tạo ra ma trận này. Trước hết ta khảo sát một giá trị riêng λi có GM(λi) < AM(λi). Nếu GM(λi) = pi , AM(λi) = mi thì ta cần tìm mi ‐ pi vec tơ độc lập tuyến tính để kết hợp với giá trị riêng này. Các vec tơ này được tạo từ các vec tơ riêng và được gọi là vec tơ riêng tổng quát hoá của [A]. Gọi λ là giá trị riêng và [x] là vec tơ riêng tương ứng. k ‐ 1 vec tơ riêng tổng quát hoá {[x1], , [xk]} được tạo ra như sau: []Ax[]11=λ [] x []Ax[]=λ [][] x + x 221 M []Ax[]kkk1=λ [][ x + x− ] {[x1], , [xk]} tạo thành chuỗi các vec tơ có vec tơ [x1] đứng đầu. Chuỗi này tương ứng với khối Jordan đơn. 78
⎡ λ 10L 0⎤ ⎢ 010λ ⎥ ⎢ M ⎥ []Ax,,x[][]1n1nKK= x,,x⎢ MOOO0 ⎥ ⎢ ⎥ ⎢ 01LL λ ⎥ ⎣⎢ 00LL λ ⎦⎥ Xuất phát từ vec tơ tổng quát hoá bậc k của [A] ứng với λ, kí hiệu là [xk] ta có: [xk] [xk‐1] = ([A] ‐ λ[E])[xk] M [xi] = ([A] ‐ λ[E])k‐i[xk] M [x1] = ([A] ‐ λ[E])k‐1[xk] Chú ý là [x1] là một vec tơ riêng của [A] vì: ([A] ‐ λ[x1]) = ([A] ‐ λ[E])([A] ‐ λ[E])k‐1[xk] Để phân tích ma trận [A] ta dùng thuật toán Filipov gồm các bước sau: ‐ Giả sử rằng kích thước không gian cột của ma trận [A] là r < n. Phải có r vec tơ độc lập tuyến tính [xi] trong không gian cột mà nó là các vec tơ riêng hay vec tơ riêng tổng quát hoá, nghĩa là [A][xi] = λ[xi] hay [A][xi] = λ[xi] + [xi‐1] ‐ Giả sử rằng không gian không và không gian cột của ma trận [A] có phần chung với kích thước p. Mỗi vec tơ [xi] trong jg guan không của [A] là một vec tơ riêng tương ứng với λ = 0, nhưvậy [A][xi] = 0. Bây giờ nếu [xi] cũng là không gian cột của [A] thì [xi] = [A][yi] với mọi [yi] ‐ Cuối cùng do kích thước của không gian không của [A] là n ‐ r và p của các vec tơ là trong cả không gian không lẫn không gian cột nên có n ‐ r ‐ p vec tơ [zi] ở trong không gian không mà không ở trong không gian cột. Các vec tơ [xi], [yi] và [zi] tìm được là độc lập tuyến tính. Chúng tạo nên các cột của [M] và [J] = [M][A][M]‐1 là dạng Jordan. Ta xây dựng hàm jordandecom() thực hiện thuật toán trên: function [M, J] = jordandecom(a) %Tinh phan tich Jordan cua ma tran A % sao cho A*M = M*J small = 2*sqrt(eps); [r, c] = size(a); 79
if r ~= c error(ʹMa tran A phai la ma tran vuong!ʹ) end n = r; if n == 1 J = a; M = 1; return end if n<1 J = []; M = []; return end [m, d] = eig(hess(a)); d = sort(diag(d)); tiny = norm(a)*eps; %lam cac gia tri bang zero p = find(abs(d)<=tiny); if ~isempty(p) d(p) = 0; end %A*M=M*J [M, J] = jord(a, d, small); function [M, D] = jord(a, d, small) %Tinh phan tich Jordan cua ma tran A norma = sqrt(mean(mean(abs(a)))); tiny = norma*eps; if nargin<3 small = (1 + norma)*sqrt(eps); end [r, c] = size(a); if r~=c error(ʹA phai la ma tran vuong!ʹ) 80
end n = r; I = eye(n); if r == 1 D = a; M = 1; return end if r 1e6 Condition_number_of_A = condofa warning(ʹSo dieu kien cua A qua lon!ʹ); end d = d(:); if size(d,1) ~= n d = d error(ʹGia tri rieng khong dung!ʹ) end da = det(a); dp = prod(d); e = abs(abs(da) ‐ abs(dp)); if e>sqrt(eps) disp(ʹ ʹ) warning(ʹCac gia tri rieng co the khong chinh xac!ʹ) end ds = flipud(sort(d)); sds = size(ds,1); du = flipud(unique(ds)); sdu = size(du, 1); if sdu == sds 81
[M, D] = eig(a); return end M = []; for kk = 1:sdu e = du(kk); ameig = sum(ismember(ds, e)); a1 = a ‐ e*I; if ameig == 1 [u, s, v] = svd(a1); M =[M v(:, end)]; else pp = 0; ns = []; pp = pp + 1; aa = I; for k = 1:ameig aa = a1*aa; nn = size(nulld(aa, small), 2); ns(k) = nn; end nsaa = [0; ns(:)]ʹ; dns = diff(nsaa); if max(dns) ~= dns(1) Cond_of_A = cond(a) save jord M = I; D = I; error(ʹKich thuoc khong gian khong saiʹ) end clear ec; ec(1:dns(1)) = 1; for k = 2:length(dns) ec(1: dns(k)) = ec(1:dns(k)) + 1; end 82
if sum(ec) ~=a meig Cond_of_A = cond(a) save jord M = I; D = I; error(ʹKich thuoc khong gian khong saiʹ) end k = 1; clear jv; while k<= dns(1) p = find(ec == ec(k)); if isempty(p) Cond_of_A = cond(a); save jord M = I; D = I; error(ʹKich thuoc khong gian khong saiʹ); end aa = I; for m = 1:ec(k) aa = aa*a1; end pp = max(size(p)); v = nulld(aa, small); jv(:,p) = v*(rand(size(v, 2), pp) ‐ 0.5)*16; k = k + pp; end clear v; for k = 1:dns(1) v(:,1) = jv(:, k); for p = 2:ec(k) v(:, p) = a1*v(:, p‐1); end vv = fliplr(v(:, 1:ec(k))); M = [M vv]; 83
end end end k = abs(det(M))^(‐1/n); M = k*M; Mi = inv(M); D = Mi*a*M; d0 = diag(D); d1 = diag(D, 1); D = diag(d0) + diag(d1, 1); function Z = nulld(A, small) norma = sqrt(mean(mean(abs(A)))); tiny = norma*eps; if nargin<2 small = (1 + norma)*sqrt(eps); end [m, n] = size(A); if m~= n error(ʹMa tran phai vuong!ʹ) end p = find(abs(A)<tiny); if ~isempty(p) A(p) = 0; end [U,S,V] = svd(A, 0); S = diag(S); s = S; norma = max(s); smax = max(s); if smax == 0; smax = 1; end s = s/smax; snorm = s; 84
t = find(s>0); if isempty(t); Z = V; return; end p = find(s 0 r = r(end); else r = 0; end Z = V(:,r+1:n); if snorm == ones(n, 1); Z = zeros(n, 0); end if max(S) <= small; Z = V; end §13. PHÂN TÍCH MA TRẬN THEO CÁC GIÁ TRỊ KÌ DỊ Phân tích ma trận theo các giá trị kì dị (kì dị value) được thực hiện trên các ma trận vuông hay chữ nhật. Ta có: [Anp] = [Unn][Snp][Vpp] Trong đó: 85
[U]T[U] = [Enn] [V]T[V] = [Epp] nghĩa là các ma trận [U] và [V] là trực giao. Các cột của [U] là các vec tơ kì dị trái, [S] có các giá trị kì dị và là ma trận đường chéo và [V]T có các hàng là các vec tơ kì dị phải. Để tính các ma trận [U], [S] và [V] ta tìm các giá trị riêng của [A][A]T và [A]T[A]. Các vec tơ riêng của [A]T[A] tạo nên các cột của [V]. Các vec tơ riêng của [A][A]T tạo nên các cột của [U]. Các giá trị kì dị của [S] là căn bậc hai của các giá trị riêng của [A][A]T hay [A]T[A]. Các giá trị riêng trên đường chéo của [S] được sắp xếp theo thứ tự giảm dần. Để hiểu được thuật toán này ta xét ví dụ sau: Ta xây dựng hàm svddecomp() để thực hiện thuật toán này: function [U, S , V] = svddecomp(A) [m, n] = size(A); if (m > n) % ta can timcac vec to kì dị phai truoc %[V D] = eigs(Aʹ*A) [V, D] = eigs(Aʹ*A); [dummy, perm] = sort(‐diag(D)); S = diag(sqrt(diag(D(perm, perm)))); V = V(:, perm); sinv = diag(1./sqrt(diag(D))); U = (A*V)*sinv; else % ta can tim cac vec to kì dị trai truoc % [U D] = eigs(A*Aʹ) [U, D] = eigs(A*Aʹ); [dummy, perm] = sort(‐diag(D)); S = diag(sqrt(diag(D(perm, perm)))); U = U(:, perm); sinv = diag(1./sqrt(diag(D))); V = sinv*(Uʹ*A); V = Vʹ; end 86
Để phân tích một ma trận ta dùng chương trình ctsvddecomp.m: clear all, clc %a = [ 1 2 3 4 5; 6 7 8 9 0; 3 4 5 6 7; 8 9 0 1 2; 2 4 6 8 1]; a = [ 1 1; 0 1; 1 0]; [u, s, v] = svddecomp(a) §14. PHÂN TÍCH SCHUR Cho ma trận vuông [A], cấp n ta phân tích nó thành: [A] = [T][U][T]* Trong đó: [T] ‐ ma trận unita và [T]* là ma trận chuyển vị liên hợp của [T](lâý chuyển vị của [T] rồi lấy liên hợp của các phần tử). [U] = [Λ] + [N] [Λ] ‐ là ma trận đường chéo có các phần tử là các giá trị riêng của [A] [N] ‐ ma trận tam giác phải, có các phần tử thuộc đường chéo chính bằng zero. Mọi ma trận vuông đều có thể phân tích Schur. Tuy nhiên phân tích này không duy nhất. Nếu [A] là ma trận trực giao thì [U] là ma trận đường chéo và các cột của [T] là các vec tơ riêng của [A]. Phân tích Schur khi này được gọi là phân tích phổ. Nếu [A] xác định dương, phân tích Schur chính là phân tích SVD. Để phân tích ma trận theo thuật toán Schur ta thực hiện các bước sau: - Tìm giá trị riêng λ1 của [A] và vec tơ riêng tương ứng [v1] - Chọn n‐ 1 vec tơ [w1], ,[wn‐1] độc lập tuyến tính và trực giao với [v1] - Tạo [V1] là ma trận có các cột là [v1], [w1], ,[wn‐1] và tính: ∗ ⎡⎤λ∗1 []VAV11[][]= ⎢⎥ ⎣⎦0A[]1 Trong đó [A1] là ma trận (n‐1)×(n‐1) - Lặp lại quá trình với ma trận [A1] ta có: ∗ ⎡⎤λ∗2 [][][]VAV212= ⎢⎥ ⎣⎦0A[]2 Trong đó [A2] là ma trận (n‐2)×(n‐2) 87
⎡⎤λ∗1 ∗ ∗ ∗ ⎢⎥ - =λ∗ Do []TAT22[][]⎢⎥ 0 1 ⎢⎥ ⎣⎦00A[]2 ⎡10⎤ TVV= ⎡⎤⎡⎤ˆ ⎡⎤ˆ = Trong đó []212⎣⎦⎣⎦ với ⎣⎦V2 ⎢ ⎥ ⎣0V[]2 ⎦ - Tiếp tục quá trình ta tìm được [V1], ,[Vn]. Cuối cùng [U] = [T]*[A][T] ⎡ ⎤ = ⎡⎤⎡⎤⎡⎤ˆˆ ⎣TVVV212n⎦ ⎣⎦⎣⎦⎣⎦L Ta xây dựng hàm schurdecom() thực hiện thuật toán trên: function [T, U] = schurdecom(a) % Phan tich Schur ma tran A n = size(a, 1); v = zeros(n, 1); v(1) = 1; b = zeros(n, n); b(:, n) = v; for k = 2:n v = a*v; b(:, n‐k+1) = v; end c = a*v; rho = ‐b\c; rho = [1 rhoʹ]; lambda = roots(rho); n = size(lambda, 1); evec = zeros(n); c = evec; e = eye(n); for i = 1:n b = a ‐ lambda(i)*e; c = nulld(b); evec(:, i) = c(:,1); 88
end p = grams(evec); T = conj(transpose(p))*a*p; U = p; Để phân tích ma trận ta dùng chương trình ctschur.m: clear all, clc a = [ 1 2 3 5; 4 5 6 2; 4 6 8 9; 9 3 6 7]; [t, u] = schurdecom(a) §15. ĐỊNH THỨC CỦA MA TRẬN Cho một ma trận vuông cấp n. Ta cần tìm định thức của nó. Trước hết chúng ta nhắc lại một số tính chất quan trọng của định thức: - nếu nhân tất cả các phần tử của một hàng (hay cột) với k thì định thức được nhân với k - định thức không đổi nếu ta cộng thêm vào một hàng tổ hợp tuyến tính của các hàng còn lại. - nếu đổi chỗ hai hàng cho nhau thì định thức đổi dấu Trước khi đi đến định nghĩa về định thức ta tìm hiểu khái niệm về hoán vị và phép thế. Cho một dãy số, nếu ta đổi chỗ các số trong dãy cho nhau thì ta đã thực hiện một phép hoán vị. Ví dụ 123, 132, là các hoán vị của dãy số {1, 2, 3}. Trong hoán vị α1α2 αi αj αn ta nói αi làm một nghịch thế với αj nếu i αj. Ví dụ trong hoán vị 1432 số 4 làm với số 3 một nghịch thế , số 4 làm với số 2 một nghịch thế, số 3 làm với số 2 một nghịch thế. Một hoán vị gọi là chẵn nếu tổng số nghịch thế trong hoán vị đó là một số chẵn; một hoán vị gọi là lẻ trong trường hợp ngược lại. Như vậy 1432 là một hoán vị lẻ. Cho một dãy số, nếu ta tạo ra một dãy số mới bằng cách đổi chỗ các phần tử cho nhau thì ta đã thực hiện một phép thế. ⎛⎞2143 Ví dụ p = ⎜⎟ là phép thế biến 2 thành 1, 1 thành 4, 4 thành 2 và 3 ⎝⎠1423 thành 3. Một phép thế gọi là chẵn nếu tính chẵn lẻ của dòng trên và dòng dưới như nhau và lẻ trong trường hợp ngược lại. Phép thế trên là phép thể lẻ. 89
Cho ma trận vuông [A] cấp n. Các phần tử của hàng thứ i là ai,1, ai,2, ,ai,n. Các phần tử của cột thứ j là a1,j, a2,j , , an,j. Ta xem hàng thứ i là một vec tơ, kí hiệu là Ai* và cột thứ j cũng là một vec tơ, kí hiệu là A*j. Với mỗi phép thế: ⎛⎞ii12L i n p = ⎜⎟ (1) ⎝⎠jj11L j n ta lập tích: aa a (2) ij11 ij 22K i n j n Trước mỗi tích (2) ta đặt dấu + nếu và dấu ‐ nếu phép thế (1) lẻ. Sau đó ta lập tổng của n! tích có dấu như vậy, nghĩa là tổng: (1)a− t(p) a a (3) ∑ ij11 ij 22K i nn j p trong đó: t(p) = 1 nếu phép thế p lẻ t(p) = 0 nếu phép thế p chẵn Tổng (4) được gọi là định thức của ma trận vuông [A], cấp n. Ta xây dựng hàm determinant() để tính định thức của ma trận theo định nghĩa: function d = determinant(A) % DETERMINANT tinh dinh thuc theo dinh nghia. [m, n] = size(A); if ( m ~= n ) fprintf ( ʹ\nʹ ); fprintf ( ʹ Chi ma tran vuong moi co dinh thuc!\nʹ ); return end p = zeros(1, n); nf = prod([1:n]); d = 0.0; for i = 1:nf p = nextperm(p); s = permsign(p); x = diag(A([1:n],p)); 90
d = d + s*prod(x); end function psign = permsign(p) % PERMSIGN tra ve dau phep the . % +1, neu phep the chan, % ‐1, neu phep the le. n = length ( p ); psign = 1; for i = 1:n‐1 j = i; while (p(j) ~= i) j = j + 1; end if ( j ~= i ) temp = p(i); p(i) = p(j); p(j) = temp; psign = ‐ psign; end end function q = nextperm(p) n = length(p); q = p; if(n == 1) q = 1; elseif (q == 0) q = [1:n]; else i = n ‐ 1; while (q(i) > q(i+1)) i = i ‐ 1; if (i == 0) break; 91
end end if (i == 0) q = [1:n]; else j = n; while (q(j) < q(i)) j = j ‐ 1; end t = q(j); q(j) = q(i); q(i) = t; q(i+1:n) = q(n:‐1:i+1); end end Để tính định thức ta dùng chương trình ctdeterminant.m: clear all, clc %a = [1 2; 3 5]; a = [1 3 5; 3 4 6; 4 6 3]; d = determinant(a) §16. TÍNH ĐỊNH THỨC BẰNG CÁCH PHÂN TÍCH MA TRẬN Cho ma trận [A]. Nếu [A] = [B]×[C] thì det[A] = det[B]×det[C] Mặt khác với một ma trận tam giác, ví dụ: ⎡⎤b11bbb 12 13 14 ⎢⎥0b b b []B = ⎢⎥22 23 24 ⎢⎥00bb33 34 ⎢⎥ ⎣⎦000b44 thì 92
det[] B= b11 b 22 b 33 b 44 nghĩa là đối với ma trận tam giác, định thức bằng tích các phần tử trên đường chéo chính. Khi phân tích ma trận [A] theo thuật toán Doolitte, ta dùng chương trình ctdoodecmp.m để tính định thức của nó: clear all, clc a = [1 2 3 4; 3 4 5 7; 2 3 8 5; 4 9 1 4]; [l, r] = doolittle(a); d = prod(diag(l))*prod(diag(r)) Khi phân tích ma trận [A] theo thuật toán Crout, ta dùng chương trình ctcrotdecmp.m để tính định thức của nó: clear all, clc a = [1 2 3 4; 3 4 5 7; 2 3 8 5; 4 9 1 4]; [l, r] = crout(a); d = prod(diag(l))*prod(diag(r)) Khi phân tích ma trận [A] theo thuật toán Choleski, ta dùng chương trình ctcholdecmp.m để tính định thức của nó: clear all, clc a = [4 ‐2 2;‐2 2 ‐4;2 ‐4 11]; a = pascal(5); l = choleski(a); d = prod(diag(l))*prod(diag(lʹ)) §17. THUẬT TOÁN TRỤ CHIÓ Cho ma trận [A] có a01,1 ≠ . Ta xây dựng ma trận [B] có các phần tử bi,j=−aa 1,1 i,j aa i,n n,j thì: 2n− []Aa= 1,1 [] B nghĩa là: 93
⎡⎤⎡⎤⎡⎤aa11 12 aa 11 13 ⎡⎤ aa 11 1n ⎢⎥det⎢⎥⎢⎥ detL det ⎢⎥ ⎢⎥⎣⎦⎣⎦aa21 22 aa 21 23 ⎣⎦ aa 21 2n ⎢⎥⎡⎤⎡⎤aa aa ⎡⎤ aa ⎢⎥11 12 11 13 11 1n n2− det⎢⎥⎢⎥ detL det ⎢⎥ det[] A= a11 det ⎢⎥⎣⎦⎣⎦aa31 32 aa 31 33 ⎣⎦ aa 31 3n ⎢⎥ ⎢⎥MMLM ⎢⎥⎡⎤⎡⎤aa11 12 aa 11 13⎡⎤ aa 11 1n ⎢⎥det⎢⎥⎢⎥ detL det ⎢⎥ ⎣⎦⎣⎦⎣⎦aan1 n2 aa n1 n3⎣⎦ a n1nna Ta xây dựng hàm chiopivot() để thực hiện thuật toán trên: function d = chiopivot(a) %tinh dinh thuc bang thuat toan Chio pivotal condensation if a(1, 1) == 0 error(ʹKhong dung phuong phap nay de tinh dinh thuc duoc !ʹ); return; end c = a(1, 1); n = size(a, 1); if (n = 1) for i = 1:m%hang b(i, 1:m) = a(1, 1)*a(i+1, 2:m+1) ‐ a(i+1, 1)*a(1, 2:m+1); end if (m > 2) a = b; c = c*a(1,1); clear b; end m = m ‐ 1; end d = b(1, 1)*b(2, 2) ‐ b(2, 1)*b(1, 2); 94
d = d/c; Để tính định thức ta dùng chương trình ctchiopivot.m: clear all, clc a = [1 2 3 4; 3 4 5 7; 2 3 8 5; 4 9 1 4]; d = chiopivot(a) §18. THUẬT TOÁN LAPLACE Để tính định thức theo thuật toán Laplace, ta khai triển định thức theo hàng hay cột. Cho ma trận vuông [A] cấp n. Nếu bỏ đi hàng i và cột j (tức xoá hàng và cột chứa phần tử aij) thì ta có một ma trận cấp (n ‐ 1), định thức của nó gọi là định thức con cấp (n ‐ 1) ứng với phần tử aij (minor) , ký hiệu là Mij. Ta chú ý đến hàng thứ i và aij là một phần tử của hàng đó. Trong det[A] ta gộp những số hạng chứa aij lại và đặt aij làm thừa số chung, hệ số của nó kí hiệu là Aij và gọi là phần bù đại số (cofactor) của phần tử aij. Cofactor Aij của ma trận [A] là: ij+ A(1)Mij=− ij Định thức của [A] khi khai triển theo hàng là: n det[] A= ∑ aij A ij i1= Ta xây dựng hàm cofactor() để tính các phần bù đại số: function c = cofactor(A, i, j) % cac phan bu dai so cua ma tran % dung de nghich dao mt % C = cofactor(A, i, j) tra ve phan bu dai so cua %ng i, cot j cua A. if nargin == 3 M = A; M(i,:) = []; M(:,j) = []; c = (‐1)^(i+j) * det(M); else 95
[n, n] = size(A); for i = 1:n for j = 1:n c(i,j) = cofactor(A, i, j); end end end Sau khi phần bù đại số, ta xây dựng hàm cofactordet() để tính định thức của [A] function d = cofactordet(A) d = 0; for i = 1:size(A, 1) c = cofactor(A, i, 1); d = d + A(i, 1)*c; end Để tính định thức ta dùng chương trình ctcofactordet.m: clear all, clc a = [1 2 3 4; 3 4 5 7; 2 3 8 5; 4 9 1 4]; det = cofactordet(a) §19. THUẬT TOÁN DODGSON Thuật toán cô đặc Dodgson ( Dodgson condensation) dùng để tính định thức của ma trận vuông. Nó được Charles Ludwidge Dodgson đưa ra. Để tính định thức của ma trận cấp n × n, ta xây dựng các ma trận cấp (n ‐ 1) × (n ‐ 1) cho đến ma trận cấp 1 × 1 là định thức của ma trận cần tìm. Bước đầu tiên ta xây dựng ma trận cấp (n ‐ 1)×(n ‐ 1) từ các định thức của các ma trận con 2×2. Ví dụ với ma trận ⎡⎤510 ⎢⎥ ⎢⎥285 ⎣⎦⎢⎥067 96
ta có các ma trận con là: ⎡⎤⎡⎤⎡⎤⎡⎤51 10 28 85 ⎢⎥⎢⎥⎢⎥⎢⎥ ⎣⎦⎣⎦⎣⎦⎣⎦28 85 06 67 Các ma trận này sẽ tạo ra các phân tử của ma trận 2×2. Phần tử ử hàng r, cột c là định thức của ma trận con 2×2 của ma trận ban đầu với hàng r và cột c ở góc trên trái. Như vậy ma trận mới là: ⎡⎤38 5 ⎢⎥ ⎣⎦12 26 Ma trận k×k được tạo ra bằng cách lấy định thức của ma trận con 2×2 của ma trận (k+1)×(k+1) như ở trên và chia nó cho phần tử tương ứng của ma trận trung tâm, nghĩa là ma trận đã bỏ đi hàng trên cùng, hàng dưới cùng, cột bên phải và cột bên trái của ma trận (k+2)×(k+2) Ta xây dựng hàm dodgson() để thực hiện thuật toán trên: function dt = dodgson(a) if size(a, 1) ~= size(a, 2) error(ʹMa tran A phai la mt tran vuongʹ); end; n = size(a, 1); if n == 2 dt = a(1, 1)*a(2, 2) ‐ a(2, 1)*a(1, 2); return end; if n == 3; for i = 1:n‐1 b(i, 1:n‐1) = a(i, 1:n‐1).*a(i+1, 2:n) ‐ a(i+1, 1:n‐1).*a(i, 2:n); end dt = (b(1, 1)*b(2, 2) ‐ b(2, 1)*b(1, 2))/a(2,2); return end c = a; c(:, 1) = []; c(:, n‐1) = []; c(1, :) = []; 97
c(n ‐ 1, :) = []; for i = 1:n‐1 b(i, 1:n‐1) = a(i, 1:n‐1).*a(i+1, 2:n) ‐ a(i+1, 1:n‐1).*a(i, 2:n); end m = size(b, 1); while m >= 2 for i = 1:m‐1 for j = 1:m‐1 d(i, j) = (b(i, j)*b(i+1, j+1) ‐ b(i+1, j)*b(i, j+1))/c(i, j); end end if m > 3 c = b; c(:, 1) = []; c(:, m‐1) = []; c(1, :) = []; c(m ‐ 1, :) = []; b = d; end m = m ‐ 1; end dt = (d(1, 1)*d(2, 2) ‐ d(2, 1)*d(1, 2))/b(2,2); Để tính định thức ta dùng chương trình ctdodgson.m: clear all, clc; a = [1 2 3 4; 5 1 3 2; 4 9 2 2; 6 3 4 1]; dt = dodgson(a) §20. NGHỊCH ĐẢO MA TRẬN BẰNG CÁCH DÙNG MINOR Cho ma trận [A], ta có: A ()a−1 = j,i i,j det[] A Trong đó: 98