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:

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 [50]:

```
A=[0. 1. 3. 4.;
0. 1. 1. 0.]
n=size(A,2) # number of nodes
```

Out[50]:

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 [51]:

```
B=[1 2 3;
2 3 4]
m=size(B,2) # number of nodes
```

Out[51]:

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 [52]:

```
A[2,[1,3,4]],A[1:2,[1,3]],A[1,B[1:2,1]],vec(A[1,B[1:2,1]])
```

Out[52]:

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 [53]:

```
A[1:2,2:3]=[1 2; 3 4]
A
```

Out[53]:

**Exercise 2** What do you think the following returns?

In [54]:

```
A=[0. 1. 3. 4.;
0. 1. 1. 0.]
A[1:2,[1,3]]=[6 7; 8 9]
A
```

Out[54]:

Rather than working with matrices directly, it is better to work with a special type to represent structures. This way we can incorporate special behaviour directly. We are going to implement routines for the following type:

In [55]:

```
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 [56]:

```
function Structure(nodes::Matrix{Float64})
bars=zeros(Int,2,0) # creates a 0 x 2 array
Structure(nodes,bars)
end
```

Out[56]:

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 [57]:

```
numnodes(S::Structure)=size(S.nodes,2)
numbars(S::Structure)=size(S.bars,2)
```

Out[57]:

**(b)** Complete the following function that adds a bar for node `k`

to node `j`

:

In [58]:

```
function addbar!(S::Structure,k::Integer,j::Integer)
oldbars=S.bars
m=numbars(S)
newbars=zeros(2,m+1)
newbars[:,1:m]=oldbars
newbars[:,m+1]=[k,j]
S.bars=newbars # replace S.bars with the newbars
S # return S
end
```

Out[58]:

**3(c)** Verify the following code returns `[1 3; 2 4]`

In [59]:

```
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.])
addbar!(S,1,2)
addbar!(S,3,4)
S.bars
```

Out[59]:

`Tuples`

are similar to vectors but written with the notation `(x,y,z)`

instead of `[x,y,z]`

. They allow the storage of *different types*. For example:

In [60]:

```
t=(1,2.0,"hi")
```

Out[60]:

On the surface, this is very similar to a `Vector{Any}`

:

In [61]:

```
v=[1,2.0,"hi"]
```

Out[61]:

The main difference is that a `Tuple`

knows the type of its arguments:

In [62]:

```
typeof(t)
```

Out[62]:

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 [63]:

```
function mytuplereturn(x)
(cos(x),x^2)
end
```

Out[63]:

This returns a tuple:

In [64]:

```
mytuplereturn(5)
```

Out[64]:

But we can also employ the convenient syntax to create two variables at once:

In [65]:

```
x,y=mytuplereturn(5)
x
```

Out[65]:

**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 [66]:

```
function nodes(S::Structure,k)
S.nodes[:,S.bars[1,k]],S.nodes[:,S.bars[2,k]]
end
```

Out[66]:

**4(b)** Verify the following code returns `([3.,1.],[4.,0.])`

In [68]:

```
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.])
addbar!(S,1,2)
addbar!(S,3,4)
nodes(S,2)
```

Out[68]:

We now want to plot the associated structure. The following code plots the nodes:

In [70]:

```
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 [75]:

```
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 bars
for k=1:numbars(S)
a,b=nodes(S,k)
plot([a[1],b[1]],[a[2],b[2]])
end
# the following code changes the axis to
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
```

Out[75]:

**5(b)** Check that this plots correctly:

In [76]:

```
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.])
addbar!(S,1,2)
addbar!(S,3,4)
addbar!(S,1,3)
plot(S)
```

Out[76]:

The force each bar will make on the nodes will be a function of how long they are. Remember that the 2-norm gives the Euclidean length of a vector.

**Exercise 6(a)** Write a function called `barlength(S::Structure,k::Integer)`

that returns the length of the `k`

th bar.

In [77]:

```
function barlength(S::Structure,k::Integer)
a,b=nodes(S,k)
norm(a-b)
end
```

Out[77]:

**6(b)** Verify the following returns `3.1622776601683795`

:

In [78]:

```
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)
```

Out[78]:

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.

In [120]:

```
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
```

Out[120]:

In [121]:

```
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)
```

Out[121]:

**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 [134]:

```
function elongation(S::Structure,ref::Structure)
if numbars(ref)!=numbars(S)
error("The two structures have a different number of bars, so cannot be compared.")
end
m=numbars(ref)
ret=zeros(m)
for k=1:m
ret[k]=barlength(S,k)-barlength(ref,k)
end
ret
end
```

Out[134]:

**7(b)** Verify the following returns

```
3-element Array{Float64,1}:
0.0857864
0.366236
0.0950218
```

In [136]:

```
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)
```

Out[136]: