In this lecture we introduce linear algebra. We first consider the problem of calculating the equilibrium of a system of four springs and three balls, affixed to walls at 0 and 4, where 0 is at the top and 4 is at the bottom.
Here is a plot of our system, where the green circles represent the balls and the blue dotted lines the spring:
In [1]:
using PyPlot # load the plot command
n=3
plot([0.,0.],[0.,-4.];marker="^",linestyle="--") # plot triangles to represent the
# anchors at 0 and 4
p=1:n # The resting position of the balls is equi spaced, so at 1, 2, and 3
plot(zeros(n),-p;marker="o",linestyle="") # plot the three balls
Out[1]:
Denote the displacement of the balls by $u_1,u_2,u_3$ and the elongation of the springs by $e_1,e_2,e_3,e_4$. Then we have
$$u_1=e_1, u_2 - u_1 = e_2, u_3 - u_2 = e_e, -u_3 = e_4$$Or in matrix form
$$ A\mathbf{u} = \mathbf{e} \qquad\hbox{for}\qquad A = \begin{pmatrix} 1 \cr -1 & 1 \cr &-1 &1 \cr &&-1 \end{pmatrix}, \mathbf{u} = \begin{pmatrix}u_1\cr u_2\cr u_3 \end{pmatrix}, \mathbf{e} = \begin{pmatrix}e_1\cr e_2\cr e_3\cr e_4 \end{pmatrix} $$Denote the forces of the springs by $y_1,y_2,y_3,y_4$. As the spring becomes more elongated, the force becomes stronger, which we model by Hook's law: $ y_k = c_k e_k$ where $c_k$ are the stiffness of each spring. This also can be written in matrix form: $$ C {\mathbf e} = {\mathbf y} \qquad \hbox{for} \qquad C = \begin{pmatrix} c_1\cr&c_2\cr&&c_3\cr&&&c_4\end{pmatrix} $$
We finally have the external forces $f_1,f_2,f_3$ pulling down on the balls, these could be gravity or something else. The external force on each ball has to balance with the forces from the spring, giving us:
$$f_1+y_2 = y_1, f_2+y_3 =y_2, f_3 +y_4 = y_3$$Or in matrix form
$$A^\top \mathbf{y} = \mathbf{f} \qquad \hbox{for} \qquad \mathbf{f} = \begin{pmatrix}f_1\cr f_2\cr f_3\end{pmatrix}$$We thus want to find the displacement of each ball by solving: $$A \mathbf{u} = \mathbf{e}, C\mathbf{e} = \mathbf{y}, A^\top \mathbf{y} = \mathbf{f}$$ which is equivalent to $$A^\top CAu= f$$ In other words, $$K\mathbf{u} =\mathbf{f}$$ for $K=A^\top C A$.
We will use Julia to solve this equation and plot the results. To do this, we need to create Vectors and Matrices, do Matrix-Matrix multiplication, take transposes and solve linear sytems.
In [2]:
v=zeros(5)
v[2]=3.3423
v
Out[2]:
We can specify the type of the Vector by adding it as the first argument to zeros:
In [3]:
v=zeros(Int64,5)
v[2]=3
v
Out[3]:
Note: we can't assign a Float to an integer vector:
In [4]:
v=zeros(Int64,5)
v[2]=3.5
We can also create vectors with ones and rand:
In [5]:
ones(5)
Out[5]:
In [6]:
ones(Int64,5)
Out[6]:
In [7]:
rand(5)
Out[7]:
In [8]:
rand(Int,5)
Out[8]:
We have already seen another way to create vectors directly:
In [9]:
[1,2,3,4]
Out[9]:
When the elements are of different types, they are mapped to a type that can represent every entry. For example, here we input a list of one Float64 followed by three Ints, which are automatically converted to Float64s:
In [10]:
[1.0,2,3,4]
Out[10]:
In the event that the types cannot automatically be converted, it defaults to an Any vector. This is bad performancewise so should be avoided.
In [11]:
[1.0,1,"1"]
Out[11]:
We can also specify the type of the Vector explicitly by writing the desired type before the first bracket:
In [12]:
Float64[1,2,3,4]
Out[12]:
We can also create an array using brackets, a formula and a for command:
In [13]:
[k^2 for k=1:5]
Out[13]:
In [14]:
Float64[k^2 for k=1:5]
Out[14]:
In [15]:
zeros(5,5) # creates a 5x5 matrix of Float64 zeros
Out[15]:
We can also have matrices of different types:
In [16]:
zeros(Int,5,5)
Out[16]:
eye creates the identity matrix:
In [17]:
eye(5)
Out[17]:
We can also create matrices by hand. Here, spaces delimit the columns and semicolons delimit the rows:
In [18]:
[1 2; 3 4; 5 6]
Out[18]:
In [19]:
Float64[1 2; 3 4; 5 6]
Out[19]:
We can also create matrices using brackets, a formula, and a for command:
In [20]:
[k^2+j for k=1:5,j=1:5]
Out[20]:
Matrices are really Vectors in disguise. They are still stored in memory in a sequence of addresses. We can see the underlying vector using the vec command:
In [21]:
M=[1 2; 3 4; 5 6]
vec(M)
Out[21]:
The only difference between matrices and vectors from the computers perspective is that they have a size which changes the interpretation of whats stored in memory:
In [22]:
size(M)
Out[22]:
Matrices can be manipulated easily on a computer. We can easily take determinants:
In [23]:
M=rand(3,3)
det(M)
Out[23]:
Or multiply:
In [24]:
M*M
Out[24]:
In [25]:
[1 2; 3 4]*[4 5; 6 7]
Out[25]:
If you use .*, it does entrywise multiplication:
In [26]:
[1 2; 3 4].*[4 5; 6 7]
Out[26]:
Vectors are thought of as column vectors, and so * is not defined:
In [27]:
a=[1,2,3]
b=[4,5,6]
a*b
Whereas entry-wise multiplication works fine:
In [28]:
a.*b
Out[28]:
Transposing a Vector gives a row vector, which is represented by a 1 x n matrix:
In [29]:
a'
Out[29]:
Thus we can do dot products as follows:
In [30]:
a'*b
Out[30]:
This is a vector with one entry, because matrix-vector multiplication always returns a vector. If we use dot, we get the dot product as a scalar:
In [31]:
dot(a,b)
Out[31]:
One important note: a vector is not the same as an n x 1 matrix:
In [32]:
v=zeros(Int,3,1)
v[1:3,:]=1:3 # the : notation means all columns
v
Out[32]:
In [33]:
a
Out[33]:
In [34]:
v==a
Out[34]:
Finally, we can solve linear systems use \:
In [35]:
A=rand(5,5)
b=ones(5)
u=A\b
Out[35]:
In [36]:
A*u-b
Out[36]:
In [37]:
n=3 # the number of balls
A=zeros(n+1,n)
for k=1:n
A[k,k]=1
A[k+1,k]=-1
end
A
Out[37]:
We now create $C$, where we assume the stiffness is equal. We then can create $K$ and $\mathbf{f}$ and solve the system using \ to find $\mathbf{u}$:
In [38]:
C=eye(n+1)
K=A'*C*A
f=[0,1,0.]
u=K\f # solves for the displacement u
Out[38]:
The masses are now shifted by u. We can plot the new locations in red versus the original locations in green as follows:
In [39]:
plot([0.,0.],[0.,-4.];marker="^",linestyle="--")
p=1:n
newp=p+u
plot(zeros(n),-p;marker="o",linestyle="")
plot(zeros(n),-newp;marker="o",linestyle="")
Out[39]:
Let's make this example interactive!
In [40]:
using Interact # let's you use @manipulate
In [41]:
n=3 # the number of balls
# set up A
A=zeros(n+1,n)
for k=1:n
A[k,k]=1
A[k+1,k]=-1
end
# the original locations of the springs
p=1:n
# creates a figure to modify
fig=figure()
# @manipulate creates two sliders: one for F and one for c. F represents the force
# on the middle ball and c the stiffness of the springs
@manipulate for F=-2.:.01:2., c=1.:10.
withfig(fig) do
C=c*eye(n+1) # this is the stiffnest matrix
K=A'*C*A
f=[0,F,0.]
u=K\f
newp=p+u
plot([0.,0.],[0.,-4.];marker="^",linestyle="--")
plot(zeros(n),-newp;marker="o",linestyle="")
end
end
Out[41]:
In [ ]: