Starting from:
$31.99

$25.99

Assignment 2: LU Decomposition Solution

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

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

 

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 .

 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 , P4    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,

 1   0   0 

e.g.,  ∗    1   0  where ∗ is a potentially nonzero element.

∗    ∗    1

 

•  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 ith step,  for each row k below the ith 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 ith column of L.

 

(b)  Reduce  the  kth  row  of U  using  the  pivot  multiplier,  i.e.,  Uk =

Uk − (Uk,i /Ui,i )Ui where Ui is the ith row.

 

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 correspond- ing 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:

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    LU  Decomposition with Partial Pivoting

(4  points)

 

 

Based on your my_lu, you will write numerically  stable  LU decomposition  with partial pivoting.  At the  ith step  of LU decomposition (ith   pivot  column),  you will find the row that  has the largest absolute value in the pivot column (say row j), and swap the ith and jth rows of U as usual.  Simultaneously, you will swap the partial  entries  of the ith and jth 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 ith  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 ex- change  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 ith 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.

 

 

 

 

 

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

ith 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])

  1       0    0


4   −2      2  


0   0   1

L =  0.5     1   0, U = 0    4       8 , and P = 0   1   0

−0.5   1   1

 

•  Early  termination:


0     0     −4


1   0   0

[L,U,P] = my_lup([4,-2,2;-2,5,3;4+5E-13,-2+2E-13,2+7E-13])

  1       0    0


4   −2    2


0   0   1

L = −0.5    1   0, U = 0    4     4, and P = 0   1   0

1       0   1


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

 

 

 

 

 

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.

 

1.5


1.5

 

1                                                                                                                                               1

0.5


0.5

 

0                                                                                                                                               0

-0.5


 

-0.5         0          0.5          1          1.5          2          2.5          3          3.5


-0.5


 

 

-0.5         0          0.5          1          1.5          2          2.5          3          3.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 : R8  → R8 , and there exists some 8 × 8 matrix  A such that  f = Ad.

 

To visualize this relationship,  we have provided a function draw(), which can ei- ther 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  choos- ing a vector  d ∈ R8   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

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 ∈  R8 , 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 ∈ R8 , 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                                 y           g

PAd = Pf       ⇐⇒            L zU}|d{ = zP}|f{ .

 

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

 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:

−23.5355    −3.5355    · · ·

A ≈ 


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

 0.0185 

 0.1432 

              

−0.0269



              
 
 0.1212 

d ≈    0.0149    .

              

              

 0.1950 

              

 0.0177 

              

0.2063

More products