CSCI 421

Numerical Computing, Spring 2025

Homework Assignment 4

1.  (Taken from AG Ex 4.9, reworded for clarity).  Let

Carry out all the steps of GEPP (Gauss Elimination with Partial Piv- oting) on this matrix, writing out the matrices M(i), P(i) , i = 1, 2, 3, giving the upper triangular matrix

Also determine the matrices M(˜)(i) , i = 1, 2, 3, which should, like the M(i), be unit lower triangular (see AG, p. 121).  Then compute

and L =  M(˜)—1, which should also be unit lower triangular.   Finally check that, as desired,

PA = LU.

You can use a mix of hand calculations and matlab calculations, as you find convenient.  And you can double check that you got the right final answer with the matlab command  [L,U,P]=lu(A).

For the remainder of the homework, see next page

A Truss Problem

A typical task in structural engineering is to design a bridge to be strong  enough  to  withstand  a  certain  load.   Consider  the  following  plane  truss,  which is a set of metal bars connected by frictionless pin joints.  (“Plane” refers to the fact that the truss is two-dimensional, not three-dimensional  as it would be in reality.)  The symbol at the left end of the truss indicates  that it is fixed at that end, while the symbol at the right end indicates that  the truss is free to move horizontally, but not vertically.  The three arrows  pointing down represent loads on the truss.

The  problem  is to  solve  a  linear  system  of equations  for the  internal forces in the bars.  A positive internal force indicates that the bar is being extended (pulled apart a little), by the load, while a negative internal force indicates that the bar is being  compressed.  It is assumed that, as long as the internal forces are not too big, bars will not be stretched or compressed more than a tiny amount:  thus the structure does not collapse, but remains in  equilibrium.   By  computing  the  internal  forces,  an  engineer  has  more information as to whether the truss is indeed strong enough to withstand the load.

There are two linear equations for each internal joint in the truss, repre- senting forces in the horizontal and vertical direction which must balance at the joints.  Let us denote the internal forces by 儿1 , 儿2 , . . . , 儿13, correspond- ing to the numbers on the bars in the illustration.  The balancing of forces at joint C in the horizontal direction gives the equation

x4 = x8

while the balancing of forces at joint C in the vertical direction gives simply

x7  = 0.

The balancing of forces at joint B in the horizontal direction gives

x2 = x6

while the vertical direction at joint B gives

x3  = 10.

The  “10”  comes  from the  10 ton vertical  load  at joint  B.  The balancing of forces at joint A is a little more complicated, since it involves two bars oriented at an angle of 45 degrees as well as a horizontal and a vertical bar. Let  Q =  cos(π/4)  =  sin(π/4)  = √2/2. Then the  balancing of horizontal forces at joint A gives the equation

αx1  = x4 + αx5

and the balancing of vertical forces at joint A gives

αx1 + x3 + αx5  = 0.

There are also horizontal and vertical force equations at joints D, E and F which can be derived using the same ideas.  These amount to 12 equations altogether.  The 13th equation comes from the right end point G: since this end point is free to move horizontally, but not vertically, there is just one force equation, balancing the forces horizontally:

x13 + αx12  = 0.

Thus, we have a total of 13 linear equations defining the 13 internal force variables.

2.  Derive and solve the 13 linear equations in 13 variables.  Write the equations  using  matrix  notation,  as  Ax  =  b,  and  enter  the  matrix A and right-hand  side vector  b in  a  matlab function.   You  should start by allocating the space for A like this:  A=zeros(13,13).   Or- der the variables corresponding to bars 1, 2, 3, . . . and the equations corresponding to nodes A, B, C, . . . .  Then solve the system of linear equations using the matlab backslash operator:  x = A\b.  You  can put this in the function too, and return x, the vector of internal forces, as an output parameter of the function.  Print the solution vector x (in any format).  Mark the bars on a copy of the truss picture given above indicating which internal force is an extension force  (xj   > 0) and which is a compression force  (xj  < 0). Leave the bar unmarked if the internal force is zero. It’s fine to do this by hand, although you can write a program to do this if you want.

