Description
In this assignment, you will implement a Matlab function to decompose a
matrix into lower and upper triangular matrices (L and U), i.e., PA = LU
where P is a row permutation matrix, and apply it to solve a computational
physics problem.
1 Download
For Section 6, we provide codes that can compute forces in a simple physical
system and visualize the results. This code has recently been updated.
Download the file sec6code.zip from the Moodle submission page, and extract
it into the folder where you will be working on the assignment. You should
obtain the files get_forces.m, draw.m, and a directory +defo.
2 Submission Guidelines
You will submit a zip file that contains the following .m files on Moodle:
• replacement_lu.m
• interchange_lup.m
• my_lu.m
• my_lup.m
• build_matrix.m
• get_displacements.m
Note 1: Each file must be named as specified above.
Note 2: Each function must return the specified value.
Note 3: You may use Matlab’s high-level array manipulation syntax, e.g. size(A),
A(i,:), and [A,B], but the built-in linear algebra functions such as inv, lu,
rref, and A\b are not allowed (except in Section 6). Please contact the TAs for
further questions.
Note 4: No collaboration is allowed. Plagiarism cases will be directly
reported to the University.
1
3 Elementary Row Operations (4 points)
You will implement two simple modular functions that will be useful for LU
decomposition. These are very similar to the functions you implemented in
Assignment 1.
Notation: for subsequent sections, Ai
indicates the i
th row of A, and Aij indicates
the (i, j) entry of A.
Specification:
function [U_new, L_new] = replacement_lu(U, i, j, s, L)
Input: two square matrices U and L, two integers i and j, and a scalar s.
Output:
• Unew: updated U by subtracting s times row j from row i, i.e. performing
the row replacement operation Ui ← Ui − sUj .
• Lnew: updated L by filling in its (i, j) entry with s, i.e., Lij ← s.
Algorithm 1 replacement lu
1: n ← the number of columns of U
2: Lij ← s
3: for 1 ≤ k ≤ n do
4: Uik ← Uik − sUjk . Row replacement
5: end for
6: return L, U
Warning! If you are adapting your code from Assignment 1, be careful that we
are now subtracting s times row j from row i instead of adding it.
2
function [U_new, L_new, P_new] = interchange_lup(U, i, j, L, P)
Input: three same size square matrices U, L, and P, and two integers i and j.
Output:
• Unew: updated U by swapping rows i and j, i.e. Ui ↔ Uj .
• Lnew: updated L by swapping only the first (i − 1) entries of rows i and j,
i.e., Li,1 ↔ Lj,1, . . . , Li,(i−1) ↔ Lj,(i−1) as shown in Figure 1.
• Pnew: updated P by swapping rows i and j, i.e. Pi ↔ Pj .
i
j
i
1 0 0 0
1
0
0 1 0
0
∗
∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
i
j
i
1 0 0 0
1
0
0 1 0
0
∗
∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
Partial row
exchange
n
n
i-1 i-1
Figure 1: Partial row exchange in L.
Algorithm 2 interchange lup
1: Ui ↔ Uj . Row interchange
2: Pi ↔ Pj . Row interchange
3: if i 1 then
4: for 1 ≤ k ≤ i − 1 do
5: Lik ↔ Ljk . Partial row interchange
6: end for
7: end if
8: return L, U, P
3
4 LU Decomposition (4 points)
In this part, you will write a function that performs LU decomposition without
pivoting. We will deal with pivoting in the next part of the assignment.
Specification:
function [L, U] = my_lu(A)
Input: an n × n square matrix A.
Output:
• L: an n × n lower triangular matrix where the diagonal entries are all one,
e.g.,
1 0 0
∗ 1 0
∗ ∗ 1
where ∗ is a potentially nonzero element.
• U: an n × n upper triangular matrix.
To get full credit, your function should handle the following case:
• Early termination: Due to round-off error, there may exist free variables
whose pivots are extremely small but not precisely zero. You should
terminate your LU decomposition if the absolute value of a pivot is less
than 10−12
.
The process of LU decomposition uses Gaussian elimination that transforms A
to an upper triangular matrix U while recording the pivot multipliers in a lower
triangular matrix L.
1. Initialize L to the identity matrix, and U to A. You can use Matlab’s
built-in function eye(n).
2. At the i
th step, for each row k below the i
th row,
(a) Record the pivot multiplier to L at (k, i), i.e., Lk,i = Uk,i/Ui,i as
shown in Figure 2. Note that this fills in the i
th column of L.
(b) Reduce the k
th row of U using the pivot multiplier, i.e., Uk =
Uk − (Uk,i/Ui,i)Ui where Ui
is the i
th row.
U = i
k
i
L = i
k
i
1 0 0 0
1
/ 1 0 Uki Uii
∗
∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
∗ ∗ ∗
0
0 0
U11
Uii
0 ∗ ∗
0 0
0 0
0
Figure 2: LU decomposition at the i
th step.
4
We provide pseudocode for LU decomposition in Algorithm 3.
Algorithm 3 LU decomposition of A
1: n ← the number of columns of A
2: L ← In where In is n × n identity matrix.
3: U ← A
4: for 1 ≤ i ≤ n − 1 do
5: if Uii is nearly zero* then
6: return L, U . Early termination
7: end if
8: for i + 1 ≤ k ≤ n do
9: p ← Uki/Uii
10: [U, L] = replacement_lu(U, k, i, p, L)
11: end for
12: end for
13: return L, U
*Consider a number to be nearly zero if its absolute value is smaller than 10−12
.
Note: When terminating early, L is not unique, i.e., the order of rows corresponding
to zero rows in U can be different. As long as the resulting decomposition
satisfies A = LU, the solution L and U are valid.
Test cases:
• [L,U] = my_lu([4,-2,2;-2,5,3;2,3,9])
L =
1 0 0
−0.5 1 0
0.5 1 1
, U =
4 −2 2
0 4 4
0 0 4
• Early termination:
[L,U] = my_lu([4,-2,2,1;-2,5,3,3;4,-2,2,1;-2,5,3,3])
L =
1 0 0 0
−0.5 1 0 0
1 0 1 0
−0.5 1 0 1
, U =
4 −2 2 1
0 4 4 3.5
0 0 0 0
0 0 0 0
.
Without early termination, it will return L and U with NaN (Not-a-Number)
entries.
• You can test your algorithm by constructing a random n×n matrix A, e.g.,
A=rand(10,10), and testing whether multiplying L and U reconstructs
A. You also need to check whether L is a lower triangular matrix with
diagonal entries equal to 1, and U is an upper triangular matrix.
5
5 LU Decomposition with Partial Pivoting
(4 points)
Based on your my_lu, you will write numerically stable LU decomposition with
partial pivoting. At the i
th step of LU decomposition (i
th pivot column), you
will find the row that has the largest absolute value in the pivot column (say row
j), and swap the i
th and j
th rows of U as usual. Simultaneously, you will swap
the partial entries of the i
th and j
th rows of L, and record the row exchange in
a permutation matrix P. For further details, please see
http://www.math.kent.edu/~reichel/courses/intr.num.comp.1/fall09/lecture9/lecture4.pdf
Specification:
function [L, U, P] = my_lup(A)
Input: an n × n square matrix A.
Output:
• L: an n × n lower triangular matrix where the diagonal entries are all 1.
• U: an n × n upper triangular matrix.
• P: an n × n permutation matrix.
The process of LU decomposition with partial pivoting needs to compute an
additional row permutation matrix P.
1. Initialize L and P to the identity matrix, and U to A. You can use
Matlab’s built-in function eye(n).
2. At the i
th step,
(a) Similar to Assignment 1, perform partial pivoting in U.
(b) Record the row exchange to the permutation matrix P by exchanging
its corresponding rows.
(c) Exchange the corresponding row entries in L to reflect the row exchange
in U. Note that this exchange only involves with the row
entries that have been recorded, i.e., not entire row exchange.
(d) For each row k below the i
th row, record the pivot multiplier to L and
replace the row in U using the pivot multiplier, like in the previous
part of this assignment.
6
We provide pseudocode for LUP decomposition in Algorithm 4.
Algorithm 4 LUP decomposition of A
1: n ← the number of columns of A
2: L ← In where In is n × n identity matrix.
3: P ← In . Permutation initialization
4: U ← A
5: for 1 ≤ i ≤ n − 1 do
6: j ← the row index of the largest absolute value among rows i to n in the
i
th column of U
7: [U, L, P] = interchange_lup(U, i, j, L, P)
8: if Uii is nearly zero* then
9: return L, U, and P . Early termination
10: end if
11: for i + 1 ≤ k ≤ n do
12: p ← Uki/Uii
13: [U, L] = replacement_lu(U, k, i, p, L)
14: end for
15: end for
16: return L, U, P
The red lines specify the major modification for partial pivoting.
*As before, consider a number to be nearly zero if its absolute value is smaller
than 10−12
.
Note: When terminating early, L and P are not unique, i.e., the order of rows
corresponding to zero rows in U can be different. As long as the resulting
decomposition satisfies PA = LU, the solution L and P are valid.
Test cases:
• Partial pivoting:
[L,U,P] = my_lup([-2,5,3;2,3,9;4,-2,2])
L =
1 0 0
0.5 1 0
−0.5 1 1
, U =
4 −2 2
0 4 8
0 0 −4
, and P =
0 0 1
0 1 0
1 0 0
• Early termination:
[L,U,P] = my_lup([4,-2,2;-2,5,3;4+5E-13,-2+2E-13,2+7E-13])
L =
1 0 0
−0.5 1 0
1 0 1
, U =
4 −2 2
0 4 4
0 0 0
, and P =
0 0 1
0 1 0
1 0 0
• You can test your algorithm by comparing its results with Matlab’s
built-in function, [l, u, p] = lu(A), on a randomly generated n × n
matrix A, e.g., A=rand(10,10).
7
6 Force-Displacement Computations (3 points)
Numerical linear algebra is used extensively in computational physics, with
applications in a wide range of fields including science, engineering, robotics, and
animation. Often, we represent a physical object as a collection of vertices, and
we have a computational model that can predict the internal forces that would
result from any given deformation of the object (i.e. any given displacement of
the vertices). Linear algebra allows us to do the reverse: apply any specified
external forces on the object, and predict how it will deform as a result.
-0.5 0 0.5 1 1.5 2 2.5 3 3.5
-0.5
0
0.5
1
1.5
-0.5 0 0.5 1 1.5 2 2.5 3 3.5
-0.5
0
0.5
1
1.5
Figure 3: Left: The rest shape of a simple system with four free vertices. Right:
A displacement of the vertices, and the resulting internal forces.
For example, consider the object above which has four free vertices (in black) that
can move around, and two constrained vertices (in red) that remain fixed. We
can collect the 2D displacements of the 4 free vertices into a single 8-dimensional
vector d. We have provided for you a function f = get_forces(d) that takes
this vector of displacements d as input and returns the resulting forces f on
the vertices, similarly packed into an 8-dimensional vector. In fact, the forces
are linear in the displacements, so the function is really a linear transformation
get forces : R
8 → R
8
, and there exists some 8 × 8 matrix A such that f = Ad.
To visualize this relationship, we have provided a function draw(), which can either
be called without arguments to draw the rest shape, or with two arguments
draw(d,f) to visualize a deformed shape with forces as arrows. Try choosing
a vector d ∈ R
8 with small entries and calling draw(d, get_forces(d))
to produce something like Fig. 3 (right). (Actually, you may have to draw
0.1*get_forces(d) because the forces are very large.)
8
Your task is to implement a function d = get_displacements(f, …) that
computes the inverse of this transformation: given the internal forces on the
free vertices, find the displacements that would give rise to them. To do this,
you will have to find the matrix A for the transformation, compute its LUP
factorization, and then use the factorization to quickly solve Ad = f to get d
for any given f.
Specification:
A = build_matrix()
Input: None.
Output: the 8 × 8 matrix A corresponding to the linear transformation
get_forces. For any vector d ∈ R
8
, the matrix-vector product A*d should
equal get_force(d).
Hint: Construct the vector e1 = [1, 0, . . . , 0]T and call get_forces on it. What
you get must be the first column of A. Do the same for the rest of the columns.
Once you have A, use either my_lup or Matlab’s lu to compute its LUP
factorization. You will need it for the following function:
d = get_displacements(f, L, U, P)
Input: a vector of forces f ∈ R
8
, and the LUP decomposition of A.
Output: the vector of vertex displacements d that result in these internal forces,
i.e. the solution of Ad = f.
Since PA = LU, we can write
PAd = Pf ⇐⇒ L
y
z}|{
Ud =
g
z}|{
Pf .
Thus we can compute d in three steps: compute the vector g = Pf, solve Ly = g
for y, and then solve Ud = y for d. In this function, you may use Matlab’s
backslash operator A\b to perform the two solves, since Matlab automatically
performs backsubstitution if the matrix is triangular.
Note: In the original code we provided, it turned out that the system
matrix you obtained from build_matrix did not require any pivoting, so
you did not actually need P to get the right answer. We have now provided
an updated version of the code which does induce pivoting. (It still models
the same system, but the ordering of the entries in f is different.) Download
the updated code and verify that your solution still works on it.
Once you have these functions working, you can compute the deformation caused
by any external forces fext: The resulting displacement is the one where the internal
and external forces cancel out, so d = get_displacements(-f_ext, …)
with a negative sign. For example, choose fext = [0, 0, 0, 0−0.2, −0.2, −0.2, −0.2]T
to apply a downward gravity force to all vertices, compute d, and visualize the
result using draw(d, f_ext).
9
Test cases:
These have been updated for the new version of the provided code.
For verification, we give the top left 2 × 2 block of the matrix A you should
obtain:
A ≈
−23.5355 −3.5355 · · ·
0 0 · · ·
.
.
.
.
.
.
.
.
.
Pulling the fourth vertex upwards with 1 unit of force, that is, f_ext = [0;0;0;0; 0;0;0;1],
should produce a displacement d = get_displacements(-f_ext, L,U,P) of
approximately
d ≈
0.0185
0.1432
−0.0269
0.1212
0.0149
0.1950
0.0177
0.2063
.
10
