Ma tran va Mang trong Matlab

Chươ ng ng 2 :Ma tr ận và M ảng trong Matlab

Biên soạn: Nguyễ n Thị H ồng Thúy

Chươ ng ng 2

MA TR ẬN VÀ MẢNG TRONG MATLAB Tất c ả mọi s ự tính toán đều có một điểm chung là có s ử dụng đến các đại l ượ ng ng vô h ướ ng ng g ọi là scalars. Phép toán có liên quan đến scalars là các phép toán cơ bản, nhưng một lúc nào đó phép toán ph ải lậ p lại nhiều lần khi tính trên nhi ều số. Ðể giải quyết vấn đề này MATLAB đưa ra các khái niệm và thao tác tính toán trên mảng và ma tr ận.

2.1.

Mảng đơ n

2.1.1. Cấu trúc các mảng cơ bản - Ð ể tạo m ảng, ta

đặt các phần t ử của mảng vào giữa hai dấu ngoặc vuông, gi ữa hai phần tử của mãng có thể là dấu cách hoặc dấu phẩy. - Vớ i mảng có số lượ ng ng phần tử ít thì ta có thể nhậ p vào tr ực tiế p, nhưng vớ i mảng có số lượ ng ng lớ n các phần tử thì ta có thể dùng các cách sau: ướ c x = first : last : tạo vectơ hàng x bắt đầu tại first, phần tử sau bằng phần tử tr ướ cộng vớ i 1, k ết thúc là ph ần tử có giá tr ị bằng hoặc nhỏ hơ n last. x = first : increment : last : tạo vectơ hàng x bắt đầu tại first, giá tr ị cộng là increment, k ết thúc là ph ần tử có giá tr ị bằng hoặc nhỏ hơ n last. x = linspace (first,last,n): tạo vectơ hàng x bắt đầu tại first, k ết thúc là last, có n phần tử. x = logspace (first,last,n): tạo vectơ hàng không gian logarithm x b ắt đầu tại first last 10 , k ết thúc tại 10 , có n ph ần tử.

Ví dụ: >> a=1:5 a= 1 2 3

4

5

4

5

>> b=[6 7 8] b= 6 7 8 >> c=[a b] c= 1

2

3

6

7

8

2.1.2. Vectơ hàng và vectơ cột ướ c, Trong các ví dụ tr ướ c, mảng chứa một hàng và nhi ều cột, thườ ng ng gọi là vectơ hàng. Mảng có một cột và nhiều hàng gọi là vectơ cột. 12

Chươ ng ng 2 :Ma tr ận và M ảng trong Matlab


Ðể tạo vectơ cột, ta dùng dấu ch ấm ph ẩy

để phân cách các phần t ử. Ngoài ra, ta cũng có thể dùng các hàm linspace, logspace, hay từ các vectơ hàng, sau đó dùng ph ươ ng ng pháp chuy ển vị. Ví dụ: >> x=linspace(0,pi,5) x = 0

0.7854

1.5708

2.3562

2.1416

>> y=x' y = 0 0.7854 1.5708 2.3562 2.1416

2.2.

Ma trận

2.2.1. Nhập một ma trận trong MATLAB 2.2.1.1. Nhập một ma trận từ một danh sách tườ ng ng minh: Một ma tr ận trong MATLAB đượ c định ngh ĩ a như một mảng nhiều chiều và theo nguyên t ắc sau: Mỗi phần tử trên từng dòng c ủa ma tr ận đượ c cách nhau bở i dấu phẩy hoặc khoảng tr ống. Mỗi hàng đượ c phân cách bở i một dấu chấm phẩy. Bao quanh ma tr ận bở i một cặ p ngoặc vuông. Ví dụ: Tạo ma tr ận 3x3:  1 2 3 4 5 6   7 8 9  Lệnh MATLAB: >> A=[1 2 3;4 5 6;7 8 9] A = 1 4

2 5

3 6

7

8

9

2.1.1.2. Tạo ma trận từ nhữ ng ng hàm có sẵn trong MATLAB: MATLAB có một thư viện các hàm cho phép t ạo ma tr ận. Sau đây là một số hàm:

zeros(n,m): tạo một ma tr ận kích thướ c n x m, vớ i các phần tử đều bằng không. eye(n): tạo một ma tr ận đơ n vị kích thướ c n x n. ones(n,m): tạo ma tr ận kích thướ c n x m, v ớ i các phần tử đều bằng một. 13

Chươ ng 2 :Ma tr ận và M ảng trong Matlab


rand(n,m): tạo ma tr ận kích thướ c n x m, vớ i các phần tử có giá tr ị ngẫu nhiên từ 0 -1