3.  Generalization:  Write  a new function that sets up and solves the equations  for  a  variable-sized  truss,   with  k  sections   exactly  like  the section ABCDEF instead of one, so that two neighboring sections share a vertical bar, with the EF bar from one section identical to the AB bar for the next section, where k is an input parameter to the function. Note that this means 4 new nodes and 8 new bars for each increase of k by 1.  This will require some careful thought.  Start by sketching the larger truss on paper and carefully writing down the relevant equations systematically.  Use orderings for the variables (i.e., the bars) and the equations  (i.e.,  the  nodes) that  are  consistent  with  the  case  k  =  1. Include plenty of comments explaining the code.  Don’t forget to start by allocating the space for A by A=zeros(n,n), where n depends on k. Make sure you recover the previous solution when k = 1, and then go on to test larger values of k, with a load vector that increases like this from left to right:  10, 20, 30, 40, 50, . . . .  Print the computed internal force vector x for k = 10.  There is no need to plot the truss.

4.  Sparsity:  The  bandwidth  of  A  (see  AG,  section  5.1)  is  defined  as p+q - 1 where p and q are the smallest nonnegative integers such that

ai,j  = 0    if    i ≥ j + p     or     j ≥ i + q.

You can visualize this using spy to display the nonzero entries in the matrix A, say, for k = 10.   What  is  the  bandwidth  of this  matrix  A? Would it be significantly di   erent if you  chose  di   erent  orderings,  say numbering  all the  top horizontal  bars first,  then  the  vertical  and diag- onal  bars,  and  then  the  bottom  horizontal  bars?  You  can  answer this either by  a clear  explanation,  or  changing the order  and displaying the results with spy. Also give the bandwidth for A—1 (the inverse of A, computed by inv(A)) and for the L and U factors obtained from [L,U] = lu(A). Remember that when you request only two output arguments  from  lu,  the  first,  L,  is  not  actually  lower  triangular  if pivoting (row interchanges) takes place; it is permuted lower triangu- lar  (or “psychologically” lower triangular).  Did pivoting occur?  How sparse are A—1 , L and U, compared with the sparsity of A?  Answer this carefully, comparing the four diferent spy plots.  You don’t need to use  matlab’s  sparse for this question,  although you  can if you want (the spy figures should be the same).

5.  Timing  Comparsions:  Execution  of  matlab code  can  be  timed using either timeit or tic . . .toc (see here for details).  Experiment with how long it takes to solve the system of equations for k that is large  enough  that  timing  comparisons   are  meaningful.   Compare  the following:

• x=A\b

• getting x by first computing A—1 and multiplying it onto b

• getting x by first computing the L and U factors with lu and then solving the two triangular systems Ly = b and Ux = y using \. (This is what is actually going on when you type x=A\b, so the timing should be about the same.)

• the same 3 again, but with A set up as a sparse matrix instead of the default “full” or dense matrix type; type help sparse for in- formation as to how this works. In particular, to initialize thema- trix to the sparse zero matrix, you need to write A=sparse(n,n) (if you use A=sparse(1,1) it should work but it will be slower because the data structure will change every time you assign a new matrix entry; do not use A=sparse(1) which will generate a sparse matrix with a11  = 1).  And if you don’t initialize A at all, it will be full, not sparse, by default.   This means editing your function that sets up A accordingly; add an input parameter that determines whether A is to be set up as a full (dense) or sparse matrix.

How do the execution times compare?  Does this relate to the sparsity displayed in the spy plots?

Finally compare the timings for larger values of k.   Do this just for solving  with  x=A\b,  but  for  two  cases:   A is  full  (dense),  and  A is sparse. Avoid making k too large as then the full case may run out of memory. Plot both running times on one plot as a function of k using semilogy, legend, title, xlabel, ylabel. Then make a second plot of the running times for the sparse case only, which should allow you to make k much larger – but don’t make it so large that you have to wait more than a few minutes for the program to run.

Submit:   listings  of the  matlab functions,  the  marked  truss  picture  re- quested in #2, the output x for the original truss and the generalized truss with k = 10, the plots generated by spy for k = 10, the derivation of the equations for both the original and the generalized truss, the results of the timing comparisons, including the plots of running times as a function of k , and answers to all the questions above.