Читать книгу Applied Numerical Methods Using MATLAB - Won Y. Yang - Страница 21

We can define a new scalar/vector/matrix or redefine any existing ones in terms of the existent ones or irrespective of them. In the MATLAB Command window, let us define A and B as

by running

>A=[1 2 3;4 5 6], B=[3;-2;1]

We can modify them or take a portion of them. For example:

>A=[A;7 8 9] A= 1 2 3 4 5 6 7 8 9 >B=[B [1 0 -1]'] B= 3 1 -2 0 1 -1

Here, the apostrophe (prime) operator ( ') takes the complex conjugate transpose and functions virtually as a transpose operator for real‐valued matrices. If you want to take just the transpose of a complex‐valued matrix, you should put a dot ( .) before ', i.e. ‘ .'’.

When extending an existing matrix or defining another one based on it, the compatibility of dimensions should be observed. For instance, if you try to annex a 4 × 1 matrix into the 3 × 1 matrix B, MATLAB will reject it squarely, giving you an error message.

>B=[B ones(4,1)] Error using horzcat Dimensions of matrices being concatenated are not consistent.

We can modify or refer to a portion of a given matrix.

>A(3,3)=0 A= 1 2 3 4 5 6 7 8 0 >A(2:3,1:2) % from 2nd row to 3rd row, from 1st column to 2nd column ans= 4 5 7 8 >A(2,:) % 2nd row, all columns ans= 4 5 6

The colon ( :) is used for defining an arithmetic (equal difference) sequence without the bracket [] as

>t=0:0.1:2

which makes

t=[0.0 0.1 0.2 … 1.9 2.0]

1 (Q8) What if we omit the increment between the left/right boundary numbers?

2 (A8) By default, the increment is 1.>t=0:2 t= 0 1 2

3 (Q9) What if the right boundary number is smaller/greater than the left boundary number with a positive/negative increment?

4 (A9) It yields an empty matrix, which is useless.>t=0:-2 t= Empty matrix: 1-by-0

5 (Q10) If we define just some elements of a vector not fully, but sporadically, will we have a row vector or a column vector and how will it be filled in between?

6 (A10) We will have a row vector filled with zeros between the defined elements.>D(2)=2; D(4)=3 D= 0 2 0 3

7 (Q11) How do we make a column vector in the same style?

8 (A11) We must initialize it as a (zero‐filled) row vector, prior to giving it a value.>D=zeros(4,1); D(2)=2; D(4)=3 D= 0 2 0 3

9 (Q12) What happens if the specified element index of an array exceeds the defined range?

10 (A12) It is rejected. MATLAB does not accept nonpositive or noninteger indices.>D(5) Index exceeds matrix dimensions. >D(0)=1; Subscript indices must either be real positive integers or logicals. >D(1.2) Subscript indices must either be real positive integers or logicals.

11 (Q13) How do we know the size (the numbers of rows/columns) of an already‐defined array?

12 (A13) Use the ‘ length()’, ‘ size()’, and ‘ numel()’ commands as indicated as follows:>length(D) ans = 4 >[M,N]=size(A) M = 3, N = 3 >Number_of_elements=numel(A) Number_of_elements = 9

MATLAB enables us to handle vector/matrix operations in almost the same way as scalar operations. However, we must make sure of the dimensional compatibility between vectors/matrices, and put a dot ( .) in front of the operator for elementwise (element‐by‐element) operations. The addition of a matrix and a scalar adds the scalar to every element of the matrix. The multiplication of a matrix by a scalar multiplies every element of the matrix by the scalar.

There are several things to know about the matrix division and inversion.

Remark 1.1 Rules of Vector/Matrix Operation

1 (1) For a matrix to be invertible, it must be square and nonsingular, i.e. the numbers of its rows and columns must be equal and its determinant must not be zero.

2 (2) The MATLAB command ‘pinv(A)’ provides us with a matrix X of the same dimension as AT such that AXA = A and XAX = X. We can use this command to get the right/left pseudo (generalized) inverse AT[AAT]−1/[ATA]−1AT for a matrix A given as its input argument, depending on whether the number (M) of rows is smaller or greater than the number (N) of columns, so long as the matrix is of full rank, i.e. rank(A) = min(M,N) [K-2, Section 6.4]. Note that AT[AAT]−1/[ATA]−1AT is called the right/left inverse since it is multiplied onto the right/left side of A to yield an identity matrix.

3 (3) You should be careful when using the ‘ pinv(A)’ command for a rank‐deficient matrix, because its output is no longer the right/left inverse, which does not even exist for rank‐deficient matrices.

4 (4) The value of a scalar function having an array value as its argument is also an array with the same dimension.

Suppose we have defined vectors a₁, a₂, b₁, b₂, and matrices A₁, A₂, B as follows:

>a1=[-1 2 3]; a2=[4 5 2]; b1=[1 -3]'; b2=[-2 0]; , , , >A1=[a1;a2], A2=[a1;[b2 1]], B=[b1 b2'] , ,

The results of various operations on these vectors/matrices are as follows (pay attention to the error message):

>A3=A1+A2, A4=A1-A2, 1+A1 % matrix/scalar addition/subtraction A3= -2 4 6 A4= 0 0 0 ans= 0 3 4 2 5 3 6 5 1 5 6 3 >AB=A1*B % matrix multiplication? Error using * Inner matrix dimensions must agree. >BA1=B*A1 % regular matrix multiplication BA1=-9 -8 -1 3 -6 -9 >AA=A1.*A2 % element-wise (term-wise) multiplication AA= 1 4 9 -8 0 2 >AB=A1.*B %AB(m, n) = A1(m, n)B(m, n) element-wise multiplication Error using .* Matrix dimensions must agree. >A1_1=pinv(A1),A1'*(A1*A1')̂-1,eye(size(A1,2))/A1 % A1_1= -0.1914 0.1399 %right inverse of a 2x3 matrix A1 0.0617 0.0947 0.2284 -0.0165 >A1*A1_1 %A1/A1=I implies the validity of A5_1 as the right inverse ans= 1.0000 0.0000 0.0000 1.0000 >A5=A1'; % a 3x2 matrix >A5_1=pinv(A5),(A5'*A5)̂-1*A5',A5\eye(size(A5,1)) % A5_1= -0.1914 0.0617 0.2284 %left inverse of a 3x2 matrix A5 0.1399 0.0947 -0.0165 >A5_1*A5 %A5\A5=I implies the validity of A5_1 as the left inverse ans= 1.0000 -0.0000 -0.0000 1.0000 >A1_li=(A1'*A1)̂-1*A1' %the left inverse of matrix A1 with M<N? Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.211163e-18. A1_li = -0.2500 0 0.2500 -0.5000 0.5000 0.2500

