## 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