# Linear system for structures¶

We return to the Structure example from Lab 3. In this lab we will construct a linear system that represents the forces from elongating bars, to determine whether or not a structure is stable.

Recall that the nodes are 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)$.

A Structure S has two fields: nodes, which are a $2 \times m$ matrix

$$S.nodes=\begin{pmatrix} x_1 & x_2 & \cdots & x_n \cr y_1 & y_2 & \cdots & y_n \end{pmatrix}$$

and bars, a $2 \times m$ matrix

$$S.bars=\begin{pmatrix} p_1 & p_2 & \cdots & p_n \cr q_1 & q_2 & \cdots & q_n \end{pmatrix}$$

## Setup code¶

Here is the necessary code from the last lab, that sets up a type called Structure:

In [1]:
using PyPlot
import PyPlot.plot

type Structure
nodes::Matrix{Float64}
bars::Matrix{Int64}
end

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

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

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

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

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

INFO: Recompiling stale cache file /Users/solver/.julia/lib/v0.4/PyCall.ji for module PyCall.
INFO: Recompiling stale cache file /Users/solver/.julia/lib/v0.4/PyPlot.ji for module PyPlot.
INFO: Recompiling stale cache file /Users/solver/.julia/lib/v0.4/LaTeXStrings.ji for module LaTeXStrings.

Out[1]:
plot (generic function with 2 methods)

## Normals¶

The force given by a bar is dictated by how elongated the bar is. But now we are living in 2D: the force is at an angle specified by the angle of the bar. This is given by $\mathbf{n}$ defined in terms of two nodes $\mathbf c$ and $\mathbf d$:

$$\mathbf{n}(\mathbf{c},\mathbf{d}) := {\mathbf{c} - \mathbf{d} \over \|\mathbf c - \mathbf d\|}$$

Exercise 1 What is $\mathbf{n}(\mathbf{c},\mathbf{d})$ corresponding to the bar connecting the node $\mathbf{c} = (0,0)$ to $\mathbf{d} = (1/\sqrt 2,1/\sqrt 2)$?

Thus a bar connecting $\mathbf{c}$ and $\mathbf{d}$ puts a force in the direction of $\mathbf{n}(\mathbf{c},\mathbf{d})$ on node $\mathbf{b}$ and in the direction of $-\mathbf{n}(\mathbf{c},\mathbf{d})$ on node $\mathbf{a}$: it is pulling each node in the opposite direction, towards each other.

We will represent the normal directions coresponding to a Structure by a $2 \times m$ matrix:

$$N:=\left(\mathbf{n}(\mathbf{a}_{p_1},\mathbf{a}_{q_1}) \ |\ \mathbf{n}(\mathbf{a}_{p_2},\mathbf{a}_{q_2})\ |\ \cdots\ |\ \mathbf{n}(\mathbf{a}_{p_m},\mathbf{a}_{q_m}) \right) = \left(\mathbf{n}_1\ |\ \mathbf{n}_2\ |\ \cdots\ |\ \mathbf{n}_m \right) = \begin{pmatrix} n_1^x & n_2^x & \cdots & n_m^x \cr n_1^y & n_2^y & \cdots & n_m^y \end{pmatrix}$$

Exercise 2(a) Complete the following function that returns the normal matrix for a given structure:

In [2]:
normal(a,b)=(a-b)/norm(a-b)
function normals(S::Structure)
N=zeros(Float64,2,numbars(S))

for k=1:numbars(S)
a,b=nodes(S,k)
N[:,k]=normal(a,b)
end
N
end

Out[2]:
normals (generic function with 1 method)

2(b) What should the following return?

In [161]:
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.],
[1 2 3;
2 3 4])

normals(S)

Out[161]:
2x3 Array{Float64,2}:
-0.707107  -1.0  -0.707107
-0.707107   0.0   0.707107

## Plotting arrows¶

The arrow(x,y,δx,δy) command in PyPlot plots an arrow from the point (x,y) to (x+δx,y+δy), where the arrow head appears after the line. The optional arguments head_width and head_length can be used to set the size and shape of the arrow:

In [162]:
arrow(0.,0.5,1.,0.5; head_width=0.05,head_length=0.1)
axis([-1,2,0,2]);


Exercise 3(a) Complete the following function that, for each bar, plot two arrows depicting the direction of the force on each node. See the successful plot below as an example:

In [163]:
function plotnormals(S::Structure)
N=normals(S)

for k=1:numbars(S)
a,b=nodes(S,k)

end
end

Out[163]:
plotnormals (generic function with 1 method)

3(b) Ensure your plotnormals routine plots the following figure correctly:

In [164]:
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.],
[1 2 3;
2 3 4])

plot(S)
plotnormals(S)


## Forces from bars¶

$\mathbf n$ only indicate the direction in which the bars will put forces on the structure. The magnitude will be given by Hooke's law: if the bar $\mathbf{b}_k$ is elongated by an amount $e_k$, then the force will be of magnitude $y_k = c_k e_k$, where $c_k$ is the stiffness of the bar. We will take $c_k = 1$, and assume that a given Structure represents the bars at rest.

Therefore, elongations will be represented by a a vector of length $m$, one for each bar: $$\mathbf{e} = \begin{pmatrix} e_1 \cr e_2 \cr \vdots \cr e_m\end{pmatrix}$$

The direction of the force is given by $N$ above, that is, the force is $\mathbf y_k$ on node $\mathbf a_{q_k}$ and $-\mathbf y_k$ on node $\mathbf a_{p_k}$, where

$$\mathbf y_k = e_k \mathbf n_k.$$

Exercise 4(a) Complete the following function that plots arrows representing the forces of each bar: this time, the length of the arrow is given by $e_k$. (Hint: Copy-and-paste and modify plotnormals.)

In [90]:
function plotbarforces(S::Structure,e::Vector{Float64})
N=normals(S)

for k=1:numbars(S)
a,b=nodes(S,k)

end
end

Out[90]:
plotbarforces (generic function with 1 method)

4(b) Check that the following picture is plotted correctly, where the first and third bar are stretched while the second bar is squeezed:

In [102]:
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.],
[1 2 3;
2 3 4])

plot(S)
plotbarforces(S,[0.1,-0.2,0.3])


4(c) Can you write a one-line redefinition of plotnormals by calling plotbarforces?

In [73]:
plotnormals(S::Structure)=plotbarforces(S,ones(numbars(S)))

Out[73]:
plotnormals (generic function with 1 method)

## External forces¶

The external forces will be given by vectors $\mathbf{f}_1$, $\mathbf{f}_2$, …, $\mathbf{f}_n$, one external force per node. They will be represented by the $2 \times n$ matrix

$$F := \left(\mathbf f_1\ |\ \cdots\ |\ \mathbf f_n \right) = \begin{pmatrix} f_1^x & f_2^x & \cdots & f_n^x \cr f_1^y & f_2^y & \cdots & f_n^y \end{pmatrix}$$

Exercise 5(a) Complete the following function that plots the external forces given by the $2 \times n$ matrix F.

In [80]:
function plotexternalforces(S::Structure,F::Matrix)
for k=1:size(F,2)
a=S.nodes[:,k]

end
end

Out[80]:
plotexternalforces (generic function with 1 method)

5(b) Ensure that the following is plotted correctly:

In [101]:
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.],
[1 2 3;
2 3 4])

plotbarforces(S,[0.1,-0.2,0.3])
plotexternalforces(S,[.1 .2 -.3 .2;  -0.2 0.2 .4 -.1])
plot(S);


## Displacements¶

We want to model the displacement $\mathbf{u}_k$ of each node $\mathbf{a}_k$ under forces, or as a matrix

$$U =\left(\mathbf u_1\ | \ \cdots \ |\ \mathbf u_n \right)= \begin{pmatrix} u_1^x & \cdots & u_n^x \cr u_1^y & \cdots &u_n^y \end{pmatrix}$$

Now consider how far the bars are elongated. The true elongation is the length of the new bar, which is found by finding the difference between the two nodes:

$$e_k^{\rm true} = \|\mathbf{b}_{p_k} - \mathbf{b}_{q_k}\|$$

where $k=1,\ldots,m$. This leads to a nonlinear equation, unfortunately.

We assume that the displacements are small, and so we can replace the true elongation with a linear model. We will use

$$e_k = \mathbf{n}_k^\top (\mathbf{u}_{p_k} - \mathbf{u}_{q_k}) \qquad \hbox{for} \qquad \mathbf{n}_k := {\mathbf{a}_{p_k} - \mathbf{a}_{q_k} \over \|\mathbf a_{p_k} - \mathbf{a}_{q_k}\|}$$

We want to convert this equation giving $e_k$ in terms of $U$ into an equation. We first construct two $m \times n$ matrices, $M_x$ and $M_y$ so that

$$\mathbf e = M_x \begin{pmatrix} u_1^x \cr \vdots \cr u_n^x \end{pmatrix} + M_y \begin{pmatrix} u_1^x \cr \vdots \cr u_n^x \end{pmatrix}$$

Exercise 6(a) Why are these matrices given by the property that

