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
