In this lab, we will investigate solving and plotting a 2D version of the spring system introduced in Lecture 6, where now the springs are replaced with bars that can be elongated but not bent, and the bars are attached at points to create structures. We will use linear algebra to understand which structures are stable: that is, won't fall down.
Remark You will not be examined on physics or modelling in this course. The physical systems studied are only to facilitate understanding the computational aspects of the problem.
We must first set up a way to interact. A structure consists of nodes and bars, where each bar is always attached to two nodes. We define the nodes by a sequence of $n$ points in 2D:
$$\mathbf{a}_1=(x_1,y_1),\mathbf{a}_2=(x_2,y_2),\ldots,\mathbf{a}_n=(x_n,y_n).$$The bars are then given by a sequence of $m$ pairs of nodes:
$$\mathbf{b}_1=(p_1,q_1),\mathbf{b}_2=(p_2,p_2),\ldots,\mathbf{a}_n=(p_n,q_n),$$where $p_j$ and $q_j$ are integers between 1 and $n$. For example, the bar attaching $\mathbf{a_2}$ to $\mathbf{a_4}$ is represented as $(1,4)$.
We will represent the nodes by a $2 \times m$ matrix
$$A=\begin{pmatrix} x_1 & x_2 & \cdots & x_n \cr y_1 & y_2 & \cdots & y_n \end{pmatrix}$$and the bars by a $2 \times m$ matrix
$$B=\begin{pmatrix} p_1 & p_2 & \cdots & p_n \cr q_1 & q_2 & \cdots & q_n \end{pmatrix}$$In Julia, the nodes $A$ will be represented by a Matrix{Float64}, for example, the following matrix represents nodes at $(0,0),(1,1),(3,1),(4,0)$:
In [ ]:
A=[0. 1. 3. 4.;
0. 1. 1. 0.]
n=size(A,2) # number of nodes
Note that size(A,1) gives the number of rows of A and size(A,2) gives the number of columns of A.
Theb bars $B$ will be represented by a Matrix{Int}. For example, the following represents the bars connecting nodes 1 and 2, 2 and 3, and 3 and 4:
In [ ]:
B=[1 2 3;
2 3 4]
m=size(B,2) # number of nodes
We will use the following notation to get at the columns and rows of matrices:
A[a:b,k] # returns the a-th through b-th rows of the k-th column of A as a Vector of length (b-a+1)
A[k,a:b] # returns the ath through bth columns of the k-th row of A as a 1 x (b-a+1) Matrix
A[:,k] # returns all rows of the k-th column of A as a Vector of length size(A,1)
A[k,:] # returns all columns of the k-th row of A as a 1 x size(A,2) Matrix
A[a:b,c:d] # returns the a-th through b-th rows and c-th through d-th columns of A
# as a (b-a+1) x (d-c+1) Matrix
The ranges a:b and c:d can be replaced by any AbstractVector{Int}. For example:
A[[1,3,4],2] # returns the 1st, 3rd and 4th rows of the 2nd column of A
The function vec turns a matrix into a vector:
vec(A[k,:]) # returns all columns of the k-th row of A as a size(A,2) Vector
Exercise 1 Can you guess what A[2,[1,3,4]] returns, using the definition of A as above? What about A[1:2,[1,3]]? And A[1,B[1:2,1]]? And vec(A[1,B[1:2,1]])?
In [ ]:
We can also use this notation to modify entries of the matrix. For example, we can set the 1:2 x 2:3 subblock of A to [1 2; 3 4] as follows:
In [ ]:
A[1:2,2:3]=[1 2; 3 4]
A
Exercise 2 What do you think the following returns?
In [ ]:
A=[0. 1. 3. 4.;
0. 1. 1. 0.]
A[1:2,[1,3]]=[6 7; 8 9]
A
In [ ]:
type Structure
nodes::Matrix{Float64}
bars::Matrix{Int64}
end
We want to be able to create a Structure just with nodes, and then add bars later. Last lecture we learned about constructors that verify arguments. But we can also create functions with the same name as a type: these are called external constructors. The following makes a structure with only nodes, and no bars:
In [ ]:
function Structure(nodes::Matrix{Float64})
bars=zeros(Int,2,0) # creates a 0 x 2 array
Structure(nodes,bars)
end
We want to be able to add bars one at a time. The convention in Julia is to end a function which modifies its input arguments with !.
Exercise 3(a) Add two functions, numnodes(S::Structure) and numbars(S::Structure), that return the number of nodes and bars, respectively.
In [ ]:
(b) Complete the following function that adds a bar for node k to node j:
In [ ]:
function addbar!(S::Structure,k::Integer,j::Integer)
oldbars=S.bars
m=numbars(S)
newbars=zeros(2,m+1)
## TODO: Populate newbars
S.bars=newbars # replace S.bars with the newbars
S # return S
end
3(c) Verify the following code returns [1 3; 2 4]
In [ ]:
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.])
addbar!(S,1,2)
addbar!(S,3,4)
S.bars
In [ ]:
t=(1,2.0,"hi")
On the surface, this is very similar to a Vector{Any}:
In [ ]:
v=[1,2.0,"hi"]
The main difference is that a Tuple knows the type of its arguments:
In [ ]:
typeof(t)
This is useful for the Julia compiler, but that is a topic beyond the scope of this course.
The main benefit of tuples for us is that they provide a convenient way to return multiple arguments from a function. For example, the following returns both cos(x) and x^2 from a single function:
In [ ]:
function mytuplereturn(x)
(cos(x),x^2)
end
This returns a tuple:
In [ ]:
mytuplereturn(5)
But we can also employ the convenient syntax to create two variables at once:
In [ ]:
x,y=mytuplereturn(5)
x
Exercise 4(a) Write a function nodes(S::Structure,k) that returns a tuple of the two nodes attached to bar k, each represented by a length 2 Vector.
In [ ]:
4(b) Verify the following code returns ([3.,1.],[4.,0.])
In [ ]:
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.])
addbar!(S,1,2)
addbar!(S,3,4)
nodes(S,2)
In [ ]:
using PyPlot
A=[0. 1. 3. 4.;
0. 1. 1. 0.]
x=vec(A[1,:]) # A vector of [x_1,x_2,…,x_n]
y=vec(A[2,:]) # a Vector of [y_1,y_2,…,y_n]
plot(x,y;marker="o",linestyle="") # plot dots at each coordinate specified by x and y
axis([-1,5,-1,1.5]) # change the axis so that x ranges between -1 and 5 and y ranges between 0 and 1.5;
Exercise 5(a) Complete the following function to plot the nodes as dots as above, and the bars as lines. Hint: It is easiest to call plot multiple times in a for loop, once for each row of B. Excercise 1 might help simplify the code. It is also easier to design functions line-by-line outside of a function definition.
In [ ]:
using PyPlot
import PyPlot.plot
function plot(S::Structure)
A=S.nodes
# plot the nodes
x=vec(A[1,:])
y=vec(A[2,:])
plot(x,y;marker="o",linestyle="")
## TODO: plot the bars
# the following code changes the axis so all nodes are visible
minx=minimum(A[1,:])
maxx=maximum(A[1,:])
miny=minimum(A[2,:])
maxy=maximum(A[2,:])
axis([minx-1,maxx+1,miny-1,maxy+1])
end
5(b) Check that this plots correctly:
In [80]:
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.])
addbar!(S,1,2)
addbar!(S,3,4)
addbar!(S,1,3)
plot(S);
In [ ]:
6(b) Verify the following returns 3.1622776601683795:
In [ ]:
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.])
addbar!(S,1,2)
addbar!(S,3,4)
addbar!(S,1,3)
barlength(S,3)
To calculate the forces on a structure, we need to compare it with a reference structure: this is the structure where the bars are at their natural length, and so no forces are acting on the structure. We have two quantities that we are interested in: the displacement of the nodes, and the elongation of the vectors.
The following code returns a 2 x n matrix containing the displacement of each of the n-nodes:
In [ ]:
function displacement(S::Structure,ref::Structure)
if numnodes(ref)!=numnodes(S)
error("The two structures have a different number of nodes, so cannot be compared.")
end
n=numnodes(ref)
ret=zeros(2,n)
for k=1:n
ret[:,k]=S.nodes[:,k]-ref.nodes[:,k]
end
ret
end
ref=Structure([0. 1. 3. 4.;
0. 1. 1. 0.])
S=Structure([0. 1.2 3.1 4.5;
0.1 1. 1.1 0.])
displacement(S,ref)
Exercise 7(a) Write a function elongation(S::Structure,ref::Structure) that returns a Vector of length numbars(S) containing how far the bars of S have been elongated when compared to the bars of ref.
In [ ]:
7(b) Verify the following returns
3-element Array{Float64,1}:
0.0857864
0.366236
0.0950218
In [ ]:
ref=Structure([0. 1. 3. 4.;
0. 1. 1. 0.])
addbar!(ref,1,2)
addbar!(ref,3,4)
addbar!(ref,1,3)
S=Structure([0. 1.2 3.1 4.5;
0.1 1. 1.1 0.])
addbar!(S,1,2)
addbar!(S,3,4)
addbar!(S,1,3)
elongation(S,ref)
Next lecture we will construct a linear system to determine whether a given structure is stable or not.