functions about matrix on MatLab

Ricardo

The code should be stand-alone. No credit will be given if the function does not comply with the expected input and output.

Late submission policy: 25% o↵ up to 24 hours late; 50% o↵ up to 48 hours late; No point for more than 48 hours late.

2 Elementary row operations (30 points)

As this may be your first experience with serious programming in Matlab, we will ease into it by first writing some simple functions that perform the elementary row operations on a matrix: interchange, scaling, and replacement.

In this exercise, complete the following files:

function B = interchange(A, i, j) Input: a rectangular matrix A and two integers i and j. Output: the matrix resulting from swapping rows i and j, i.e. performing the row operation Ri $ Rj .

function B = scaling(A, i, s) Input: a rectangular matrix A, an integer i, and a scalar s. Output: the matrix resulting from multiplying all entries in row i by s, i.e. per- forming the row operation Ri sRi.

function B = replacement(A, i, j, s) Input: a rectangular matrix A, two integers i and j, and a scalar s. Output: the matrix resulting from adding s times row j to row i, i.e. performing the row operation Ri Ri + sRj .

Implementation tips:

In each of these functions, you should check that the input indices i and j are in range, i.e. 1  i  m and 1  j  m, where m is the number of rows in the matrix (which may not be the same as the number of columns!). If any index is out of range, print an error using the built-in function disp, and return the original matrix. This could help you diagnose problems when you write your RREF function in the next part.

Testing:

Test your code on rectangular m⇥n matrices of various sizes, both when m < n and when m > n. You can generate some simple matrices for testing on using the functions eye(m,n), ones(m,n), and randi(10,m,n).

2

3 RREF (40 points)

Next, you will use these row operations to write a function that performs Gauss- Jordan elimination and compute the reduced row echelon form of any matrix. We will call the function my_rref, because the rref function already exists in Matlab.

Specification:

function R = my_rref(A) Input: a rectangular matrix A. Output: the reduced row echelon form of A.

For full credit, your function should handle the following:

• Partial pivoting: At each step, you should swap the current row with the one whose entry in the pivot column has the largest absolute value.

• Free variables: Due to numerical error, the entries in a column corre- sponding to a free variable may be extremely small but not precisely zero. Therefore, you should consider an entry to be zero if its absolute value is smaller than 10�12.

We suggest first implementing the algorithm without considering these two issues, then adding code to deal with them one at a time.

Implementation tips:

There are two di↵erent ways one can implement Gauss-Jordan elimination.

• In Section 1.2 under “The Row Reduction Algorithm”, the book describes it in two phases: first do Gaussian elimination (Steps 1-4), then perform row operations equivalent to back-substitution (Step 5).

• Gauss-Jordan elimination can also be done in a single phase: every time you find a pivot, perform scaling so the pivot entry becomes 1, then perform elimination on all the other rows, both above and below the pivot row.

You may use either approach in your implementation. Below, we provide pseudocode for the latter approach.

In either case, since we want to be able to handle free variables, the pivot entry won’t necessarily be on the di- agonal. Instead, you’ll need to keep track of both the pivot row, say k, and the pivot column, l, as you go along; see the illustration on the right.

3

Algorithm 1 RREF

1: initialize pivot row k = 1, pivot column l = 1 2: while 1  k  m and 1  l  n do 3: among rows k to m, find the row p with the biggest* entry in column l 4: Rk $ Rp 5: if akl is zero** then 6: l l + 1 7: else 8: Rk (1/akl)Rk 9: for i = 1, . . . , k � 1, k + 1, . . . ,m do

10: Ri Ri � (ail/akl)Rk 11: end for 12: k k + 1, l l + 1 13: end if 14: end while

*Here “biggest” means having the largest absolute value (abs in Matlab). **Consider a number to be zero if its absolute value is smaller than 10�12.

Test cases:

my_rref([2,3,4; 1,1,0]) should return the matrix [1,0,-4; 0,1,4], i.e.  2 3 4 1 1 0

� !

 1 0 �4 0 1 4

� .

Partial pivoting: my_rref([0,6; -3,0]) should return the matrix [1,0; 0,1]:  0 6 �3 0

� !

 1 0 0 1

� .

Numerical error: my_rref([3,0,0,1; 0,3,0,1; 0,0,3,1; 1,1,1,1]): 2

664

3 0 0 1 0 3 0 1 0 0 3 1 1 1 1 1

3

775!