diag(V): nếu V là một vectơ sẽ tạo ra ma tr ận đườ ng chéo, vớ i các phần t ử của vectơ V nằm trên đườ ng chéo. Ví dụ: >> zeros(2,3) ans = 0 0

0 0

0 0

>> diag([1 2 3]) ans = 1 0 0

0 2 0

0 0 3

2.1.1.3. Nhập ma trận từ một file Lệnh Load trong MATLAB dùng để đọc file chứa ma tr ận t ạo ra từ những l ệnh MATLAB tr ướ c đó. L ệnh này còn dùng để đọc text file chứa nh ững d ữ liệu s ố. Text file phải đượ c t ổ chức như một bảng số mà các phần tử đượ c cách nhau bở i các khoảng tr ống, mỗi hàng của ma tr ận chiếm mỗi hàng của text file. Số phần tử của mỗi hàng phải bằng nhau. 2.1.1.4. Sử dụng M-file Trên màn hình MATLAB, ch ọn File → New → M-file. Trên màn hình so ạn thảo nhậ p vào các l ệnh sau: A=[1 2 3 4 5 6 7 8 9] Save file vớ i tên matran. Sau đó từ màn hình MATLAB đánh vào: matran K ết quả là: A = 1 2 3 4 7

5 8

6 9

2.2.2. Các thao tác đối vớ i ma trận 2.2.2.1. Sự móc nối ma trận MATLAB cho phép k ết hợ p các ma tr ận con để tạo nên một ma tr ận lớ n hơ n.

Ví dụ: >> b=ones(3,3); >> c=zeros(3,3); >> a=[b c;c b] 14



a= 1 1 1 0

1 1 1 0

1 1 1 0

0 0 0 1

0 0 0 1

0 0 0 1

0 0

0 0

0 0

1 1

1 1

1 1

2.2.2.2. Xóa dòng và cột của ma trận MATLAB cho phép xóa dòng ho ặc cột của ma tr ận bằng cách gán các giá tr ị r ỗng cho hàng ho ặc cột của ma tr ận. Một giá tr ị r ỗng đượ c kí hiệu bở i [].

Ví dụ: >> a=[1 2 3;4 5 6;7 8 9] a= 1 4

2 5

3 6

7 8 9 >> a(2,:)=[] a= 1 2 3 7 8 9 >> a(:,3)=[] a= 1

2

7

8

2.2.2.3. Ma trận chuyển vị Ma tr ận chuyển v ị của ma tr ận A là một ma tr ận mà các hàng của ma tr ận A là các cột của ma tr ận này.

Ví dụ: >> a=[1 2 3;4 5 6;7 8 9] a= 1 4

2 5

3 6

7 8 >> b=a' b=

9

1 2 3

4 5 6

7 8 9 15



2.2.2.4. Lệnh diag Dùng để tạo ma tr ận đườ ng chéo và rút ra đườ ng chéo của ma tr ận Cú pháp: diag(v,k): nếu v là một vectơ gồm n phần t ử thì k ết quả là một ma tr ận vuông bậc n + |k|. Trong đó các phần tử của v nằm trên đườ ng chéo thứ k. Nếu k = 0, đườ ng chéo là đườ ng chéo chính, k > 0 là đườ ng chéo thứ k nằm trên đườ ng chéo chính, k < 0 là đườ ng chéo thứ k nằm dướ i đườ ng chéo chính.

diag(X,k): nếu X là một ma tr ận thì k ết quả là một vectơ cột hình thành từ những phần tử của đườ ng chéo thứ k. diag(X): tr ả về một vectơ là đườ ng chéo chính c ủa ma tr ận. diag(diag(X)): tr ả về một ma tr ận đườ ng chéo. Ví dụ: >> v=[1 2 3]; >> diag(v,1) ans = 0 0 0

1 0 0

0 2 0

0 0 3

0 0 0 0 >> X=[1 2 3;4 5 6;7 8 9]; >> diag(X,0) ans = 1 5 9

2.2.2.5. Lệnh sum Tính tổng các hàng hay các cột của ma tr ận. Cú pháp: sum(X) hay sum(X,1): tr ả về một vectơ mà mỗi phần tử là tổng của từng cột trong ma tr ận. sum(X,2): tr ả về một vectơ mà mỗi phần tử là tổng của từng hàng trong ma tr ận. Ví dụ: >> a=[1 2 3;4 5 6;7 8 9] a= 1 2 3 4

5

6 16


7

8


9

>> tong_cot=sum(a) tong_cot = 12 15 18 >> tong_hang=sum(a,2) tong_hang = 6 15 24

