# User:Eml4500.f08.ramrod.D/HW4

## October 10th, 2008

Note: Consider case where $\tilde{d}_{4}^{(e)} \ne 0$

$\tilde{d}_{1}^{(e)}=\tilde{d}_{2}^{(e)}=\tilde{d}_{3}^{(e)}= 0$

$\tilde{\mathbf{f}}_{4x1}^{(e)}=\tilde{\mathbf{k}}_{4x4}^{(e)}\tilde{\mathbf{d}}_{4x1}^{(e)}=\mathbf{0}_{4x1}$

The $\mathbf{0}_{4x1}$ is the 4th column of $\mathbf{k}^{(e)}$

Interpretation of the the transverse degrees of freedom

$\tilde{\mathbf{d}}^{(e)}=\tilde{\mathbf{T}}^{(e)}{d}^{(e)}$

Similarly: $\tilde{\mathbf{f}}^{(e)}=\tilde{\mathbf{T}}^{(e)}{f}^{(e)}$

Homework: Derive $\tilde{\mathbf{f}}^{(e)}=\tilde{\mathbf{T}}^{(e)}{f}^{(e)}$

Using the tranformation matrix we can show when $f^{(e)}$ is multiplied by $\tilde{\mathbf{T}}^{(e)}$ we will get $\tilde{\mathbf{f}}^{(e)}$

$\begin{Bmatrix} \tilde{f}_{1}^{(e)}\\ \tilde{f}_{2}^{(e)}\\ \tilde{f}_{3}^{(e)}\\ \tilde{f}_{4}^{(e)}\end{Bmatrix} = \begin{bmatrix} \mathbf{R}_{2x2}^{(e)} & \mathbf{0}_{2x2}\\ \mathbf{0}_{2x2} & \mathbf{R}_{2x2}^{(e)} \end{bmatrix} \begin{Bmatrix} f_{1}^{(e)}\\ f_{2}^{(e)} \\f_{3}^{(e)} \\f_{4}^{(e)}\end{Bmatrix}$

Also:

$\tilde{\mathbf{k}}^{(e)}{\tilde{\mathbf{d}}^{(e)}}=\tilde{\mathbf{f}}^{(e)}$

Therefore:

$\tilde{\mathbf{k}}^{(e)}\tilde{\mathbf{T}}^{(e)}\mathbf{d}^{(e)}=\tilde{\mathbf{T}}^{(e)}\mathbf{f}^{(e)}$

If $\tilde{\mathbf{T}}^{(e)}$ is invertible, then:

$[\tilde{\mathbf{T}}^{(e)-1}{\tilde{\mathbf{k}}^{(e)}}\tilde{\mathbf{T}}^{(e)}]\mathbf{d}^{(e)}=\mathbf{f}^{(e)}$

$\tilde{\mathbf{T}}^{(e)}$ block diagonal matrix: $\begin{bmatrix} \mathbf{R}_{2x2}^{(e)} & \mathbf{0}_{2x2} \\ \mathbf{0}_{2x2} & \mathbf{R}_{2x2}^{(e)} \end{bmatrix}$

Consider a general block diagonal matrix: $\mathbf{A}=\begin{bmatrix} \mathbf{D}_{1}& \mathbf{0} \\ \mathbf{0} & \mathbf{D}_{s} \end{bmatrix}$

Question: What is $\mathbf{A}^{-1}$?

Here is a simple example

$\mathbf{B}=\begin{bmatrix} \mathbf{d}_{11}& \mathbf{0}& \mathbf{0} \\ \mathbf{0} & \mathbf{d}_{22} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{d}_{nn} \end{bmatrix}$

$\mathbf{B}=Diag[d_{11}, d_{22}, ..., d_{nn}]$ (This is the shortened notation for $\mathbf{B}$)

$\mathbf{B}^{-1}=Diag[1/d_{11}, 1/d_{22}, ..., 1/d_{nn}]$

Assumming $d_{ii}\ne 0$ for $i=1, ..., n$

For a block diagonal matrix A:

$\mathbf{A}=Diag[d_{1} ..., d_{s}]$

$\mathbf{A}^{-1}=Diag[d_{1}^{-1} ..., d_{s}^{-1}]$

$\tilde{\mathbf{T}}^{(e)-1}=Diag[\mathbf{T}^{(e)-1} ..., R^{(e)-1}]$

We know:

$\mathbf{R}^{(e)T}=\begin{bmatrix} l^{(e)}& -m^{(e)} \\ m^{(e)} & -l^{(e)} \end{bmatrix}$

Also

$\mathbf{R}_{2x2}^{(e)}\mathbf{R}_{2x2}^{(e)T}=\begin{bmatrix} 1& 0 \\ 0 & 1 \end{bmatrix}_{2x2}=\mathbf{I}_{2x2}$ (This is the Identity matrix)

In conclusion: $\mathbf{R}^{(e)-1}=\mathbf{R}^{(e)T}$

So $\tilde{\mathbf{T}}^{(e)-1}=(Diag[\mathbf{R}^{(e)} ,\mathbf{R}^{(e)}])^{T}$ where $Diag[\mathbf{R}^{(e)} ,\mathbf{R}^{(e)}]=\tilde{\mathbf{T}}^{(e)}$

This shows that: $\tilde{\mathbf{T}}^{(e)-1}=\tilde{\mathbf{T}}^{(e)T}$

Homework: Verify $[\tilde{\mathbf{T}}^{(e)-1}{\tilde{\mathbf{k}}^{(e)}}\tilde{\mathbf{T}}^{(e)}]= \mathbf{k}^{(e)}$

We know the definitions of:

$\tilde{\mathbf{T}}^{(e)T} = \begin{bmatrix} l^{(e)} & -m^{(e)} & 0 & 0\\ m^{(e)} & l^{(e)} & 0 & 0\\ 0 & 0 & l^{(e)}& -m^{(e)}\\0 & 0 & m^{(e)} & l^{(e)}\\\end{bmatrix}$

$\tilde{\mathbf{k}}^{(e)} = k^{(e)} \begin{bmatrix} 1 & 0 & -1 & 0\\ 0 & 0 & 0 & 0\\ -1 & 0 & 1 & 0\\0 & 0 & 0 & 0\\\end{bmatrix}$

$\mathbf{T}^{(e)} = \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0\\ -m^{(e)} & l^{(e)} & 0 & 0\\ 0 & 0 & l^{(e)}& -m^{(e)}\\0 & 0 & m^{(e)} & l^{(e)}\\\end{bmatrix}$

So:

$\tilde{\mathbf{k}}^{(e)} \tilde{\mathbf{T}}^{(e)T} \mathbf{T}^{(e)} = \begin{bmatrix} 1 & 0 & 0& 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1& 0\\0 & 0 & 0 & 1\\\end{bmatrix}\tilde{\mathbf{k}}^{(e)} = k^{(e)}$

Therefore:

$\begin{bmatrix} 1 & 0 & 0& 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1& 0\\0 & 0 & 0 & 1\\\end{bmatrix}\tilde{\mathbf{k}}^{(e)}= k^{(e)}= \mathbf{k}^{(e)}$

## MATLAB

### 5 Truss MATLAB code

The 5 bar truss system in example 1.4 in the book can be solved by the below MATLAB code. analysis and explanation of the MATLAB code is shown below. This code and the the functions were downloaded from the books accompanying website:

http://bcs.wiley.com/he-bcs/Books?action=index&itemId=0471648086&bcsId=2256

% Five bar truss example
e1=200*10^3; e2=70*10^3; a1=40*100; a2=30*100; a3=20*100;
P = 150*10^3;
nodes = 1000*[0, 0; 1.5, 3.5; 0, 5; 5, 5];
conn = [1,2; 2,4; 1,3; 3,4; 2,3];
lmm = [1,2,3,4; 3, 4, 7, 8; 1, 2, 5, 6; 5, 6, 7, 8; 3, 4, 5, 6];
K=zeros(8);
% Generate stiffness matrix for each element and assemble it.
for i=1:2
lm=lmm(i,:);
con=conn(i,:);
k=PlaneTrussElement(e1, a1, nodes(con,:));
K(lm, lm) = K(lm, lm) + k;
end
for i=3:4
lm=lmm(i,:);
con=conn(i,:);
k=PlaneTrussElement(e1, a2, nodes(con,:));
K(lm, lm) = K(lm, lm) + k;
end

lm=lmm(5,:); con=conn(5,:);
k=PlaneTrussElement(e2, a3, nodes(con,:));
K(lm, lm) = K(lm, lm) + k