2

664

1 0 0 1/3 0 1 0 1/3 0 0 1 1/3 0 0 0 0

3

775

Free variables: my_rref([4,4,4,4; 2,2,4,4; 2,2,4,2]): 2

4 4 4 4 4 2 2 4 4 2 2 4 2

3

5!

2

4 1 1 0 0 0 0 1 0 0 0 0 1

3

5

You can also obtain a random m ⇥ n matrix with entries in {�1, 0, 1} by calling A = randi([-1,1], m, n), then compare your result my_rref(A) with Matlab’s built-in rref(A).

4

Note: Solving linear systems (no points)

Once we have an RREF function, we can use it to solve linear systems Ax = b by computing the RREF of the augmented matrix [A | b]. For simplicity, we will assume that the system has a unique solution. In this case, the RREF will be of the form 2

64 1 x1

. . . ...

1 xn

3

75

and the solution vector is just its last column. This is easy to implement in Matlab (which is already provided to you as in solve.m):

function x = solve(A, b) augmented_matrix = [A b]; R = my_rref(augmented_matrix); x = R(:, end);

end

Test this function on some simple examples of linear systems which you know to have unique solutions.

4 GPS localization (30 points)

One real-life application of solving linear systems is GPS localization. For simplicity, let us work in a 2D world first. Suppose your cellphone receives signals from three GPS satellites at known positions A = (a1, a2), B = (b1, b2), and C = (c1, c2), and can measure its distances rA, rB, rC from all of them, as shown below.

1 2 4 5( , ) ( , )a a= =A ( , )x y

5r =A

1 2 5 2( , ) ( , )b b= = −B

5r =B

1 2 11 6( , ) ( , )c c= = −C

13r =C

This gives us three quadratic equations for our own position P = (x, y):

r2A = (a1 � x)2 + (a2 � y)2, (1a) r2B = (b1 � x)2 + (b2 � y)2, (1b) r2C = (c1 � x)2 + (c2 � y)2. (1c)

5

Subtracting equation (1a) from equation (1b) and equation (1b) from equa- tion (1c), then simplifying, gives us a 2⇥ 2 linear system in x and y:

r2B � r2A = 2(a1 � b1)x+ 2(a2 � b2)y � a21 � a22 + b21 + b22, (2a) r2C � r2B = 2(b1 � c1)x+ 2(b2 � c2)y � b21 � b22 + c21 + c22. (2b)

Since this is a system of linear equations, it can be expressed in matrix form,  ⇤ ⇤ ⇤ ⇤

�  x y

� =

 ⇤ ⇤

� , (3)

for some matrix and some vector on the right-hand side that you should be able to figure out. Your task in this part of the assignment is to implement this method in Matlab (in gps2D.m and gps3D.m). That is, given the points A, B, C and distances rA, rB, rC, you will have to:

1. construct the matrix and the right-hand side vector in (3) corresponding to the linear system (2a)-(2b),

2. pass this matrix and vector to the solve function (already provided in the zip file) to obtain the point P. Note: if you did not get your my_rref function to work, you can change it to the Matlab built-in rref in the solve function in solve.m.

Exactly the same approach also works in 3D, except we will now need our distances from four known points A = (a1, a2, a3), B = (b1, b2, b3), C = (c1, c2, c3), and D = (d1, d2, d3). Your second task is to work out the analogous equations to (2a)-(2b), and implement another program that works in 3D.

Specification:

function p = gps2d(a, b, c, ra, rb, rc) Input: The 2D coordinates of three points A, B, C, and their distances rA, rB, rC from an unknown point P. Output: The 2D coordinates of the point P.

function p = gps3d(a, b, c, d, ra, rb, rc, rd) Input: The 3D coordinates of four points A, B, C, D, and their distances rA, rB, rC, rD from an unknown point P. Output: The 3D coordinates of the point P.

Test cases:

gps2d([4;5], [5;-2], [-11;6], 5, 5, 13) should return the vector [1;1].

gps3d([6;-3;3], [-1;-6;5], [-5;4;7], [6;8;4], 5, 7, 9, 9) should re- turn [2;0;3].

In general, you can create a random 2D vector with entries between say 0 and 10 by using a = 10*rand(2,1). Do this four times for a, b, c, and p, and set ra = norm(a - p) and so on. (The norm function returns the length of a vector.) Then check whether gps2d(a, b, c, ra, rb, rc) gives you back p.

6