2.2.2.6. Ma trận symbolic Có hai cách định ngh ĩ a một ma tr ận symbolic: - Từ các tham số - Từ các số thực Ðể định ngh ĩ a ma tr ận symbolic, hai l ệnh sym và syms th ườ ng đượ c sử dụng: sym(‘a’): tr ả về k ết quả là một biến symbolic tên là a. sym(‘[…;…;…]’):tr ả về một ma tr ận symbolic. sym(A): vớ i A là một số thực hay ma tr ận số thực sẽ tr ả về một biến hay ma tr ận symbolic.

syms arg1 arg2 tươ ng đươ ng vớ i arg1 = sym(‘arg1’) ; arg2 = sym(‘arg2’) 2.2.2.7. Lệnh det Dùng để tính định thức của ma tr ận. Cú pháp: det(A): k ết quả là biểu thức symbolic nếu A là ma tr ận symbolic, là một giá tr ị số nếu A là một ma tr ận số. Ví dụ: >> syms a b c d >> A=[a b;c d]; >> r=det(A) r= a*d-b*c Ðối vớ i một số ma tr ận đặc biệt ta có một số k ết quả sau: - Ðịnh thức của một ma tr ận đơ n vị bằng một. - Ðịnh thức của một ma tr ận đườ ng chéo đơ n giản là tích của của các phần tử đườ ng chéo. Một ma tr ận mà định thức của nó có giá tr ị bằng không ng ườ i ta gọi đó là ma tr ận suy bi ến. Ngoài vi ệc dùng định thức để giải hệ phươ ng trình tuyến tính, ng ườ i ta còn dùng nó

để xác định điều kiện có nghi ệm hay không c ủa hệ. 17



2.2.2.8. Các toán hạng ma trận Trong MATLAB t ồn tại các toán hạng sau: A+B

Cộng ma tr ận A và B. A và B ph ải có cùng kích th ướ c, ngoại tr ừ một trong hai là một giá tr ị vô hướ ng.

A–B

Tr ừ ma tr ận A và B. A và B ph ải có cùng kích thướ c, ngoại tr ừ một trong hai là một giá tr ị vô hướ ng.

A*B

Nhân ma tr ận A và B. S ố cột của ma tr ận A phải bằng số hàng của ma tr ận B, ngoại tr ừ một trong hai là một giá tr ị vô hướ ng.

Nhân từng phần t ử của ma tr ận A v ớ i t ừng phần tử của ma tr ận B. K ết quả là A .* B một ma tr ận. A và B ph ải có cùng kích th ướ c, ngoại tr ừ một trong hai là một giá tr ị vô hướ ng Chia trái ma tr ận. X = A\B t ươ ng

A\B

đươ ng v ớ i vi ệc gi ải h ệ phươ ng trình tuyến

tính A*X = B

A .\ B

Chia trái mảng. A .\ B t ươ ng đươ ng vớ i B(i,j)/A(i,j). A và B phải có cùng kích thướ c, ngoại tr ừ một trong hai là m ột giá tr ị vô hướ ng.

A/B

Chia phải ma tr ận. X = A/B tươ ng đươ ng vớ i việc giải hệ phươ ng trình tuyến tính B*X = A.

A ./ B

Chia phải mảng. A ./ B t ươ ng đươ ng vớ i A(i,j)/B(i,j). A và B ph ải có cùng kích thướ c, ngoại tr ừ một trong hai là một giá tr ị vô hướ ng.

đều là ma tr ận. Lũy th ừa mảng. K ết qu ả là một ma tr ận mà các số hạng là A(i,j)^B(i,j). A và A .^ B B phải có cùng kích thướ c, ngoại tr ừ một trong hai là m ột giá tr ị vô hướ ng. A^B

2.3.

Lũy thừa ma tr ận. Lỗi sẽ phát sinh nếu A và B

Gải hệ phươ ng trình tuyến tính Một hệ phươ ng trình tuyến tính có dạng tổng quát sau: a11x1 + a12x2 + a1nxn = b1 a21x1 + a22x2 + a2nxn = b2 M

M

am1x1 + am2x2 + amnxn = bm Vớ i: A = [aij]m x n là ma tr ận hệ số. *

A = [A b]m x (n +1) là ma tr ận đầy đủ. Một số phươ ng pháp để giải hệ này: - Nghịch đảo ma tr ận - Phươ ng pháp kh ử Gauss - Phươ ng pháp kh ử Gauss - Jordan - Phươ ng pháp phân rã ma tr ận (LU) Một trong những ứng dụng của MATLAB là để giải hệ phươ ng trình đại số tuyến tính. Trong MATLAB có m ột số hàm đã phươ ng pháp này. 18

đượ c xây dựng để sử dụng cho các



2.3.1. Nghịch đảo ma trận Xét hệ phươ ng trìng tuyến tính. Dướ i dạng ma tr ận hệ có dạng sau: AX = B ⇒

-1

X=A B

-1

Vớ i A là ma tr ận nghịch đảo của ma tr ận hệ số A.

 a 11 a 12 a a 21 22 A=  M M  a n1 a n2

L a 1n L a 2n M

O

K a nn

     

 x1  x   2 X=  M x n 

 b1   b   2 B=  M  b n 

2.3.1.1. Lệnh inv inv(A): dùng để tính ma tr ận nghịch đảo. Ví dụ: giải hệ phươ ng trình tuyến tính A*X = B sau: 1 1  A= 1 3  2 1

- 2

 x1    X = x 2  x   3

  - 5 -1

 4   B = 7  7   

>> A=[1 1 -2;1 3 -1;2 1 -5]; >> B=[4;7;7]; >> A_inv=inv(A) A_inv = 14.0000 -3.0000 -5.0000 -3.0000 1.0000 1.0000 5.0000 -1.0000 -2.0000 >> X=A_inv*B X = 0 2.0000 -1.0000

2.3.1.2. Lệnh pinv pinv(A): dùng để tính giả nghịch đảo c ủa ma tr ận m x n, vớ i m ≠ n. L ệnh pinv không sử dụng đượ c vớ i phươ ng pháp symboic. Ví dụ: giải hệ phươ ng trình tuyến tính A*X = B sau: 1 2 A= 0  2

4

-2

-1

5

4

-2

0

-2

-6

1

0

0

3

-2

1

4

-8

19

  0  5 2

 x1  x   2  x 3  X=  x 4  x 5    x 6 

− 6  − 2    B=  0  − 3



>> A=[1 -2 -1 5 4 4;2 -2 0 -2 -6 2;0 1 0 0 3 0;2 -2 1 4 -8 5]; >> B=[-6;-2;0;-3]; >> X=pinv(A)*B X = -0.4031 0.8008 0.3984 -0.0157 -0.2669 -0.6126 Vì A không ph ải là ma tr ận vuông nên m ột thông báo l ỗi sẽ hiện ra khi ta thay lệnh inv(A) bằng lệnh pinv(A). Có thể giải lại hệ phươ ng trình tuyến tính trên bằng phươ ng pháp symbolic: >> A=sym([1 -2 -1 5 4 4;2 -2 0 -2 -6 2;0 1 0 0 3 0;2 -2 1 4 -8 5]); >> B=sym([-6;-2;0;-3]); >> X=A\B X = [ 0] [ 0] [ 2] [ 0] [ 0] [ -1]

2.3.1.3. Lệnh invhilb Dùng để tính nghịch đảo ma tr ận Hilbert. Cú pháp: invhilb(n) n là kích thướ c của ma tr ận Hilbert. Ma tr ận Hilbert đượ c xem là điều ki ện y ếu “ill condition”, có ngh ĩ a là định thức của ma tr ận có giá tr ị r ất nhỏ. Ma tr ận Hilbert n xn có d ạng như sau: 1 1/2  H = 1/3   M 1/n

1/2

1/3

K

1/3

1/4

L

1/4

1/5

L

M

M

O

1/(n + 1)

1/(n + 2)

L

  1/(n + 1)  1/(n + 2)  M  1/(2n)  1/n

Ví dụ: Giải phươ ng trình H*X = B b ằng hai phươ ng pháp s ố, symbolic và dùng l ệnh invhilb. Trong đó H là một ma tr ận Hilbert 8 x 8 và B là vect ơ cột [1 1 1 1 1 1 1 1]. 20



Phươ ng pháp số: >> H=hilb(8) H = 1.0000 0.5000

0.5000 0.3333

0.3333 0.2500

0.2500 0.2000

0.2000 0.1667

0.1667 0.1429

0.1429 0.1250

0.1250 0.1111

0.3333 0.2500

0.2500 0.2000

0.2000 0.1667

0.1667 0.1429

0.1429 0.1250

0.1250 0.1111

0.1111 0.1000

0.1000 0.0909

0.2000 0.1667

0.1667 0.1429

0.1429 0.1250

0.1250 0.1111

0.1111 0.1000

0.1000 0.0909

0.0909 0.0833

0.0833 0.0769

0.1429 0.1250 0.1250 0.1111 >> B=ones(8,1)

0.1111 0.1000

0.1000 0.0909

0.0909 0.0833

0.0833 0.0769

0.0769 0.0714

0.0714 0.0667

B = 1 1 1 1 1 1 1 1 >> X=inv(H)*B X = 1.0e+005 * -0.0001 0.0050 -0.0756 0.4620 -1.3860 2.1622 -1.6817 0.5148 Phươ ng pháp symbolic: >> H=sym(hilb(8)) H = [ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8] [ 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9] [ 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10] [ 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11] 21