1 (Q14) Does the left inverse of a matrix having rows fewer than columns exist?

2 (A14) No. There is no N × M matrix that is premultiplied on the left of an M × N matrix with M < N to yield a nonsingular matrix, far from an identity matrix. In this context, MATLAB should have rejected the above case on the ground that is singular and so its inverse does not exist. But, because the round‐off errors make a very small number appear to be a zero or a real zero appear to be a very small number (as will be mentioned in Remark 2.3), it is not easy for MATLAB to tell a near‐singularity from a real singularity. That is why MATLAB dares not to declare the singularity case and instead issues just a warning message to remind you to check the validity of the result so that it will not be blamed for a delusion. Therefore, you must be alert for the condition mentioned in Remark 1.1(2), which says that, in order for the left inverse to exist, the number of rows must not be less than the number of columns.>A1_li*A1 %No identity matrix, since A1_li isn't the left inverse ans = 1.2500 0.7500 -0.2500 -0.2500 0.5000 0.7500 1.5000 3.5000 2.5000 >det(A1'*A1) %A1 is not left-invertible for A1'*A1 is singular ans = 0

3 (cf) Let us be nice to MATLAB as it is to us. From the standpoint of promoting mutual understanding between us and MATLAB, we acknowledge that MATLAB tries to show us apparently good results to please us like always, sometimes even pretending not to be obsessed by the demon of ‘ill‐condition’ in order not to make us feel uneasy. How kind MATLAB is! But, we should be always careful not to be spoiled by its benevolence and not to accept the computing results every inch as it is. In this case, even though the matrix [A1'*A1] is singular and so not invertible, MATLAB tried to invert it and that's all. MATLAB must have felt something abnormal as can be seen from the ominous warning message prior to the computing result. Who would blame MATLAB for being so thoughtful and loyal to us? We might well be rather touched by its sincerity and smartness.

In the aformentioned statements, we see the slash( /)/backslash( \) operators. These operators are used for right/left division, respectively; B/A is the same as B*inv(A) and A\B is the same as inv(A)*B when A is invertible and the dimensions of A and B are compatible. Noting that B/A is equivalent to (A'\B')', let us take a close look at the function of the backslash( \) operator.

>X=A1\A1 % an identity matrix? X= 1.0000 0 -0.8462 0 1.0000 1.0769 0 0 0

1 (Q13) It seems that A1\A1 should have been an identity matrix, but it is not, contrary to our expectation. Why?

2 (A13) We should know more about the various functions of the backslash( \), which can be seen by typing ‘ help slash’ into the MATLAB Command window. Let Remark 1.2 answer this question in cooperation with the next case.>A1*X-A1 %zero if X is the solution to A1*X=A1? ans= 1.0e-015 * 0 0 0 0 0 -0.4441

Remark 1.2 The Function of Backslash(\) Operator

Overall, for the command ‘ A\B’, MATLAB finds a solution to the equation A*X=B. Let us denote the row/column dimension of the matrix A by M and N.

1 (1) If matrix A is square and upper/lower‐triangular in the sense that all of its elements below/above the diagonal are zero, then MATLAB finds the solution by applying backward/ forward substitution method (Section 2.2.1).

2 (2) If matrix A is square, symmetric (Hermitian), and positive definite, then MATLAB finds the solution by using Cholesky factorization (Section 2.4.2).

3 (3) If matrix A is square and has no special feature, then MATLAB finds the solution by using lower‐upper (LU) decomposition (Section 2.4.1).

4 (4) If matrix A is rectangular, then MATLAB finds a solution by using QR factorization (Section 2.4.2). In case A is rectangular and of full rank with rank( A) = min( M, N), it will be the least‐squares (LSs) solution (Eq. (2.1.10)) for M> N (over‐determined case) and one of the many solutions that is not always the same as the minimum‐norm solution (Eq. (2.1.7)) for M< N (under‐determined case). But for the case where A is rectangular and has rank deficiency, what MATLAB gives us may be useless. Therefore, you must pay attention to the warning message about rank deficiency, which might tell you not to count on the dead‐end solution made by the backslash(\) operator. To find an alternative in the case of rank deficiency, you had better resort to the singular value decomposition (SVD) (see Problem 2.8 for details).

For the moment, let us continue to try more operations on matrices.

>A1./A2 % termwise right division ans= 1 1 1 -2 Inf 2 >A1.\A2 % termwise left division ans= 1 1 1 -0.5 0 0.5 >format rat, B̂-1 %represent the numbers (of B^‐1<![CDATA[) in fractional form ans= 0 -1/3 -1/2 –1/6 >inv(B) % inverse matrix, equivalently ans= 0 -1/3 -1/2 –1/6 >B.̂-1 % element-wise inversion(reciprocal of each element) ans= 1 -1/2 -1/3 Inf >B̂2 % square of B, i.e., B²<![CDATA[=B*B ans= 7 -2 -3 6 >B.̂2 % element-wise square(square of each element) ans= 1() 4() 9() 0() >2.̂B % 2 to the power of each number in B ans= 2 () 1/4() 1/8() 1() >A1.̂A2 % element of A1 to the power of each element in A2 ans= -1 () 4() 27() 1/16() 1() 2() >format short, exp(B) %elements of eB with 4 digits below the dp ans= 2.7183() 0.1353() 0.0498() 1.0000()

There are more useful MATLAB commands worthwhile to learn by heart.

Remark 1.3 More Useful Commands for Vector/Matrix Operations

1 (1) We can use the commands ‘ zeros()’, ‘ ones()’, and eye()’ to construct a matrix of specified size or the same size as an existing matrix that has only zeros, only ones, or only ones/zeros on/off its diagonal.>Z=zeros(2,3) % or zeros(size(A1)) yielding a 2x3 zero matrix Z= 0 0 0 0 0 0 >E=ones(size(B)) % or ones(3,2) yielding a 3x2 one matrix E= 1 1 1 1 1 1 >I=eye(2) % yielding a 2x2 identity matrix I= 1 0 0 1

2 (2) We can use the ‘ diag()’ command to make a column vector composed of the diagonal elements of a matrix or to make a diagonal matrix with on‐diagonal elements taken from a vector given as the input argument.>A1=[-1 2 3; 4 5 2] A1= -1 2 3 4 5 2 >a=diag(A1) % The column vector consisting of diagonal elements a= -1 5 >diag(a) % The column vector consisting of diagonal elements ans= -1 0 0 5

3 (3) We can use the commands ‘ sum()’/‘ prod()’ to get the sum/product of elements in a vector or a matrix, column‐wisely first (along the first nonsingleton dimension).>sa1=sum(a1) % sum of all the elements in vector a1 sa1= 4 %∑a1(n) = - 1 + 2 + 3 = 4 >sA1=sum(A1) % sum of all the elements in each column of matrix A1 sA1= 3 7 5 % >SA1=sum(sum(A1)) % sum of all elements in matrix A1 SA1= 15 % >pa1=prod(a1) % product of all the elements in vector a1 pa1= -6 %∏a1(n) = (-1) × 2 × 3 = - 6 >pA1=prod(A1) % product of all the elements in each column of matrix A1 pA1= -4 10 6 % >PA1=prod(prod(A1))% product of all the elements of matrix A1 PA1= -240 %

4 (4) We can use the commands ‘ max()’/‘ min()’ to find the first maximum/minimum number and its index in a vector or in a matrix given as the input argument.>[aM,iM]=max(a2) aM= 5, iM= 2 %means that the max. element of vector a2 is a2(2)=5 >[AM,IM]=max(A1) AM= 4 5 3 IM= 2 2 1 %means that the max. elements of each column of A1 are A1(2,1)=4, A1(2,2)=5, A1(1,3)=3 >[AMx,J]=max(AM) AMx= 5, J= 2 %implies that the max. element of A1 is A1(IM(J),J)=A1(2,2)=5

5 (5) We can use the commands ‘ rot90()’/‘ fliplr()’/‘ flipud()’ to rotate a matrix by an integer multiple of 90° and to flip it left–right/up–down.>A1, A3=rot90(A1), A4=rot90(A1,-2) A1= -1 2 3 4 5 2 A3= 3 2 %90o rotation 2 5 -1 4 A4= 2 5 4 %90 ox(-2) rotation 3 2 -1 >A5=fliplr(A1) %flip left-right A5= 3 2 -1 2 5 4 >A6=flipud(A1) %flip up-down A6= 4 5 2 -1 2 3

6 (6) We can use the ‘ reshape()’ command to change the row–column size of a matrix with its elements preserved (column‐wisely first).>A7=reshape(A1,3,2) A7= -1 5 4 3 2 2 >A8=reshape(A1,6,1), A8=A1(:) %makes super-column vector A8= -1 4 2 5 3 2

7 (7) We can use the ‘ repmat(A,M,N)’ command to repeat a matrix A to make a large matrix consisting of an M× N tiling of copies of A.>A9=repmat(A1,2,3) A9= -1 2 3 -1 2 3 -1 2 3 4 5 2 4 5 2 4 5 2 -1 2 3 -1 2 3 -1 2 3 4 5 2 4 5 2 4 5 2

8 (8) We can use the ‘ find’ command to find the row/column indices of the entries satisfying the condition specified as its input argument.>[ir,ic]=find(2<A1&A1<5) ir = 2 ic = 1 1 3This means that the matrix A1 has two entries, A1(2,1) and A1(1,3), satisfying the two inequalities, ‘greater than 2’ and ‘smaller than 5’.