R = zeros(8,1); R(4)=-P

% Nodal solution and reactions
[d, reactions] = NodalSoln(K, R, [1,2,7,8], zeros(4,1))
results=[];
for i=1:2
results = [results; PlaneTrussResults(e1, a1, ...
nodes(conn(i,:),:), d(lmm(i,:)))];
end
for i=3:4
results = [results; PlaneTrussResults(e1, a2, ...
nodes(conn(i,:),:), d(lmm(i,:)))];
end
format short g
results = [results; PlaneTrussResults(e2, a3, ...
nodes(conn(5,:),:), d(lmm(5,:)))]


### Functions

The functions PlaneTrussElement.m, NodalSoln.m, PlaneTrussResults.m were all called in the code for the 5 truss system. They are all necessary for the completion of the truss problem.

PlaneTrussElement.m

This function take the Young's Modulus (e), the Area of the element (A) and the coordinates of the element ends (coord) and generates the stiffness matrix for a plane truss element. The function first calculates the lengths of the the elements. Once those are calculated the stiffness matrix is calculated.

function k = PlaneTrussElement(e, A, coord)
% PlaneTrussElement(e, A, coord)
% Generates stiffness matrix for a plane truss element
% e = modulus of elasticity
% A = area of cross-section
% coord = coordinates at the element ends