[ 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12] [ 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13] [ 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14] [ 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15] >> B=ones(8,1); >> X=inv(H)*B X = [ -8] [ 504] [ -7560] [ 46200] [ -138600] [ 216216] [ -168168] [ 51480] Dùng l ệnh invhilb: >> B=ones(8,1); >> X=invhilb(8)*B X = -8 504 -7560 46200 -138600 216216 -168168 51480

2.3.2. Phươ ng pháp khử Gauss - Jordan rref(A): tr ả về ma tr ận là bướ c cuối cùng trong phươ ng pháp kh ử Gauss - Jordan. Trong đó A là ma tr ận vuông hay ch ữ nhật. Lệnh rref cho phép sử dụng vớ i phươ ng pháp symbolic. Ví dụ: giải hệ phươ ng trình tuyến tính: 1 0  0  1

-2

-2

1

-1

1

-1

2

0

0  x 1 

18       3 x 2  18  . =   1  x 3  54  4  x 4  36

>> A=[1 -2 -2 0;0 1 -1 3; 0 1 -1 1;1 2 0 4]; 22



>> B=[18;18;54;36]; >> G_J=rref([A B]) G_J = 1 0 0 0 30 0 1 0 0 39 0 0

0 0

1 0

0 -33 1 -18

>> X=G_J(:,length(A)+1) X = 30 39 -33 -18 Khi sử dụng phươ ng pháp kh ử Gauss - Jordan s ẽ dẫn t ớ i một bất tiện là ta phải tiến hành l ại từ đầu thủ tục Gauss - Jordan cho t ừng vectơ cột B. Một phươ ng pháp cho phép tiết kiệm đượ c số lần tính toán mà phươ ng pháp phân rã ma tr ận.

đạt cùng hiệu quả là dùng

2.3.3. Phươ ng pháp phân rã ma trận [L,u] = lu(A): tr ả về ma tr ận tam giác trên U, ma tr ận tam giác dướ i L. Phân rã ma tr ận A thành các ma tr ận tam giác: A = L*U L: ma tr ận tam giác dướ i cỡ n x n, các ph ần tử đườ ng chéo đều bằng 1. U: ma tr ận tam giác trên. α 11 α 21 L=  M  α n1

0

L

α 22

L

M

O

α n2

K

0 

  M   α nn  0

và

 β 11 0 U=  M  0

β 1n 

β 12

L

β 22

L

β 2n 

M

O

M

0

L

   β nn 

Như vậy hệ phươ ng trình đượ c viết lại như sau:

⇒

A*X = B

(LU)*X = B

Ðặt U*X = Y thì: A*X = B

⇔

L*Y = B  U * X = Y

Cả hai phươ ng trình trong h ệ

đều d ễ dàng tìm ra nghi ệm vì các ma tr ận L và U

đều ở dạng tam giác. Bằng cách thế ngượ c một lần nữa để tìm X. Như vậy nghiệm của hệ A*X =B là: X = U\(L\B)

23



Ví dụ: giải hệ phươ ng trình A*X = B, trong đó: 3 - 3 A= 6  - 9

-7

-2

5

1

-4

0

5

-5

2

  - 5  12 0

 - 9 5 B=  7   11

>> A=[3 -7 -2 2;-3 5 1 0;6 -4 0 -5;-9 5 -5 12]; >> B=[-9;5;7;11]; >> [L,U]=lu(A) L = -0.3333 1.0000 0 0.3333 -0.6250 -0.1304 -0.6667 1.0000 U= -9.0000

0.1250 0

1.0000 0 0

0 1.0000 0

5.0000 -5.0000 12.0000

0 -5.3333 -3.6667 6.0000 0 0 -2.8750 2.2500 0 0 0 0.0435 >> X=U\(L\B) X = 3.0000 4.0000 -6.0000 -1.0000

2.4.

Hạng của ma trận và điều kiện có nghiệm của hệ A*X = B Hạng của ma tr ận A là số hàng khác không có trong d ạng rút gọn của A. Kí hi ệu: r (A) Ðiều kiện có nghi ệm của hệ phươ ng trình tuyến tính A*X = B, có n ẩn số: - r (A) = r (A*) = n thì h ệ có nghi ệm duy nh ất. - r (A) = r (A*) < n thì h ệ có vô số nghiệm phụ thuộc n – r (A) tham s ố. - r (A) ≠ r (A*): không t ồn tại lờ i giải của hệ phươ ng trình A*X = B. Trong toolbox c ủa MATLAB có một số lệnh liên quan

đến hạng của một ma

tr ận, không gian c ơ sở của ma tr ận.

