$23.99
The Reduced Row Echelon Form Solution
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 manipu- lation 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. 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 RREF (4 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 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.
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
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 .
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 0 0 1/3
0 1 0 1/3
0 0 1 1/3
1 1 1 1
→
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
0 0
1
1
2 2
→ 0
0
1 0
1
1
2 1
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 ),
and C = (c1 , c2 ), and can measure its distances rA , rB , rC from all of them, as shown below.
C = (c1,c2 ) = (−11,6)
rC = 13
( x, y )
rA = 5
A = (a1,a2 ) = (4,5)
rB = 5
B = (b1,b2 ) = (5,−2)
This gives us three quadratic equations for our own position P = (x, y):
r2
A = (a1 − x)2
+ (a2 − y)2 , (1a)
r2
B = (b1 − x)2
r2 2
+ (b2 − y)2 , (1b)
2
C = (c1 − x)
+ (c2 − y) . (1c)
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:
r2 2
2 2 2 2
B − rA = 2(a1 − b1 )x + 2(a2 − b2 )y − a1 − a2 + b1 + b2 , (2a)
r2 2
2 2 2 2
C − rB = 2(b1 − c1 )x + 2(b2 − c2 )y − b1 − b2 + c1 + c2 . (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:
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.
