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.