2.4.1. Lệnh rank rank(A): tr ả về một số nguyên là hạng của ma tr ận A. Ví dụ: Xét điều kiện có nghi ệm của các hệ phươ ng trình tuyến tính sau:

24


1 1 a)  1  2

1

1 

2

-1

0 1

  3   4 


1  x1  - 1     x 2  =   x  3  3  4 

>> A=[1 1 1;1 2 -1;1 0 3;2 1 4]; >> B=[1;-1;3;4]; >> r_A=rank(A) r_A = 2 >> r_AB=rank([A B]) r_AB = 2

⇒ r_A = r_AB = 2 < n =3 nên h ệ tồn tại vô số nghiệm. 1  b) 0  1

1 1 2

0  x 1 

  1 1

1     x 2  =  2  x  - 2  3  

>> A=[1 1 0;0 1 1;1 2 1]; >> B=[1;2;-2]; >> r_A=rank(A) r_A = 2 >> r_AB=rank([A B]) r_AB = 3

⇒ r_A = 2 < r_AB =3 nên h ệ tồn tại vô nghi ệm. 1  c) 1  3

2 -1 -6

5

- 9     3 = 2  25 - 1  

>> A=[1 2 5;1 -1 3;3 -6 -1]; >> B=[-9;2;25]; >> r_A=rank(A) r_A = 3 >> r_AB=rank([A B]) r_AB = 3

⇒ r_A = r_AB = n =3 nên h ệ có nghiệm duy nh ất. 25



2.4.2. Lệnh null null (A): tr ả về ma tr ận r ỗng R (n x 0) n ếu ma tr ận A không suy bi ến Ví dụ: Xét hai ma tr ận magic 3 x 3 và 4 x 4 >> magic(3) ans = 8 1 6 3 4

5 9

7 2

>> null(magic(3)) ans = Empty matrix: 3-by-0 >> magic(4) ans = 16 2 5 11

3 13 10 8

9 7 6 12 4 14 15 1 >> null(magic(4)) ans = 0.2236 0.6708 -0.6708 -0.2236 Ta có thể dùng lệnh det để kiểm tra lại.

2.4.3. Lệnh clospace colspace (A): nếu A là ma tr ận symbolic kích thướ c n x n, không suy bi ến. K ết n quả tr ả về là một ma tr ận mà các cột là vectơ cơ sở của không gian R . Ví dụ: >> colspace(sym(magic(3))) ans = [ 1, 0, 0] [ 0, 1, 0] [ 0, 0, 1] >> colspace(sym(magic(4))) ans = [ [ [ [

1, 0, 0, 1,

0, 0] 1, 0] 0, 1] 3, -3]

26


2.5.


Ðộ chính xác của lờ i giải Xác định sai số của hệ phươ ng trình tuyến tính A*X = B: Gọi x và x e là nghiệm gần đúng và nghi ệm chính xác của hệ phươ ng trình tuyến tính A*X = B. Ta có: A x − b = δ b x − x e = δ x Gọi: δ b là chuẩn của vectơ cột δ b.

b là chuẩn của vectơ cột b. δ x là chuẩn của vectơ cột δx.

x là chuẩn của vectơ cột x. K là điều kiện (condition) của ma tr ận A, kí hi ệu là cond (A). Gi ữa K và ma tr ận A có mối quan hệ: K = A −1 A Trong đó: A−1 và A là chuẩn của ma tr ận A và ma tr ận nghịch

đảo của ma

tr ận A Sai số của hệ phươ ng trình A*X = B, 1 δ b K b

≤

δ x x

≤ K

δ x

x

có thể đượ c đánh giá như sau:

δ b b

Sai số xấ p xỉ = K*ε Vớ i ε là

độ chính xác của máy tính. Trong MATLAB có một bi ến đặc bi ệt dùng để định ngh ĩ a dộ chính xác này là eps. Ðộ chính xác này dùng để phân biệt hai giá tr ị r ất g ần nhau. Ví dụ b và b’ khác nhau khi và ch ỉ khi | |b| - |b’| | ≥ ε MATLAB cung c ấ p cho chúng ta hai hàm để tính điều kiện và chuẩn của ma tr ận. cond (A): tr ả về một giá tr ị là điều kiện của ma tr ận A. norm (X): tr ả về chuẩn của X. Nếu X là vectơ sẽ tr ả về chiều dài của nó. Nếu X là ma tr ận tr ả về giá tr ị căn bậc hai của tổng bình ph ươ ng các số hạng của ma tr ận. Hàm norm không làm vi ệc vớ i các biến symbolic.

Ví dụ: xét sự ảnh hưở ng của sai số trong hệ phươ ng trình tuyến tính H*X = B. H là ma tr ận hilbert n x n, v ớ i n = 6:13 và B là ma tr ận c ột vớ i tất cả các phần t ử đều bằng 1. Lệnh của MATLAB: for n=3:13 27



H=hilb(n);

% tao ma tran hilbert

HI=invhilb(n); B=ones(n,1); x=H\B; nx=norm(x);

% ma tran nghich dao cua hilbert % tao ma tran B % tinh nghiem cua he % tinh chuan cua vecto x

dx=norm(x-HI*B); db=norm(H*x-B);

% tinh chuan cua vecto dx % tinh chuan cua vecto db

nb=norm(B); K=cond(H);

% tinh chuan cua vecto B % tinh dieu kien cua ma tran hilbert

err=dx/nx; format short g

% sai so that cua he phuong trinh

err1=K*eps; format short g err2=K*db/nb;

% sai so xap xi % sai so lon nhat

format short g disp([n err err1 err2]) end K ết quả: Sai số thật

n

2.6.

Sai số xấ p xỉ

Sai số lớ n nhất

3

9.9722e-016

1.1636e-013

2.6873e-013

4

1.0284e-013

3.4447e-012

4.7732e-011

5

2.6061e-012

1.0583e-010

5.2463e-009

6

1.541e-010

3.3198e-009

7.9498e-007

7

4.784e-009

1.0555e-007

0.00017924

8

1.731e-007

3.3879e-006

0.038709

9

6.8192e-006

0.0001095

8.7484

10

0.00019489

0.0035583

1270.3

11

0.0055945

0.11599

1.8473e+005

12

0.091251

3.9846

4.7791e+007

13

1.8316

833.64

1.4563e+010

Giá trị riêng và vectơ riêng của ma trận Ðịnh ngh ĩ a: Cho ma tr ận vuông A c ấ p n. Số λ đượ c gọi là một giá tr ị riêng của ma tr ận A nếu n tồn tại vectơ cột x ≠ 0, x ∈ R , sao cho Ax = λx. Khi đó vectơ x

đượ c gọi là

vectơ riêng của ma tr ận A ứng vớ i giá tr ị riêng λ. Cách tìm: Vớ i A là ma tr ận vuông c ấ p n đã cho thì vectơ cột x ≠ 0 là vectơ riêng của ma tr ận A khi và ch ỉ khi Ax = λx hay: 28



Ax - λx = 0 (I là ma tr ận đơ n vị có kích thướ c n x n)

x = Ix

⇒ ⇔

Ax - λIx = 0 (A - λI)x = 0

Giá tr ị riêng λ là nghiệm của định thức: |A - λI| = 0 Nếu khai triển định thức ta sẽ đượ c một đa thức vớ i biến λ. Ða thức này đượ c gọi là đa thức đặc tr ưng. Nếu A là ma tr ận n x n thì đa thức này có dạng sau: P(λ) = an λ + an-1λ n

n-1

2

+ …+ a2x + a1x + a0

2.6.1. Lệnh eig Tính giá tr ị riêng và vectơ riêng của ma tr ận vuông, s ử dụng phươ ng pháp số và symbolic.

đượ c cho cả hai

Cú pháp:

d = eig(A): tr ả về ma tr ận d mà các giá tr ị riêng là các phần t ử nằm trên đườ ng chéo chính.

[V,D] = eig(A): tr ả về vectơ riêng chứa trong ma tr ận V và giá tr ị riêng chứa trong ma tr ận D. Ví dụ: >> A=[3 2/3;2/3 2]; >> A=sym([3 2/3;2/3 2]); >> d=eig(A) d= [ 5/3] [ 10/3] >> [V,D]=eig(A) V= [ 2, 1] [ 1, -2] D = [ 10/3, [

0]

0, 5/3]

2.6.2. Lệnh poly Tr ả về đa thức của ma tr ận A. Cú pháp:

poly(A): nếu A là một ma tr ận symbolic, thì k ết quả tr ả về là một đa thức đặc tr ưng. Nếu A là một hàm ma tr ận s ố thì k ết qu ả tr ả về là m ột m ảng chứa các hệ số của đa thức này. Ví dụ: >> A=sym([3 2/3;2/3 2]); >> P=poly(A) 29



P = x^2-5*x+50/9 >> d=solve(P) d= [ 5/3] [ 10/3]

2.6.3. Tính định thứ c, nghịch đảo và lũy thừ a của ma trận thông qua ma trận giá trị riêng và vectơ riêng MATLAB cung c ấ p một số hàm để tính định thức, nghịch đảo và lũy thừa của ma tr ận như det, inv, expm. Ngoài ra, căn cứ vào các tính chất c ủa ma tr ận, ta có thể tính toán thông qua các biểu thức: |A| = |D| -1

