## Description

In this assignment, you will implement a Matlab function to compute the

RREF of a matrix, and use it to solve a real-world problem that involves linear

algebra, namely GPS localization.

For each function that you are asked to implement, you will need to create a

.m file with the same name. For example, a function called my_func must be

defined in a file my_func.m. In the end, you will zip up all your .m files and

upload the zip file to the assignment submission page on Moodle.

In this and future assignments, you may not use any of Matlab’s built-in linear

algebra functionality like rref, inv, or the linear solve function A\b, except

where explicitly permitted. However, you may use the high-level array manipulation syntax like A(i,:) and [A,B]. See “Accessing Multiple Elements” and

“Concatenating Matrices” in the Matlab documentation for more information.

1 Elementary row operations (3 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.

Specification:

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. performing 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:

1

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 RREF (4 points)

Next, you will use these row operations to write a function that performs GaussJordan 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 corresponding 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 different 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.

2

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 diagonal. 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.

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([1,2,5; -2,1,0]) should return the matrix [1,0,1; 0,1,2], i.e.

1 2 5

−2 1 0

→

1 0 1

0 1 2

.

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

0 1

−1 0

→

1 0

0 1

.

3

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

3 0 0 1

0 3 0 1

0 0 3 1

1 1 1 1

→

1 0 0 1/3

0 1 0 1/3

0 0 1 1/3

0 0 0 0

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

2 2 2 2

1 1 2 2

1 1 2 1

→

1 1 0 0

0 0 1 0

0 0 0 1

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).

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

1 x1

.

.

.

.

.

.

1 xn

and the solution vector is just its last column. This is easy to implement in

Matlab:

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. If you did not get your my_rref function to work, you

can use the built-in rref instead.

3 GPS localization (3 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),

4

and C = (c1, c2), and can measure its distances rA, rB, rC from all of them, as

shown below.

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

r = 5 A

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

r = 5 B

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

r =13 C

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

r

2

A = (a1 − x)

2 + (a2 − y)

2

, (1a)

r

2

B = (b1 − x)

2 + (b2 − y)

2

, (1b)

r

2

C = (c1 − x)

2 + (c2 − y)

2

. (1c)

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

r

2

B − r

2

A = 2(a1 − b1)x + 2(a2 − b2)y − a

2

1 − a

2

2 + b

2

1 + b

2

2

, (2a)

r

2

C − r

2

B = 2(b1 − c1)x + 2(b2 − c2)y − b

2

1 − b

2

2 + c

2

1 + c

2

2

. (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. 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 to obtain the point P.

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:

5

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 return [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