x1=coord(1,1); y1=coord(1,2);
x2=coord(2,1); y2=coord(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L;
k = e*A/L*[ls^2, ls*ms,-ls^2,-ls*ms;
ls*ms, ms^2,-ls*ms,-ms^2;
-ls^2,-ls*ms,ls^2,ls*ms;
-ls*ms,-ms^2,ls*ms,ms^2];


NodalSoln.m

The function takes the global coefficient matrix (K), the global right hand side vector (R), the list of degrees of freedom with specified values, and the specified values to determine the displacements and reactions at each node. The dof of the the system is first found by using the length command which finds the longest dimension of R. df then finds the difference between the dofs that have known values (a value of zero) and the dof that were found in the previous line. The displacements and the reactions are then calculated.

function [d, rf] = NodalSoln(K, R, debc, ebcVals)
% [nd, rf] = NodalSoln(K, R, debc, ebcVals)
% Computes nodal solution and reactions
% K = global coefficient matrix
% R = global right hand side vector
% debc = list of degrees of freedom with specified values
% ebcVals = specified values
dof = length(R);
df = setdiff(1:dof, debc);
Kf = K(df, df);
Rf = R(df) - K(df, debc)*ebcVals;
dfVals = Kf\Rf;
d = zeros(dof,1);
d(debc) = ebcVals;
d(df) = dfVals;
rf = K(debc,:)*d - R(debc);


PlaneTrussResults.m

This function computes the plane truss element results. It takes in the modulus of elasticity (e), the area of the cross-section (A), the coordinates at the element ends (coord), and the displacements at the elemtent ends (disps) and calculates the axial strain (eps), stress(sigma) and force (force). The results are stored in the variable results and sent back to the main program.

function results = PlaneTrussResults(e, A, coord, disps)
% results = PlaneTrussResults(e, A, coord, disps)
% Compute plane truss element results
% e = modulus of elasticity
% A = Area of cross-section
% coord = coordinates at the element ends
% disps = displacements at element ends
% The output quantities are eps = axial strain
% sigma = axial stress and force = axial force.

x1=coord(1,1); y1=coord(1,2);
x2=coord(2,1); y2=coord(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L;
T=[ls,ms,0,0; 0,0,ls,ms];
d = T*disps;
eps= (d(2)-d(1))/L;
sigma = e.*eps;
force = sigma.*A;
results=[eps, sigma, force];


### Results from the 5 Truss MATLAB example

The results from the 5 Truss System MATLAB code is shown below. As stated in the PlaneTrussResults.m function the results columns go in the order of axial strain, axial stress, and last axial force.

>>
K =

1.0e+005 *

Columns 1 through 7

0.3260    0.7607   -0.3260   -0.7607         0         0         0
0.7607    2.9749   -0.7607   -1.7749         0   -1.2000         0
-0.3260   -0.7607    2.4309    1.1914   -0.3300    0.3300   -1.7749
-0.7607   -1.7749    1.1914    2.4309    0.3300   -0.3300   -0.7607
0         0   -0.3300    0.3300    1.5300   -0.3300   -1.2000
0   -1.2000    0.3300   -0.3300   -0.3300    1.5300         0
0         0   -1.7749   -0.7607   -1.2000         0    2.9749
0         0   -0.7607   -0.3260         0         0    0.7607

Column 8

0
0
-0.7607
-0.3260
0
0
0.7607
0.3260

R =

0
0
0
-150000
0
0
0
0

d =

0
0
0.5390
-0.9531
0.2647
-0.2647
0
0

reactions =

1.0e+005 *

0.5493
1.5993
-0.5493
-0.0993

results =

-0.0001743      -34.859 -1.3944e+005
-3.15e-005      -6.2999       -25200
-5.2941e-005      -10.588       -31764
-5.2941e-005      -10.588       -31764
0.00032087       22.461        44922

>>


### Code for the plot of the undeformed and deformed five truss system

Below is the code for the plot of the five bar system. Using the results that were obtained from the above MATLAB solution the deformed plot was created. The dotted lines are the undeformed system and the solid lines are the deformed system. The the displacements were multiplied by a factor of 500 so that they would appear on the plot.

%Plot Five beam truss system
%
clear;
close;

%model with 2-D beam elements
dof = 2; %dof per node: axial disp x, 2= disp y
%obtain the coordinatates of all nodes
%
n_node = 4;             %number of nodes
n_elem = 5;            %number of elements
total_dof = 2 * n_node; %total dof of system

position(:, 1) = [0; 0];
position(:, 2) = [1500; 3500];
position(:, 3) = [0; 5000];
position(:, 4) = [5000; 5000]
% print the node coord.
for i = 1 : n_node;
x(i) = position (1,i);
y(i) = position (2,i);
end

node_connect (1, 1) = 1;    %element 1
node_connect (2, 1) = 2;

node_connect (1, 2) = 2;    %element 2
node_connect (2, 2) = 4;

node_connect (1, 3) = 1;    %element 3
node_connect (2, 3) = 3;

node_connect (1, 4) = 3;    %element 4
node_connect (2, 4) = 4;

node_connect (1, 5) = 2;    %element 5
node_connect (2, 5) = 3;

% connect all beam elements by connectivity array
for i = 1 : n_elem;
node_1 = node_connect(1,i);
node_2 = node_connect(2,i);
xx = [x(node_1),x(node_2)];
yy = [y(node_1),y(node_2)];
axis([ -1000 6000 -1000 6000])
plot(xx,yy,'--')
hold on
end

text(x(node_connect(1,1)),y(node_connect(1,1))-300,'Global Node 1', 'HorizontalAlignment', 'center')
text(x(node_connect(1,2))+800,y(node_connect(1,2))-800,'Global Node 2', 'HorizontalAlignment', 'center')
text(x(node_connect(2,3)),y(node_connect(2,3))+500,'Global Node 3', 'HorizontalAlignment', 'center')
text(x(node_connect(2,4)),y(node_connect(2,4))+500,'Global Node 4', 'HorizontalAlignment', 'center')

text (x(node_connect(2,1))*1.1,y(node_connect(2,1))/2 ,'Element 1', 'HorizontalAlignment', 'center')
text (x(node_connect(2,1))*2+300, y(node_connect(2,1)) ,'Element 2', 'HorizontalAlignment', 'center')
text( x(node_connect(2,3))/2-500, y(node_connect(2,3))/2 ,'Element 3', 'HorizontalAlignment', 'center')
text( x(node_connect(2,4))/2, y(node_connect(2,4))*1.1 ,'Element 4', 'HorizontalAlignment', 'center')
text( x(node_connect(1,4))*1.5+1500, y(node_connect(2,3))*.85 ,'Element 5', 'HorizontalAlignment', 'center')
hold on

position_disp(:,1) = [0;0];
position_disp(:,2) = [1769.5; 3023.45];
position_disp(:,3) = [132; 4866.5];
position_disp(:,4) = [5000; 5000]

% print the node coord.
for i = 1 : n_node;0
x(i) = position_disp (1,i);
y(i) = position_disp (2,i);
end

node_connect_disp (1, 1) = 1;    %element 1
node_connect_disp (2, 1) = 2;

node_connect_disp (1, 2) = 2;    %element 2
node_connect_disp (2, 2) = 4;

node_connect_disp (1, 3) = 1;    %element 3
node_connect_disp (2, 3) = 3;

node_connect_disp (1, 4) = 3;    %element 4
node_connect_disp (2, 4) = 4;

node_connect_disp (1, 5) = 2;    %element 5
node_connect_disp (2, 5) = 3;
% connect all beam elements by connectivity array
for i = 1 : n_elem;
node_1 = node_connect_disp(1,i);
node_2 = node_connect_disp(2,i);
xx = [x(node_1),x(node_2)];
yy = [y(node_1),y(node_2)];
axis([ -1000 6000 -1000 6000])
plot(xx,yy,'-')
hold on
end

title ('Five Bar Truss System')
xlabel('x')
ylabel('y')


File:Five bar truss.JPG
Undeformed and deformed five-bar truss system
 I, the copyright holder of this work, hereby release my work worldwide into the public domain. Where this is not legally possible, I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law.

### Three Bar Truss System

A program was also written to solve the three bar truss system. The code follows the model code that was given to us in HW 2. The functions PlaneTrussElement.m, NodalSoln.m, and PlaneTrussResults.m were once again used in the code to obtain the results. Descriptions of these functions can be found in the section for the Five bar truss system.

% Three bar truss example<br>
clear all;
e = [2 4 3]; A = [3 1 2]; P = 30;
L=[5 5 10];
alpha = [pi/6 -pi/6 pi/4];

nodes = [-L(1)*cos(alpha(1)), -L(1)*sin(alpha(1));
0, 0;
L(2)*cos(alpha(2)), -L(2)*sin(alpha(2));
-L(3)*cos(alpha(3)), -L(3)*sin(alpha(3))];

dof=2*length(nodes);

conn=[1,2; 2,3; 2,4];
lmm = [1, 2, 3, 4; 3, 4, 5, 6; 3, 4, 7, 8];
elems=size(lmm,1);
K=zeros(dof); R = zeros(dof,1);
debc = [1, 2, 5, 6, 7, 8];
ebcVals = zeros(length(debc),1);

R = zeros(dof,1); R(4) = P;

% Assemble global stiffness matrix
K=zeros(dof);
for i=1:elems
lm=lmm(i,:);
con=conn(i,:);
k_local=e(i)*A(i)/L(i)*[1 -1; -1 1]
k=PlaneTrussElement(e(i), A(i), nodes(con,:))
K(lm, lm) = K(lm, lm) + k;
end
K
R
% Nodal solution and reactions
[d, reactions] = NodalSoln(K, R, debc, ebcVals)
results=[];
for i=1:elems
results = [results; PlaneTrussResults(e, A, nodes(conn(i,:),:), d(lmm(i,:)))];

end
format short g
results


Below are the results from the above code.

>>
k_local =

1.2         -1.2
-1.2          1.2

k =

0.9      0.51962         -0.9     -0.51962
0.51962          0.3     -0.51962         -0.3
-0.9     -0.51962          0.9      0.51962
-0.51962         -0.3      0.51962          0.3

k_local =

0.8         -0.8
-0.8          0.8

k =

0.6      0.34641         -0.6     -0.34641
0.34641          0.2     -0.34641         -0.2
-0.6     -0.34641          0.6      0.34641
-0.34641         -0.2      0.34641          0.2

k_local =

0.6         -0.6
-0.6          0.6

k =

0.3         -0.3         -0.3          0.3
-0.3          0.3          0.3         -0.3
-0.3          0.3          0.3         -0.3
0.3         -0.3         -0.3          0.3

K =

Columns 1 through 5

0.9      0.51962         -0.9     -0.51962            0
0.51962          0.3     -0.51962         -0.3            0
-0.9     -0.51962          1.8      0.56603         -0.6
-0.51962         -0.3      0.56603          0.8     -0.34641
0            0         -0.6     -0.34641          0.6
0            0     -0.34641         -0.2      0.34641
0            0         -0.3          0.3            0
0            0          0.3         -0.3            0

Columns 6 through 8

0            0            0
0            0            0
-0.34641         -0.3          0.3
-0.2          0.3         -0.3
0.34641            0            0
0.2            0            0
0          0.3         -0.3
0         -0.3          0.3

R =

0
0
0
30
0
0
0
0

d =

0
0
-15.167
48.231
0
0
0
0

reactions =

-11.412
-6.5885
-7.6077
-4.3923
19.019
-19.019

results =

Columns 1 through 5

2.1962       4.3923       8.7846       6.5885       13.177
-2.1962      -4.3923      -8.7846      -6.5885      -13.177
4.4829       8.9658       17.932       13.449       26.897

Columns 6 through 7

8.7846       13.177
-8.7846      -13.177
17.932       26.897

>>


### Plot of undeformed and deformed three bar truss system

Below is the code and the plot of the three bar truss system. The plot was made from the results that were found in the above code. The results demonstrate that only node two had a displacement of (-15.167, 48.231). The undeformed system is plotted with the dotted lines and the deformed system is plotted with solid lines. The deformed displacements have been scaled down by a tenth in order to enhance in the plot.

%Plot three beam truss system
%
clear;
close;

%model with 2-D beam elements
dof = 2; %dof per node: axial disp x, 2= disp y
%obtain the coordinatates of all nodes
%
n_node = 4;             %number of nodes
n_elem = 3;            %number of elements
total_dof = 2 * n_node; %total dof of system

position(:, 1) = [-4.33; -2.5];
position(:, 2) = [0; 0];
position(:, 3) = [4.33; -2.5];
position(:, 4) = [-7.07106; -7.07106]
% print the node coord.
for i = 1 : n_node;
x(i) = position (1,i);
y(i) = position (2,i);
end

node_connect (1, 1) = 1;    %element 1
node_connect (2, 1) = 2;

node_connect (1, 2) = 2;    %element 2
node_connect (2, 2) = 3;

node_connect (1, 3) = 2;    %element 3
node_connect (2, 3) = 4;

% connect all beam elements by connectivity array
for i = 1 : n_elem;
node_1 = node_connect(1,i);
node_2 = node_connect(2,i);
xx = [x(node_1),x(node_2)];
yy = [y(node_1),y(node_2)];
axis([ -10 6 -10 10])
plot(xx,yy,'--')
hold on
end

text(x(node_connect(1,1)),y(node_connect(1,1))-.5,'Global Node 1', 'HorizontalAlignment', 'center')
text(x(node_connect(2,2)),y(node_connect(2,2))-1,'Global Node 2', 'HorizontalAlignment', 'center')
text(x(node_connect(1,3)),y(node_connect(1,3))+.5,'Global Node 3', 'HorizontalAlignment', 'center')
text(x(node_connect(2,3)),y(node_connect(2,3))-1,'Global Node 4', 'HorizontalAlignment', 'center')

text (x(node_connect(1,1))/2,y(node_connect(1,1))/2 +1,'Element 1', 'HorizontalAlignment', 'center')
text (x(node_connect(2,2))/2 -1, y(node_connect(2,2))/2 -.5,'Element 2', 'HorizontalAlignment', 'center')
text( x(node_connect(2,3))/2, y(node_connect(2,3))*.75 ,'Element 3', 'HorizontalAlignment', 'center')

%Scaled by 1/10

position_disp(:,1) = [-4.33; -2.5];
position_disp(:,2) = [-1.5167; 4.8231];
position_disp(:,3) = [4.33; -2.5];
position_disp(:,4) = [-7.07106; -7.07106]

% print the node coord.
for i = 1 : n_node;0
x(i) = position_disp (1,i);
y(i) = position_disp (2,i);
end

node_connect_disp (1, 1) = 1;    %element 1
node_connect_disp (2, 1) = 2;

node_connect_disp (1, 2) = 2;    %element 2
node_connect_disp (2, 2) = 3;

node_connect_disp (1, 3) = 2;    %element 3
node_connect_disp (2, 3) = 4;

% connect all beam elements by connectivity array
for i = 1 : n_elem;
node_1 = node_connect_disp(1,i);
node_2 = node_connect_disp(2,i);
xx = [x(node_1),x(node_2)];
yy = [y(node_1),y(node_2)];
axis([ -10 6 -10 10])
plot(xx,yy,'-')
hold on
end

title ('Five Bar Truss System')
xlabel('x')
ylabel('y')

File:Three bar truss.JPG
Undeformed and deformed three-bar truss system
 I, the copyright holder of this work, hereby release my work worldwide into the public domain. Where this is not legally possible, I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law.