-1

T

A = VD V A D T e = Ve V n n T A = VD V

Vớ i D và V là ma tr ận các giá tr ị riêng và các vectơ riêng đã đượ c chuẩn hóa của ma tr ận A. Vì ma tr ận các giá tr ị riêng là ma tr ận đườ ng chéo, có ngh ĩ a là các phần tử khác không chỉ nằm trên đườ ng chéo chính. Do đó, việc thực hiện nghịch

đảo ma tr ận, tính định thức ma tr ận và lấy lũy thừa ma tr ận đượ c đơ n giản đi r ất nhiều nếu chúng ta thực hiện trên ma tr ận đườ ng chéo (chỉ thực hiện trên các phần tử đườ ng chéo) λ 1 0 V =  M  0

0

L

λ 2

L

M

O

0

L

0

  M   λ n  0

Nên: V = λ 1λ 2 λ 3 L λ n

λ 1−1  0 −1 V =   M  0 e λ  0 V e = M  0

1

0

L

0

λ -21

L

0

M

O

0

L



  λ -n1  M

0

L

0 

e λ

L

0 

M

O

0

L

2



  e λ  M

n

30


λ 1n  0 n V =  M  0


0

L

λ n2

L

M

O

0

L

0

  M   λ nn  0

→

→

Khái niệm về chuẩn hóa các vectơ riêng: vectơ chuẩn hóa xc của x bất k ỳ bằng: → →

xc =

x →

x →

→ →

Trong đó: x = x . x là chiều dài của vectơ x.

eig của MATLAB không cung c ấ p cho chúng ta ma tr ận vectơ riêng đượ c chuẩn hóa. Do đó, để sử dụng các tính chất của ma tr ận ở trên chúng ta cần chuẩn hóa chúng tr ướ c. Hàm

Ví dụ: A=[3 2/3;2/3 2]; [v,d]=eig(A) vc=v/norm(v) det_A=det(d) ham_det=det(A) d1=d; d1(1,1)=1/d(1,1); d1(2,2)=1/d(2,2); inv_A=vc*d1*vc' ham_inv=inv(A) d1(1,1)=exp(d(1,1)); d1(2,2)=exp(d(2,2)); exp_A=vc*d1*vc' ham_exp=expm(A) v= 0.4472 -0.8944 -0.8944 -0.4472 d= 1.6667 0 0

3.3333

vc = 0.4472 -0.8944 -0.8944 -0.4472 det_A = 31



5.5556 ham_det = 5.5556 inv_A = 0.3600 -0.1200 -0.1200 0.5400 ham_inv = 0.3600 -0.1200 -0.1200 0.5400 exp_A = 23.4842

9.0949

9.0949 9.8419 ham_exp = 23.4842 9.0949 9.0949

9.8419

Ví dụ: Tính tần số riêng và vectơ riêng của hệ dao động đều hòa:

Biết: k2=40 N/m

k1=10 N/m

m1 = m2 = 1 kg Tại vị trí ban đầu: x1(0) = 0.1, x2 = 0.2 Phươ ng trình dao động của hai khối lượ ng: ..

..

m1 x1 = − k 1 x1 − k 2 ( x1 − x2 )

m x1 = −(k 1 + k 2 ) x1 + k 2 x 2

⇔

..

..

m2 x2 = k 2 x1 − (k 2 + k 1 ) x2 ..

x1 =

− (k 1 + k 2 )

x1 +

m x 2 = k 2 x1 − ( k 2 + k 1 ) x2

k 2

x 2 m m .. k (k + k 1 ) m x 2 = 2 x1 − 2 x 2 m m Dướ i dạng ma tr ận các phươ ng trình này đượ c viết lại như sau:

 − (k 1 + k 2 )  x..   m  .. 1  =  k 2  x 2     m 

  m − (k 1 + k 2 )   m  k 2

 x1   x   2 32

⇔

Chươ ng 2 :Ma tr ận và M ảng trong Matlab →

&& = − Hay: x


K → x M

Vì hệ dao động điều hòa nên: →

→

2 && = −ω x x

⇒

→

− ω 2 x = −

K → x ⇒ M

→ K → x = λ x M

(Ðặt λ = ω , vớ i ω là tần số dao động riêng của hệ). 2

− 50 = M  40

K

40 

− 50

Dùng MATLAB ta gi ải hệ trên như sau: >> K_M=[-50 40;40 -50]; >> [V,D]=eig(K_M) V= 0.7071 -0.7071 D =

0.7071 0.7071

-90 0 0 -10

33

→

ω 2 x =

K → x M

⇒

Ma tran va Mang trong Matlab

Recommend Documents