$$M_x[k,p_k] = n_k^x$$$$M_x[k,q_k] = -n_k^x$$

$$M_y[k,p_k] = n_k^y$$ $$M_y[k,q_k] = -n_k^y$$

with zero entries otherwise, where $M_x[k,j]$ denotes the $(k,j)$ entry of $M_x$?

6(b) Complete the following routine that returns two matrices Mx and My that give $M_x$ and $M_y$ as defined above:

In [165]:
function incidentmatrices(S::Structure)
m=numbars(S)
n=numnodes(S)
B=S.bars

Mx=zeros(m,n)
My=zeros(m,n)
N=normals(S)
for k=1:m   # populate e_k
nx,ny=N[:,k]
Mx[k,B[1,k]]=nx
My[k,B[1,k]]=ny
Mx[k,B[2,k]]=-nx
My[k,B[2,k]]=-ny
end

Mx,My
end

Out[165]:
incidentmatrices (generic function with 1 method)

6(c) Ensure that Mx returns the correct value here:

In [166]:
S=Structure([0. 1. 3. 4.;
0. 1. 1. 0.],
[1 2 3;
2 3 4])

Mx,My=incidentmatrices(S)
Mx

Out[166]:
3x4 Array{Float64,2}:
-0.707107   0.707107   0.0       0.0
0.0       -1.0        1.0       0.0
0.0        0.0       -0.707107  0.707107

6(d) Without evaluating it, predict what My is. Are you correct?

In [167]:
My

Out[167]:
3x4 Array{Float64,2}:
-0.707107  0.707107   0.0        0.0
0.0       0.0       -0.0        0.0
0.0       0.0        0.707107  -0.707107

## Reducing a matrix of unknowns to a vector of unkowns¶

Are definition

$$\mathbf e = M_x \begin{pmatrix} u_1^x \cr \vdots \cr u_n^x \end{pmatrix} + M_y \begin{pmatrix} u_1^x \cr \vdots \cr u_n^x \end{pmatrix}$$

is not in the correct form for a linear system. We want to reduce it to a single vector of unknowns:

$$\mathbf u = \begin{pmatrix} u_1^x \cr u_1^y \cr u_2^x \cr u_2^y \cr \vdots \cr u_n^x \cr u_n^y \end{pmatrix}$$

We then combine $M_x$ and $M_y$ into a single $m \times 2n$ matrix:

$$\mathbf{e} = M \mathbf u$$

by interlacing:

$$M = \begin{pmatrix} M_x[1,1] & M_y[1,1] & M_x[1,2] & M_y[1,2] & \cdots & M_x[1,n] & M_y[1,n] \cr \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \cr M_x[m,1] & M_y[m,1] & M_x[m,2] & M_y[m,2] & \cdots & M_x[m,n] & M_y[m,n] \end{pmatrix}$$

Exercise 7(a) What Julia function reduces a Matrix $U$ to the Vector $\mathbf u$?

7(b) Complete the following function that combines Mx and My to form a single matrix M:

In [137]:
function incidentmatrix(S::Structure)
Mx,My=incidentmatrices(S)

M=zeros(size(Mx,1),2size(Mx,2))
for k=1:size(Mx,2)
M[:,2k-1]=Mx[:,k]
M[:,2k]=My[:,k]
end

M
end

Out[137]:
incidentmatrix (generic function with 1 method)

7(c) Check the following:

In [169]:
M=incidentmatrix(S)
U=[.1 .2 -.3 .2;  -0.2 0.2 .4 -.1]
M*vec(U)

Out[169]:
3-element Array{Float64,1}:
0.353553
-0.5
0.707107

7(d) Use plotbarforces along with M above to plot the forces given by $U$. What happens if $U$ is changed to $0.01U$?

In [170]:
plotbarforces(S,M*vec(U))
plot(S);


## Forces sum to zero¶

goal is to find displacements $\mathbf{u}_k$ so the forces on each node — the arrows in the above picture! — sum to zero: on node $k$, we want

$$0 = \mathbf{f}_k + \sum_{j=1}^m \begin{cases} -y_j \mathbf{n}_j & if p_j = k \cr y_j \mathbf{n}_j & if q_j = k \cr 0 & otherwise \end{cases}$$

Solving this will be the topic of the next lab. If you want to jump ahead, consider the following two questions:

Exercise 8(a) (Advanced) What is the equation relating $y_j$ and $vec(F)$, in terms of the matrix $M$?

8(b) (Advanced) Evaluate nullspace(M'*M). What issue does this raise about the problem we are trying to solve?

nullspace(M'*M)