APPM 4650 — HOMEWORK # 4The fall of the Eiffel tower
This is a translation of the original Eiffeltornet written by Professor Olof Runborg at
Royal Institute of Technology, Stockholm, Sweden.
Aside from the course book I recommend the Linear Equations chapter in Numerical Computing
with MATLAB by Cleve Moler. https://www.mathworks.com/moler/chapters.html
In many applications large systems of equations with thousands of unknowns must be solved
and the use of efficient algorithms becomes important. As an example your assignment is to study
the deformation of a complicated truss: a model of the Eiffel tower. A truss consists of a set of
rods connected in a set of nodes. The deformation is the result of an external force and is governed
by Newton’s law of motion under the assumption that the deformation is small and that Hooke’s
law is applicable. By approximation of the equations a system of equations Ax = b is obtained. If
there are N nodes there will be 2N unknowns (the x and y displacements of the N nodes) thus,
A ∈ R2N×2N. The matrix A is usually referred to as the stiffness matrix. The right hand side b
consists of the external forces
b = parenleftbigFx1 ,Fy1 ,Fx2 ,Fy2 ,...,FxN,FyNparenrightbig, b∈ R2N,
where (Fxj ,Fyj ) is the force in the jth node. The
The solution x contains the (unknown) displacements,
x= (∆x1,∆y1,∆x2,∆y2,...,∆xN,∆yN), x∈ R2N.
That is, (∆xj,∆yj)T, is the displacement of the jth node when the truss is loaded with the forces
in the vector b.
From the course homepage (D2L) you can download the files eiffel1.mat, eiffel2.mat,
eiffel3.mat, eiffel4.mat, eiffel5.mat and eiffel6.mat. They contain the matrix A, the
nodes, xnod, ynod (nod is the Swedish spelling of node) and the bar indices bars for four different
trusses of increasing fidelity (N = 261, 399, 561, 1592, 2290, 4330).
1. Load one of the models into Matlab using load. Download the program trussplot.m
from the course homepage and run trussplot(xnod,ynod,bars) to plot the tower. Pick
one of the nodes and load it with a force of magnitude 1/5 in the horizontal direction,
i.e. start with b = 0 and set Fxj = 1/5 for some j. In Matlab this would look like:
j=58; b=zeros(2*N,1); b(j*2-1)=1/5. Solve the system Ax = b using backslash and
compute the new positions of the nodes using xloadj = xnodej + ∆xj for all j. In Matlab
xload = xnod + x(1:2:end); yload = ynod + x(2:2:end);
Plot the loaded tower on top of the unloaded tower. Mark the loaded node.
2. Matlabs backslash x=A\b solves the system using Gauss elimination. Investigate how the time
it takes to obtain a solution depends on the size of the system Ax = b by solving for all four
models. Use the command cputime to measure time (NOTE that you should remove the plot
statements when you time the computation!). In order to get a reliable estimate for the time
it takes to solve the system (this is most important for small systems) you should measure the
time it takes to solve it many times and then use the average time. Plot or tabulate the time
it takes to to compute the solution as a function of N in a loglog plot. What is expected
from theory? Does the theory match your experiments?
3. Next you will estimate how sensitive the tower is to horizontal loading at different nodes.
Start with the smallest model eiffel1.mat. For each of the nodes load with the force
from above and compute the displacements. Record the size of the displacements bardblxbardbl for the
different nodes and find the node that gives the largest and smallest displacement (help max).
Plot the original tower and the most displaced tower with trussplot. Repeat for the larger
problems and tabulate the time it takes for the different models. Use a for loop to make the
computations:
dxnorm =[];
for j=1:N
b = zeros(2*N,1); b(2*j-1) = 1;
x = A\b;
dxnorm = [dxnorm norm(x)];
end
[fmax jmax] = max(dxnorm);
[fmin jmin] = min(dxnorm);
trussplot(xnod,ynod,bars); hold on
plot(xnod(jmax),ynod(jmax),’o’,xnod(jmin),ynod(jmin),’*’);
As the above problem solves the same system but with many different right hand sides it is
more efficient to first factor the matrix A and perform. back and forward substitution. Opti-
mize your program by using lu(A) and perform. forward / backward substitution. Estimate
how much shorter this makes your computation. Can you analyze the larger models now?
Tabulate the time it takes.
4. When the matrix is sparse it is possible to use methods that are far more efficient than Gauss
elimination. Use the command spy to investigate the structure of the matrix A. Let Matlab
know that A is sparse A = sparse(A); and perform. the computations in (c) again. How
much time do you gain by using the sparse methods in Matlab? Study both the case with
and without LU-decomposition. Tabulate the time it takes.