MACM 316 – Computing Assignment 5
Submission Instructions: You must upload one .pdf file in Crowdmark that consists of two pages:
page 1 is your report which should fit all discussions, data and figures into a single page; and page 2 is a
listing of your code. The deadline is 11:00pm on the due date. The actual due time is set to 11:05pm
and if Crowdmark indicates that you submitted late, you will be assigned a grade of 0 on this assignment.
Your TA has emailed you a Crowdmark link that you should save since it will allow you to upload your
completed assignments.
❼ Please review the Guidelines for Assignments carefully.
❼ Acknowledge any collaborations or assistance from colleagues/TAs/instructor.
❼ If you have any questions about Matlab or aspects of this assignment, then you are strongly encouraged to attend tutorials and drop-in workshops.
Computing Assignment – Banded Matrices From Solid Mechanics
Consider a horizontal cantilevered beam of length L that is clamped at one end but free to move vertically
along the rest of its length. Divide the the beam into n equal sections of length h =
and separated
by points xi = ih, with i = 1, 2, . . . , n. A discrete model for the forces and solid deformations along the
Then the n-vector b = [bi
] is the given load force acting along the bar (which includes its own weight)
and z = [zi
] is the n-vector specifying the vertical deflection (away from z = 0) caused by the load force;
both bi and zi are measured at location xi
. We take the bar to be uniformly loaded, which means that
there is a constant upward force F distributed along the entire beam and so each component of the load
vector is a constant, bi = F · h
. For this problem, assume that the values of the physical constants are
L = 10 m and F = 0.005 N/m.
(a) Solve the linear system Az = b using the following three methods:
I. The GE+PP algorithm for sparse (banded) linear systems, which is the default algorithm used
by Matlab’s “\” operator when the matrix (call it Asparse) is of sparse type. You may find
it easiest to set up the matrix using the spdiags command.
II. The GE+PP algorithm for dense linear systems, again using “\”. Here, you need a dense
version of A which you can obtain either with the diag command or more simply by typing
Adense = full(Asparse)
III. The Gauss-Seidel iterative algorithm, also using the matrix Asparse. You may make use of the
gs2 code from lectures, setting the parameters tol = 1e-8 and maxiter = 1e5, and taking
initial guess z0 = (1, 1, . . . , 1)T
For each method, compute the solution with n = 100 and n = 1000 points. How do the three
methods compare in terms of cost? (use elapsed time from Matlab’s tic / toc as a proxy for cost).
How well do your three answers agree with each other? (use the vector 2-norm to compare the
differences). Which result do you think is most accurate?
(b) Next, verify that the matrix A has a UL factorization (yes, that’s UL and not LU) of the form
For the same two values of n as in the previous part (100, 1000) solve the linear system using the
UUT
factorization via a sequence of two triangular solves, corresponding to forward and backward
substitution. You can use the “\” command for both solves, but make sure your matrix U is set
up as a sparse matrix. Compare your answer to the corresponding solution from Method I (sparse
GE+PP) in terms of both accuracy and cost.
(c) Compute the condition number for both A and U, as well as the spectral radius ρ(T) for the GaussSeidel iteration matrix T. Use this information to help explain your results from parts (a) and (b).
The condition number for a large sparse matrix can only be estimated, and for this purpose you
should use the Matlab function condest.
(d) Choose what you consider to be your most accurate solution from parts (a)–(c), and use it to generate a plot of the deflection z versus x. In any introductory course in engineering solid mechanics,
students learn about deflection of elastic beams and derive the following formula for the maximum
deflection:
zmax = z(L) = F L46EI
where EI is a material property of the beam called flexural rigidity or bending stiffness (with units
N · m2). Use your solution to estimate the bending stiffness.
Extra reading: You can find out more about models for beam deflection at
https://en.wikipedia.org/wiki/Deflection_(engineering)