Introduction to Matlab |
Although Matlab is a useful calculator, its main power is that it gives a simple way of working with matrices and vectors. Indeed we have already seen how vectors are used in graphs.
Remember to keep typing in the commands as they appear here, and observe and understand the Matlab response. If you mistype, it is easy to correct using the arrows and the [Del] key. Try to use the help facility to find out about unfamiliar commands. Otherwise ask another student or a tutor.
Most people will have seen systems of equations from school. For example, we
may need to find x1, x2, and x3 so that
| |||||||||||||||
In first year Mathematics the problem is rewritten in matrix-vector
notation. We introduce a matrix A and a vector b by
|
|
% Set up a system
A = [ 1 2 -1; -2 -6 4 ; -1 -3 3 ] % of equations.
b = [ 1; -2; 1 ]
x = A\b % Find x with A x = b.
Matlab should give the solution
|
A sound idea from manual computations is to substitute the computed solution back into the system, and make sure all equations are indeed satisfied. In Matlab we can do this by checking that the matrix vector product Ax gives the vector b, or better still, we check that Ax-b is exactly the zero vector.
A*x % Check that Ax and b are the same.
b
A*x - b ; % So that we do not miss small differences,
% it is better to check that A x - b = 0 .
The vector Ax-b is known as the residual.
(a.) Change the middle element of A from -6 to 5 (i.e. a22 = 5.) What is the solution to this new system? What is the new residual? Why is the residual not exactly 0? (Hint: In Matlab the number 1.23×10-5 is displayed as 1.23e-005. The column vector [1.23×10-5;4.44×10-6] could be displayed as
| 1.0e-004 * | ||||
| 1.2300e-005 | ||||
| 0.1234 | or | 4.4400e-006 | ||
| 0.0444 |
The user can control which form is used by the FILE® PREFERENCES® NUMERIC FORMAT item on the Matlab Menu bar.)
(b.) Suppose the middle element is changed from -6 to -5. Can Matlab solve this third system? Explain what has happened.
Much scientific computation involves solving very large systems of equations with many millions of unknowns. This section gives practice with the commands that are needed to work larger matrices.
The following code sets up a 4×4 matrix, A, and then finds some of its important properties. When [ ] is used to set up matrices, either blanks or commas (,) can be used to separate entries in a row. Semicolons (;) are used to begin a new row.
A = [1 2 3 4 ; 1 4 9 16 ; 1 8 27 64 ; 1 16 81 256 ] A' det(A) eig(A) inv(A)
Even in the simple 4×4 case, it is tedious to evaluate determinants and inverses. But in Matlab they can be quickly calculated and used.
Indices can be used to show parts of a larger matrix. For example try
A(2,3), A(1:2,2:4), A(:,2), A(3,:)
In general A( i:j, k:l ) means the square sub-block of A between rows i to j and columns k to l . The ranges can be replaced by just `:' if all rows or columns are to be included.
(a.) How would you display the bottom left 2×3 corner of A? How would you find the determinant of the upper left 3×3 block of A?
(b.) Find the Matlab command to calculate the row echelon form of a matrix. The row echelon form is never useful in practice, but it is handy for checking homework in other subjects!!
Matlab also provides quick ways to create special matrices and vectors.
c = ones(4,3)
d = zeros(20,1)
I = eye(5)
D = diag( [2 1 0 -1 -2] )
L = diag( [1 2 3 4], -1 )
R = randn(5,5) % R is a 5 x 5 matrix
% of random numbers.
(More information on each of these commands can be found with help.)
Matrix-matrix and matrix-vector multiplication work as expected, provided the dimensions agree. Matrices and vectors can both be multiplied by scalars. In the next example, B is another 4×4 matrix and c is a column vector of length 4 (i.e. a 4×1 matrix).
B = [1 1 0 0
0 2 1 0
0 0 3 1
0 0 0 4 ] % Rows can be separated by
% making a new line as well as using `;`.
c = [1; 0; 0; -1]
5*B % Multiply scalars and matrices.
B*c % Multiply matrices and vectors.
A*B % Matrix by matrix multiplication.
B*A % Note: AB is not the same as BA !
Some of the following commands do not work. There is no harm in trying them.
c*c % Wrong dimensions for matrix multiplication. c*A % Wrong dimensions for matrix multiplication. c' % The transpose of c, i.e. a row vector. c'*A % Now matrix multiplication is permitted. c*c' % This multiplies a 4 x 1 by a 1 x 4 matrix. c'*c % What is this quantity usually called?
(a.) Calculate inv(A)*A and A*inv(A). Do they give the expected result? (Use the help command if you can't guess what the inv command does.)
(b.) Verify for the matrices A and B above that (AB)-1 = B-1A-1.
(c.) Verify that the Matlab command A^(-1) can also be used to form the inverse of A.
Consider again the problem of solving the system of equations Ax = b. One obvious way that appeals to Mathematicians is to calculate A-1b. As in the previous 3x3 case, we should also check that the solution is correct. We do this by calculating the residual Ax-b, for our computed solution x. Because of roundoff errors this residual may not be exactly zero, but all its components should be small; say less than 1×10-15.
x = A^(-1)*b % Solve the equation A x = b
A*x , b % Check x is correct by making
% sure it satisfies the equations.
resid = A*x - b % Check the residual is zero,
% at least to within roundoff.
In fact it is far quicker, and usually more accurate, to solve equations using the backslash operator ( \ ) introduced in the first section; rather than calculating with the inverse A-1. The next lines use the backslash command to the solve the equations Ax = b, and check that it gives the same answers as those that were calculated using A-1.
x1 = A \ b % This is the fastest way to solve A x = b x - x1 % Here both answers are the same.
(a.) Solve the system of equations
|
(b.) Use rand to generate a 700 x 700 matrix, A, and a column vector b of length 700. Solve the system of equations Ax = b using the commands x = A \ b and x = inv( A ) * b , timing each command with a watch. Which command is the faster and why? (Hint: before working with these large matrices use the clear command to remove all variables from Matlab's memory. This frees up enough memory to squeeze in this largish problem. If you have problems, try a smaller system of equations that fits onto the computer you are using.)
(c.) Use the Matlab command A = hilb(10) to set up the 10×10 Hilbert matrix. Create a right hand side vector b by b = A(:,1). We are now going to solve the equations Ax = b. This is a common test problem in linear algebra. (What is the true solution x to the equations Ax = b?)
(i) Solve the equations Ax = b using the backslash command, and then check that the computed solution satisfies the equations by calculating the residual, Ax-b.
(ii) Now solve the equations by the command A^(-1)*b and calculate the residual for this second solution.
(iii) Which gives the better (i.e. smaller) residual?
(iv) Which method gives the smaller error? (Explanation: As b is the first column of A, you can work out the true value of x without computation. (Explain this briefly.) The error is the difference between Matlab's computed solution and the true solution. (The true solution, the one you work out by abstract thought, can be typed directly into Matlab.) The size of the error is the largest absolute value of the components of the error vector.)?
(d.) Extend part (c.) by solving a slightly different problem. Let the right hand side vector b be the first column of A, divided by 3; i.e.b = A(:,1)/3. (What is the solution to the equation Ax = b with this new right hand side?)
(i) What is the size of the error and the size of the residual when the \ is used to solve the system with the new right hand side?
(ii) Extend the computation further by letting the right hand side be the j th column of A divided by 3; i.e. b=A(:,j)/3 for j = 1,2,3,¼,10. Complete the following table.
| Column j | Size of the Residual | Size of the Error |
| 1 | . | . |
| 2 | . | . |
| : | : | : |
| 10 | . | . |
This exercise should show that the backslash operator \ is both faster and more accurate than the inverse. Clearly it is better to use the backslash and avoid the calculating inverses in numerical work!
Knowing about vectors, we can use the Matlab procedures dealing with polynomials. The polynomial amXm+am-1xm-1+... +a1X+a0 is represented by the Matlab vector [am,am-1,... a1,a0], a row vector of length m+1. In the following example, X2-3X-4 is represented by the vector of coefficients [1 -3 -4].
v = [ 1 -3 -4 ] % Represents X^2 - 3 X -4 z = roots( v ) z(1)^ 2 -3*z(1) -4 % Two ways to evaluate the poly. polyval(v,z(1)) % at the first root.
(a.) The polynomial X2-3X-4 has real roots, but X2-3X+4 has complex roots. Find these complex roots with roots.
(b.) Use help to find out how to use Matlab to calculate the derivative of a polynomial in Matlab. Test this Matlab command by calculating the derivative of X4+X3+X2+X+1.