Bài giảng Matlab - Sử dụng symbolic math toolbox

0.9 1 x 1/(5+4 cos(x)) -6 -4 -2 0 2 4 6 -5 -4 -3 -2 -1 0 1 2 x 32/(5+4 cos(x))3 sin(x)2+4/(5+4 cos(x))2 cos(x) 39 Từ đồ thị ta thấy rằng giá trị của f”(x) nằm trong khoảng [-4 , 1].Giá trị max vàmin của f”(x) xuất hiện tại f”’(x)=0.Phát biểu: f3 = diff(f2); cho 32/(5+4*cos(x))^3*sin(x)^2+4/(5+4*cos(x))^2*cos(x) và : pretty(f3) cho: 234 3 ))xcos(45( )xsin( 4 ))xcos(45( )xcos()xsin( 96 ))xcos(45( )xsin( 384 +−+++ Ta rút gọn f3 và viết lại d−ới dạng dễ đọc: f3 = simple(f3); pretty(f3) Kết quả là: 4 22 ))xcos(45( )25)xcos(80)xcos(80)xsin(96)(xsin( 4 + −++ Bây giờ ta tìm các giá trị zero cuả f3 bằng lệnh: z = solve(f3) kết quả cho ta ma trận: z = [ 0] [ atan((-255-60*19^(1/2))^(1/2) , 10+3*19^(1/2))] [ atan(-(-255-60*19^(1/2))^(1/2), 10+3*19^(1/2))] [ atan((-255+60*19^(1/2))^(1/2)/(10-3*19^(1/2)))+pi] [ -atan((-255+60*19^(1/2))^(1/2)/(10-3*19^(1/2)))-pi] Mỗi hàng là một nghiệm của f”’(x).Lệnh: format; % Default format of 5 digits zr = double(z) converts the zeros to double form. zr = 0 0 0 2.4483 -2.4483 Nh− vậy ta đã tìm đ−ợc 5 nghiệm.Tuy nhiên đồ thị của f3 cho thấy ta ch−a tìm đủ nghiệm của nó. -6 -4 -2 0 2 4 6 -3 -2 -1 0 1 2 3 x Zero cua f3 40 ezplot(f3) hold on; plot(zr,0*zr,'ro') plot([-2*pi,2*pi], [0,0],'g-.'); title('Zeros of f3') Điều này xảy ra do f”’(x) chứa số hạng sinx,bằng 0 tại các giá trị nguyên lần π nh−ng hàm solve(sin(x)) lại chỉ đ−a ra giá trị 0 tại x = 0.Chúng ta có thể nhận đ−ợc tất cả các nghiệm bằng cách biến đổi zr = [0 zr(4) pi 2*pi -zr(4)] bằng cách nhân 2π và có zr = [zr-2*pi zr zr+2*pi] Bây giờ ta vẽ zr đã biến đổi lên đồ thị của f3: plot(zr,0*zr,'kx') ylabel('f2'); title('Ve do thi f2 = f''''(x)') hold on plot(0,double(f20),'ro') text(-1,-0.25,'Local minimum') Kết quả đồ thị nh− sau: Từ đồ thị ta thấy rằng điểm cực tiểu xảy ra tại x gần ±π.Ta có thể tính chính xác là diểm cực tiểu đúng tại ±π bằng cách dùng các lệnh theo trình tự sau.Tr−ớc hết ta thay ±π vào f”’(x): simple([subs(f3,x,-pi),subs(f3,x,pi)]) -6 -4 -2 0 2 4 6 -3 -2 -1 0 1 2 3 x Zero cua f3 Điểm 0 đầu tiên của f”’(x) tìm bởi solve là tại x = 0.Chúng ta thay thế 0 vào biến chữ trong f2: f20 = subs(f2,x,0) để tìm giá trị t−ơng ứng của f”(0).Kết quả là: f20 = 0.0494 Trên đồ thị của f”(x) giá trị này chỉ là cực tiểu địa ph−ơng.Ta thể hiển điều này trên đồ thị bằng các lệnh: clf ezplot(f2) axis([-2*pi 2*pi -4.25 1.25]) -6 -4 -2 0 2 4 6 -4 -3 -2 -1 0 1 x Ve do thi f2 = f''(x) f2 Local minimum 0 1 Ve do thi f2 = f''(x) Local minimum 41 Kết quả: ans = [ 0, 0] Nh− vậy x = ±π là điểm đặc biệt của f”’(x).Ta thấy rằng x = ±π là điểm cực tiểu toàn cục của f2. m1 = double(subs(f2,x,-pi)); m2 = double(subs(f2,x,pi)); plot(-pi,m1,'go',pi,m2,'go') text(-1,-4,'Global minima') Giá trị cực tiểu đó là: [m1 m2] ans = -4 -4 Các phân tích trên cho thấy là phạm vi giá trị của f”(x) là từ [ -4 ,1].Ta tiếp tục kiểm tra các điểm 0 khác cho bởi solve.Tr−ớc hết ta tách nghiệm thứ 4 trong z và gán nó cho một biến riêng: s = z(4) và nhận đ−ợc kết quả: s = atan((-255+60*19^(1/2))^(1/2)/(10-3*19^(1/2)))+pi Thực hiện: sd = double(s) để nhận đ−ợc giá trị số của s sd = 2.4483 Ta vẽ điểm (s,f2(s) theo f2: M1 = double(subs(f2,x,s)); plot(sd,M1,'ko') text(-1,1,'Global maximum') để thấy đ−ợc là s là điểm max.Giá trị max này là M1 = 1.0051 Bây giờ ta tích phân f”(x) hai lần bằng lệnh: g = int(int(f2)) và có kết quả: g = -8/(tan(1/2*x)^2+9) 42 Đây không phải là hàm f(x) ta xét ban đầu.Sai khác giữa g(x) và f(x) là: d = f - g cho ta: d = 1/(5+4*cos(x))+8/(9+tan(1/2*x)^2) pretty(d) 2)x2/1tan(9 8 )xcos(45 1 +++ Ta có thể rút gọn d bằng lệnh simple(d) hay simplify(d).Cả hai cho kết quả: ans = 1 Điều này minh hoạ cho khái niệm là đạo hàm hàm f(x) hai lần và rồi tích phân kết quả hai lần ta nhận đ−ợc một hàm khác với f(x) bởi một hàm tuyến tính của x. Cuối cùng tích phân f(x) một lần ta có: F = int(f) F = 2/3*atan(1/3*tan(1/2*x)) bao gồm cả hàm arctan.Nh− vậy F(x) là nguyên hàm của một hàm liên tục nh−ng bản thân lại là hàm không liên tục mà có đồ thị nh− sau: ezplot(F) Hàm F(x) gián đoạn tại ±π. Đ4. Rút gọn và thay số 1. Rút gọn biểu thức:Ta xét 3 biểu thức khác nhau: syms x f = x^3-6*x^2+11*x-6 g = (x-1)*(x-2)*(x-3) -6 -4 -2 0 2 4 6 -1 -0.5 0 0.5 1 x 2/3 atan(1/3 tan(1/2 x)) 43 h = x*(x*(x-6)+11)-6 Thực hiện các lệnh: pretty(f), pretty(g), pretty(h) ta nhận đ−ợc: f = x3 - 6x2 +11x-6 g = (x-1)(x-2)(x-3) h = x(x(x-6)+11)-6 Cả 3 biểu thức này là các dạng biểu diễn toán học khác nhau của cùng một hàm toán học-đó là đa thức bậc 3 theo x.Mỗi một dạng thích hợp với một dạng tính toán.Dạng thứ nhất f là dạng chung nhất th−ờng đ−ợc dùng biểu diễn đa thức.Nó đơn giản là một tổ hợp tuyến tính của các số mũ của x.Dạng thứ 2,hàm g,là dạng phân tích thành thừa số.Nó biểu diễn nghiệm của đa thức.Tuy nhiên không phai đa thức nào cũng có nghiệm,nghĩa là có thể phân tích thành thừa số.Dạng thứ 2 là dạng Horner của đa thức.Nó rất tiện dùng để tính trị số của đa thức tại một giá trị nào đó của x. Symbolic Math Toolbox cung cấp một số hàm dùng để biến đổi các biểu thức đại số và l−ợng giác thành các biểu thức đơn giản hơn.Chúng gồm: collect,expand, horner, factor, simplify, và simple. a.collect:Phát biểu: collect(f) xem f nh− một đa thức gồm các biến chữ x và gộp tất cả các hệ cùng bậc của x.Đối số thứ 2 của chỉ rõ biến định gộp nếu có nhiều iến trong biểu th−c.Sau đây là một số ví dụ: f collect(f) (x-1)(x-2)(x-3) x^3-6*x^2+11*x-6 x*(x*(x-6)+11)-6 x^3-6*x^2+11*x-6 (1+x)*t + x*t 2*x*t+t b.expand:Phát biểu: expand(f) khai triển biểu thức.Sau đây là một số ví dụ: f expand(f) a*(x+y) a*x+a*y (x-1)*(x-2)*(x-3) x^3-6*x^2+11*x-6 x*(x*(x-6)+11)-6 x^3-6*x^2+11*x-6 exp(a+b) exp(a) + exp(b) cos(x+y) cos(x)*cos(y)-sin(x)*sin(y) cos(3*acos(x)) 4*x^3-3*x c.horner:Phát biểu: horner(f) biến đổi một đa thức thành dạng Horner hay biểu diễn lồng nhau.Ví dụ: f horner(f) x^3-6*x^2+11*x-6 -6+(11+(-6+x)*x)*x 1.1+2.2*x+3.3*x^2 11/10+(11/5+33/10*x)*x 44 d.factor:Nếu f là đa thức hệ số hữu tỉ,phát biểu: factor(f) biểu diễn f nh− là tích của các đa thức có bậc thấp hơn với hệ số hữu tỷ.Ví dụ: f factor(f) x^3-6*x^2+11*x-6 (x-1)*(x-2)*(x-3) x^3–6*x^2+11*x–5 x^3–6*x^2+11*x–5 x^6+1 (x^2+1)*(x^4–x^2+1) Đây là một ví dụ khác về phân tích đa thức xn +1 thành thừa số: syms x; n = 1:9; x = x(ones(size(n))); p = x.^n + 1; f = factor(p); [p; f].' trả về ma trận với các đa thức ở cột thứ nhất và các thừa số ở côth thứ 2: [ x+1, x+1 ] [ x^2+1, x^2+1 ] [ x^3+1, (x+1)*(x^2-x+1) ] [ x^4+1, x^4+1 ] [ x^5+1, (x+1)*(x^4-x^3+x^2-x+1)] [ x^6+1, (x^2+1)*(x^4-x^2+1) ] [ x^7+1, (x+1)*(1-x+x^2-x^3+x^4-x^5+x^6) ] [ x^8+1, x^8+1 ] [ x^9+1, (x+1)*(x^2-x+1)*(x^6-x^3+1) ] Hàm factor có thể phân tích các đối t−ợng chữ có chứa số nguyên thành thừa số.Ví dụ: one = '1' for n = 1:11 N(n,:) = sym(one(1,ones(1,n))); end [N factor(N)] cho kết quả: [ 1, 1 ] [ 11, (11) ] [ 111, (3)*(37) ] [ 1111, (11)*(101) ] [ 11111, (41)*(271) ] [ 111111, (3)*(7)*(11)*(13)*(37) ] [ 1111111, (239)*(4649) ] [ 11111111, (11)*(73)*(101)*(137) ] [ 111111111, (3)^2*(37)*(333667) ] [ 1111111111, (11)*(41)*(271)*(9091)] [ 11111111111, (513239)*(21649) ] 45 e.simplify:Hàm simplify là một hàm mạnh,dùng rút gọn các biểu thức.Sau đây là một số ví dụ: f simplify(f) x*(x*(x-6)+11)-6 x^3-6*x^2+11*x-6 (1-x^2)/(1-x) x + 1 (1/a^3+6/a^2+12/a+8)^(1/3) ((2*a+1)^3/a^3)^(1/3) syms x y positive log(x*y) log(x) + log(y) exp(x) * exp(y) exp(x+y) besselj(2,x) - ... 2*besselj(1,x)/x + besselj(0,x) 0 gamma(x+1)-x*gamma(x) 0 cos(x)^2 + sin(x)^2 1 f.simple: Hàm simple đ−a ra dạng ngắn nhất có thể có của một biểu thức.Hàm này có nhiều dạng,mỗi dạng trả về kết quả khác nhau.Dạng: simple(f) hiển thị dạng ngắn nhất.Ví dụ: syms x simple(cos(x)^2 + sin(x)^2) ans = 1 Trong một số tr−ờng hợp,áp dụng simple 2 lần để nhận đ−ợc hiệu quả rút gọn cao hơn.Ví dụ: syms a f = (1/a^3+6/a^2+12/a+8)^(1/3); simple(simple(f)) cho ta: 1/a+2 Trong khi lệnh : syms a simple(f) cho ta: (2*a+1)/a Hàm simple đặc biệt có hiệu quả trên các biểu thức l−ợng giác.Sau đây là một số ví dụ: f simple(f) cos(x)^2+sin(x)^2 1 2*cos(x)^2-sin(x)^2 3*cos(x)^2-1 cos(x)^2-sin(x)^2 cos(2*x) cos(x)+(-sin(x)^2)^(1/2) cos(x)+i*sin(x) cos(x)+i*sin(x) exp(i*x) cos(3*acos(x)) 4*x^3-3*x 46 2. Thay số: Có hai hàm dùng để thay trị là subexpr và subs a. subexpr:Lệnh : syms a x s = solve(x^3+a*x+1) giải ph−ơng trình : x^3+a*x+1 = 0 theo x. Kết quả: [ 1/6*(- 108+12*(12*a^3+81)^(1/2))^(1/3)-2*a/(-108+12*(12*a^3+81)^(1/2))^(1/3)] [ -1/12*(-108+12*(12*a^3+81)^(1/2))^(1/3)+a/(- 108+12*(12*a^3+81)^(1/2))^(1/3)+1/2*i*3^(1/2)*(1/6*(- 108+12*(12*a^3+81)^(1/2))^(1/3)+2*a/(-108+12*(12*a^3+81)^(1/2))^(1/3))] [ -1/12*(-108+12*(12*a^3+81)^(1/2))^(1/3)+a/(-108+12*(12*a^3+81)^(1/2))^(1/3)- 1/2*i*3^(1/2)*(1/6*(-108+12*(12*a^3+81)^(1/2))^(1/3)+2*a/(- 108+12*(12*a^3+81)^(1/2))^(1/3))] Dùng lệnh pretty để nhận đ−ợc dạng dễ đọc hơn: [ 1/3 a ] [ 1/6 %1 - 2 ----- ] [ 1/3 ] [ %1 ] [ ] [ 1/3 a 1/2 / 1/3 a \ ] [- 1/12 %1 + ----- + 1/2 i 3 | 1/6 %1 + 2 ----- | ] [ 1/3 | 1/3 | ] [ %1 \ %1 / ] [ ] [ 1/3 a 1/2 / 1/3 a \ ] [- 1/12 %1 + ----- - 1/2 i 3 | 1/6 %1 + 2 ----- | ] [ 1/3 | 1/3 | ] [ %1 \ %1 / ] 3 1/2 %1 := -108 + 12 (12 a + 81) Lệnh pretty thừa kế khái niệm %n(n là một số nguyên) từ Maple để định nghĩ biểu thức con gặp hiều lần trong đối t−ợng chữ.Hàm subexpr cho phép ta l−u các biểu thức con này cũng nh− các đối t−ợng chữ đ−ợc viết trong biểu thức con.Các biểu thức con đ−ợc l−u trong một ma trận cột gọi là sigma. Tiếp tục ví dụ của ta: r = subexpr(s) cho ta sigma = -108+12*(12*a^3+81)^(1/2) r = [ 1/6*sigma^(1/3)-2*a/sigma^(1/3)] [ -1/12*sigma^(1/3)+a/sigma^(1/3)+1/2*i*3^(1/2)*(1/6*sigma^(1/3)+2*a/sigma^(1/3))] 47 [ -1/12*sigma^(1/3)+a/sigma^(1/3)-1/2*i*3^(1/2)*(1/6*sigma^(1/3)+2*a/sigma^(1/3))] ta thấy rằng subexpr tạo biến signma rg vùng làm việc của MATLAB. b.subs:Ta tìm giá trị riêng và vec tơ riêng của ma trận vòng A: syms a b c A = [a b c; b c a; c a b]; [v,E] = eig(A) v = [ 1, -(a+(b^2-b*a-c*b-c*a+a^2+c^2)^(1/2)-b)/(a-c), -(a-(b^2-b*a- c*b-c*a+a^2+c^2)^(1/2)-b)/(a-c)] [ 1, -(b-c-(b^2-b*a-c*b-c*a+a^2+c^2)^(1/2))/(a-c), -(b-c+(b^2-b*a- c*b-c*a+a^2+c^2)^(1/2))/(a-c)] [ 1, 1, 1] E = [ b+a+c, 0, 0] [ 0, (b^2-b*a-c*b-c*a+a^2+c^2)^(1/2), 0] [ 0, 0, -(b^2-b*a-c*b-c*a+a^2+c^2)^(1/2)] Giả sử ta muốn thay biểu thức khá dài: (b^2-b*a-c*b-c*a+a^2+c^2)^(1/2) trong v và E.Tr−ớc hết ta dùng subexpr: v = subexpr(v,'S') cho ta kết quả: S = (b^2-b*a-c*b-c*a+a^2+c^2)^(1/2) v = [ -(a+S-b)/(a-c), -(a-S-b)/(a-c), 1] [ -(b-c-S)/(a-c), -(b-c+S)/(a-c), 1] [ 1, 1, 1] Sau đó thay S vào E: E = subs(E,S,'S') E = [ S, 0, 0] [ 0, -S, 0] [ 0, 0, b+c+a] Bây giờ giả sử ta muốn tính v khi a = 10.Ta dùng lệnh sau: subs(v,a,10) sẽ thay các biến a trong v bằng số 10: [ -(10+S-b)/(10-c), -(10-S-b)/(10-c), 1] [ -(b-c-S)/(10-c), -(b-c+S)/(10-c), 1] [ 1, 1, 1] Chú ý là các biểu thức có S không bị ảnh h−ởng gì cả,nghĩa là biến a trong S không đ−ợc thay bằng 10.Hàm subs là hàm hữu ích để thay thế nhiều giá trị của nhiều biến trong một biểu thức.Ta xem S.Giả sử ngoài việc thay a =10 ta cũng muốn thay giá trị b = 2 và c = 10 vào biểu 48 thức.Cách đơn giản nhất là đặt giá trị a,b,c trong vùng làm việc của MATLAB.Sau đó subs sẽ tính kết quả. a = 10; b = 2; c = 10; subs(S) ans = 8 Lệnh subs có thể kết hợp với lệnh double để tính trị số của một iểu thức chữ.Giá sử ta có: syms t M = (1-t^2)*exp(-1/2*t^2); P = (1-t^2)*sech(t); và muốn xem trên đồ thị P và M khác nhau nh− thế nào. Ta dùng các lệnh: ezplot(M); hold on; ezplot(P) và đồ thị nh− hình vẽ.Tuy nhiên ta vẫn khó hình dung đ−ợc sự sai khác giữa hai đđ−ờng cong. Vì vậy tốt hơn chúng ta kết hợp subs,double lại T =-6:0.05:6; MT = double(subs(M,t,T)); PT = double(subs(P,t,T)); plot(T,MT,'b',T,PT,'r-.') title(' ') legend('M','P') xlabel('t'); grid để tạo ra đồ thị nhiều màu. -6 -4 -2 0 2 4 6 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 t (1-t2) sech(t) -6 -4 -2 0 2 4 6 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 t M P 49 Đ5. Dại số tuyến tính 1. Các toán tử đại số cơ bản:Các toán tử đại số cơ bản trong các đối t−ợng chữ cũng là các toán tử trong các đối t−ợng số của MATLAB.điều này đ−ợc minh hoạ trong ví dụ sau.Biến đổi Givens tạo ra một mặt phẳng quay một góc t.Các phát biểu: syms t; G = [cos(t) sin(t); -sin(t) cos(t)] tạo ra ma trận biến đổi: G = [ cos(t), sin(t) ] [ -sin(t), cos(t) ] Dùng biến đổi Givens 2 lần đơn giản là quay một góc gấp đôi.Ma trận t−ơng ứng có thể tính bằng cách nhân G với chính nó hay nâng G lên luỹ thừa 2.Cả hai cáh tạo ra cùng một kết quả: A = G*G hay A = G^2 cho ta: A = [cos(t)^2-sin(t)^2, 2*cos(t)*sin(t)] [ -2*cos(t)*sin(t), cos(t)^2-sin(t)^2] Hàm simple A = simple(A) tạo ra kết quả: A = [ cos(2*t), sin(2*t)] [-sin(2*t), cos(2*t)] Một phép quay Givens là ma trận trực giao nên chuyển vị của nó chính là nghịc đảo của nó. Thật vậy: I = G.' *G tạo ra ma trận I: I = [cos(t)^2+sin(t)^2, 0] [ 0, cos(t)^2+sin(t)^2] và: I = simple(I) I = [1, 0] [0, 1] 2. Các toán tử đại số tuyến tính:Ta xét các toán tử đại số tuyến tính cơ bản.Lệnh: H = hilb(3) tạo ra ma trận Hilbert 3 ì 3.Với format short,MATLAB in ra: H = 1.0000 0.5000 0.3333 0.5000 0.3333 0.2500 0.3333 0.2500 0.2000 50 Các phần tử đ−ợc tính toán của H là số dấu chấm động mà nó là tỉ số của các số nguyên nhỏ nhất.Thật vậy.H là ma trận MATLAB lớp double.Biến đổi H thành ma trận chữ: H = sym(H) cho kết quả: H= [ 1, 1/2, 1/3] [1/2, 1/3, 1/4] [1/3, 1/4, 1/5] Điều này cho phép các toán tử kí hiệu sym(hilb(3)) trên H tạo ra kết quả t−ơng ứng với ma trận Hilbert có độ chính xác vô cùng,chứ không phải với độ chính xác của số thực dấu chấm động.Do đó: inv(H) tạo ra ma trận: [ 9, -36, 30] [-36, 192, -180] [ 30, -180, 180] và det(H) cho kết quả: 1/2160 Chúng ta có thể dùng dấu \ để giaie hệ ph−ơng trình tuyến tính.Lệnh: b = [1 1 1]' x = H\b % Solve Hx = b tạo ra nghiệm: x = [ 3] [ -24] [ 30] Tất cả kết quả này,nghịch đảo,định thức và nghiệm của hệ ph−ơng trình tuyến tính là kết quả chính xác t−ơng ứng với độ chính xác vô cùng,hữu tỉ của ma trận Hilbert.Mặt khác dùng digits(16) lệnh: V = vpa(hilb(3))%vpa là kết quả tính chính xác theo Maple cho kết quả: [ 1. , .5000000000000000, .3333333333333333] [.5000000000000000, .3333333333333333, .2500000000000000] [.3333333333333333, .2500000000000000, .2000000000000000] Kết quả đ−ợc làm tròn tới 16 chữ số.Khi nghich đảo ma trận sai số tăng lên.Do vậy: inv(V) cho ta: [ 9.000000000000082, -36.00000000000039, 30.00000000000035] [-36.00000000000039, 192.0000000000021, -180.0000000000019 ] [ 30.00000000000035, -180.0000000000019, 180.0000000000019] và: det(V) cho: .462962962962958e-3 51 và V\b which is [ 3.000000000000041] [-24.00000000000021] [ 30.00000000000019] Vì ma trận H không suy biến nên không gian null của H null(H) ans = [ empty sym ] và không gian cột của H: colspace(H) ans = [ 0, 0, 1] [ 1, 0, 0] [ 0, 1, 0] tạo ra ma trận rỗng và hoán vị của ma trận đơn vị. Ta tìm giá trị s của H(1,1) để làm cho H suy biến.Lệnh: syms s H(1,1) = s Z = det(H) sol = solve(Z) tạo ra các kết quả: H = [ s, 1/2, 1/3] [1/2, 1/3, 1/4] [1/3, 1/4, 1/5] Z = 1/240*s-1/270 sol = 8/9 Vậy: H = subs(H,s,sol) thay thế giá trị sol đã tính vào H: H = [8/9, 1/2, 1/3] [1/2, 1/3, 1/4] [1/3, 1/4, 1/5] Bây giờ lệnh: det(H) cho: ans = 0 và: inv(H) tạo ra thông báo lỗi: 52 ??? error using ==> inv Error, (in inverse) singular matrix vì H là ma trận suy biến.Với ma trận này , Z = null(H) và C = colspace(H) là không tầm th−ờng: Z = [ 1] [-4] [10/3] C = [ 0, 1] [ 1, 0] [6/5, -3/10] 3. Các giá trị riêng:Các giá trị riêng và véc tơ riêng bằng chữ của ma trận vuông A đ−ợc tính bằng lệnh: E = eig(A) [V,E] = eig(A) Nếu dùng độ chính xác có thể biến đổi: E = eig(vpa(A)) [V,E] = eig(vpa(A)) Giá trị riêng của ma trận A là nghiệm của đa thức đặc tính của A,det(Ax*I),đ−ợc tính bằng: poly(A) Ta xét ma trận H ở ví dụ tr−ớc: H = [8/9, 1/2, 1/3] [1/2, 1/3, 1/4] [1/3, 1/4, 1/5] Ma trận này suy biến nên một trong các giá trị riêng phải bằng 0.Phát biểu: [T,E] = eig(H) tạo ra ma trận T và E.Các cột của T là các vec tơ riêng của H. T = [ 1, 28/153+2/153*12589^(1/2), 28/153-2/153*12589^(12)] [ -4, 1 , 1] [ 10/3, 92/255-1/255*12589^(1/2) , 292/255+1/255*12589^(12)] T−ơng tự,các phần tử trên đ−ờng chéo của E là giá trị riệng của H E = [0, 0, 0] [0, 45+1/180*12589^(1/2), 0] [0, 0, 32/45-1/180*12589^(1/2)] Để dễ thấy cấu trúc của ma trận T và E ta tính cụ thể trị của nó: Td = double(T) Ed = double(E) Kết quả là: Td = 1.0000 1.6497 -1.2837 -4.0000 1.0000 1.0000 53 3.3333 0.7051 1.5851 Ed = 0 0 0 0 1.3344 0 0 0 0.0878 Giá trị riêng đầu tiên là 0.Vec tơ riêng t−ơng ứng(cột đầu tiên của Td) là cơ sở của không gian null tìm đ−ợc trong phần tr−ớc.Hai giá trị riêng khác là kết quả của việc áp dụng công thức cầu ph−ơng đối với biểu thức x^2-64/45*x+253/2160 mà kết quả chứa trong factor(poly(H)). syms x g = simple(factor(poly(H))/x); solve(g) Dạng gần với biểu thức của giá trị riêng chí có thể có khi đa thức đặc tính có thể biểu diễn bằng tích của dâ thức hữu tỉ bậc 4 hay nhỏ hơn.Ma trận Rosser là ma trận thử minh hoạ cho yêu cầu này.Phát biểu: R = sym(gallery('rosser')) tạo ra: R = [ 611, 196, -192, 407, -8, -52, -49, 29] [ 196, 899, 113, -192, -71, -43, -8, -44] [ -192, 113, 899, 196, 61, 49, 8, 52] [ 407, -192, 196, 611, 8, 44, 59, -23] [ -8, -71, 61, 8, 411, -599, 208, 208] [ -52, -43, 49, 44, -599, 411, 208, 208] [ -49, -8, 8, 59, 208, 208, 99, -911] [ 29, -44, 52, -23, 208, 208, -911, 99] Lệnh: p = poly(R); pretty(factor(p)) tạo ra: [x (x - 1020) (x2 - 1020 x + 100)(x2- 1040500) (x - 1000)2 ] Đa thức đặc tính (bậc 8) đ−ợc phân tích thành tích của 2 số hạng tuyến tính và 3 số hạng bình ph−ơng.Chúng ta có thể thấy ngay rằng 4 giá trị riêng là 0,1020,và giá trị riêng kép 1000.Bốn giá trị riêng khác nhận đ−ợc từ các biểu thức bậc 2 còn lại.Dùng lệnh: eig(R) ta tìm đ−ợc tất cả các các giá trị riêng: [ 0] [ 1020] [ 510+100*26^(1/2)] [ 510-100*26^(1/2) ] [ 10*10405^(1/2)] [ -10*10405^(1/2)] [ 1000] [ 1000] 54 Ma trận Rosser không phải là ví dụ điển hình.Rất hiếm ma trận bậc 8 có thể phân tích thành các thừa số đơn giản.Ta thay đổi phần tử ở góc từ 29 thnàh 30 bằng lệnh: S = R; S(1,8) = 30; S(8,1) = 30; và rồi dùng tiếp lệnh: p = poly(S) ta có: p = 40250968213600000+51264008540948000*x- 1082699388411166000*x^2+4287832912719760*x^-3- 5327831918568*x^4+82706090*x^5+5079941*x^6- 4040*x^7+x^8 và nếu factor(p) ta đ−ợc chính nó,nghĩa là đa thức đặc tính không thể phân tích thành các thừa số là đa thức hữu tỉ.Với ma trận Rosser đã biến đổi: F = eig(S) cho ta: F = [ -1020.0532142558915165931894252600] [ -.17053529728768998575200874607757] [ .21803980548301606860857564424981] [ 999.94691786044276755320289228602] [ 1000.1206982933841335712817075454] [ 1019.5243552632016358324933278291] [ 1019.9935501291629257348091808173] [ 1020.4201882015047278185457498840] Chú ý là các giá trị này gần với các giá trị của ma trận Rosser ban đầu.Hơn nữa,giá trị số của F là kết quả của phép tính dấu chấm động của Maple. Ta cũng có thể tính các giá trị riêng của ma trận chữ,nh−ng nghiệm dạng gần nhau là rất hiếm.Biến đổi Givens đ−ợc tạo ra nh− là luỹ thừa của ma trận ban đầu. Các lệnh Symbolic Math Toolbox: syms t A = sym([0 1; 1 0]); G = expm(t*A) cho: [ cos(t), sin(t)] [ -sin(t), cos(t)] Tiếp theo lệnh: g = eig(G) tạo ra: g = [ cos(t)+(cos(t)^2-1)^(1/2)] [ cos(t)-(cos(t)^2-1)^(1/2)] Ta có thể dùng lệnh simple để rút gọn dạng này của g.Thay vì lặp lại simple: for j = 1:4 [g,how] = simple(g) 55 end tạo ra kết quả tốt nhất: g = [ cos(t)+(-sin(t)^2)^(1/2)] [ cos(t)-(-sin(t)^2)^(1/2)] how = simplify g = [ cos(t)+i*sin(t)] [ cos(t)-i*sin(t)] how = radsimp g = [ exp(i*t)] [ 1/exp(i*t)] how = convert(exp) g = [ exp(i*t)] [ exp(-i*t)] how = combine Ta thấy rằng ứng dụng đầu tiên của simple dùng simplify để tạo ta tổng của các số hạng sin và cosin.Tiếp đó simple gọi radsim để tạo ra cos(t) + i*sin(t) đối với vec tơ riêng đầu tiên.Lần thứ 3,simple dùng convert(exp) để thay đổi sin và cosin thành số mũ phức.Lần cuối simple dùng combine để có kết quả cuối cùng. 4. Dạng Jordan chính tắc:Dạng Jordan chính tắc xuất phát từ ý t−ởng đ−ờng chéo hoá ma trận bằng biến đổi đồng dạng.Với một ma trận đã cho A,tìm một ma trận không suy biến V,sao cho inv(V)*A*V hay ngắn gọn J = V\A*V là gần với với đ−ờng chéo.Đối với phần lớn các ma trận dạng Jordan chính tắc là ma trận đ−ờng chéo của các giá trị riêng và các cột của ma trận biến đổi là các vec tơ riêng.Một vài ma trận không đối xứng không thể đ−ờng chéo hoá.Dạng Jordan có các giá trị riêng trên đ−ờng chéo nh−ng một số phần tử phía trên là 1 thay vì 0.Phát biểu: J = jordan(A) tính ma trận chính tắc Jordan của ma trận A.Phát biểu: [V,J] = jordan(A) cũng tính biến đổi đồng dạng.Các cột của ma trận V là các vec tơ riêng tổng quát hoá của A. Ví dụ: A = sym([12,32,66,116;25,76,164,294; 21,66,143,256;6,19,41,73]) A = [ 12, 32, 66, 116] [ -25,-76, -164, -294] [ 21, 66, 143, 256] 56 [ -6, -19, -41, -73] và: [V,J] = jordan(A) tạo ra: V = [ 4, -2, 4, 3] [ -6, 8, -11, -8] [ 4, -7, 10, 7] [ -1, 2, -3, -2] J = [ 1, 1, 0, 0] [ 0, 1, 0, 0] [ 0, 0, 2, 1] [ 0, 0, 0, 2] Nh− vậy A có hai giá trị riêng là 1,với khối Jordan đơn và hai giá trị riêng là 2 cũng với khối Jordan đơn.Ma trận chỉ có 2 vec tơ riêng V(:,1) và V(:,3).Chúng thoả mãn: A*V(:,1) = 1*V(:,1) A*V(:,3) = 2*V(:,3) Hai cột khác của V là các vec tơ riêng tổng quát hoá cấp 2.Chúng thoả mãn: A*V(:,2) = 1*V(:,2) + V(:,1) A*V(:,4) = 2*V(:,4) + V(:,3) Về khái niệm toán học,với vj=v(:,j),cột của V và các giá trị riêng thoả mãn quan hệ: (A-λ2I)v4 = v3 (A- λ1I)v2 = v1 5. Phân tích giá trị kì dị:Nếu A là ma trận kí tự hay dấu chấm động thì : S = svd(A) tính các giá trị kì dị của A với độ chính xác xác định bởi lựa chọn số chữ số hiệnh hành.Và: [U,S,V] = svd(A); tạo ra 2 ma trận trực giao U và V và một ma trận đ−ờng chéo S sao cho: A = U*S*V'; Ta xem ma trận n ì n với các phần tử đ−ợc xác định bằng: A(i,j) = 1/(i-j+1/2) Có nhiều cách tạo ra ma trận A.Chúng đ−ợc mô tả sau.Với n = 5 ta có: A = [ 2, -2, -2/3, -2/5, -2/7] [ 2/3, 2, -2, -2/3, -2/5] [ 2/5, 2/3, 2, -2, -2/3] [ 2/7, 2/5, 2/3, 2, -2] [ 2/9, 2/7, 2/5, 2/3, 2] Có nhiều giá trị kì dị của ma trận này gần với π.Khi n = 16,ma trận tính svd(A) với dấu chấm động cho kết quả 3.14159265358979 3.14159265358979 3.14159265358979 3.14159265358979 3.14159265358976 57 3.14159265358767 3.14159265349961 3.14159265052655 3.14159256925492 3.14159075458606 3.14155754359918 3.14106044663470 3.13504054399745 3.07790297231120 2.69162158686066 1.20968137605669 Bốn giá trị kì dị đầu tiên xuất hiện bằng π với độ chính xác cao.Cách th−ờng thấy để tạo ma trận A là: for i=1:n for j=1:n A(i,j) = sym(1/(i-j+1/2)); end end Tuy nhiên mô hình tính toán dựa trên ma trận của MATLAB cho phép một cách khác,hiệu quả hơn gọi là “mẹo của Tony” để tạo A là: m = ones(n,1); i=(1:n)'; j=1:n; A = sym(1./(i(:,m)-j(m,:)+1/2)); Cách hiệu quả nhất để tạo ma trận này là phát biểu thuần tuý số: [J,I] = meshgrid(1:n); A = sym(1./(I - J+1/2)); Do các phần tử của ma trận A là tỉ số của hai số nguyên nhỏ nhất,vpa(A) tạo ra một biểu diễn có độ chính xác thay đổi.Do vậy: S = svd(vpa(A)) tính toán các giá trị kì dị với độ chính xác rất cao.Với n = 16 và digits(3) ta có kết quả: S = [ 1.20968137605668985332455685357 ] [ 2.69162158686066606774782763594 ] [ 3.07790297231119748658424727354 ] [ 3.13504054399744654843898901261 ] [ 3.14106044663470063805218371924 ] [ 3.14155754359918083691050658260 ] [ 3.14159075458605848728982577119 ] [ 3.14159256925492306470284863102 ] [ 3.14159265052654880815569479613 ] [ 3.14159265349961053143856838564 ] [ 3.14159265358