Code covered by the BSD License

# Uniform Beam Analysis

Find deflections and rotations for all nodes along the beam

Beam_Equation_Numerical_Substitution.m
```% This script solves the deflection and rotation for all nodes along the
% beam assuming that they are equally spaced along the beam's length
clear;clc
format short e
%%
% Inputs Parameters
n = input('How many element''s beam?  ');
L = input('Enter the length of beam element (in)  ');
E = input('Enter the modulus of elasticity (psi)  ');
I = input('Enter the mass moment of inertia (in^4)  ');
%%
% Stiffness Matrix Definition
k = E*I/L^3*[12   6*L    -12   6*L;
6*L  4*L^2  -6*L  2*L^2;
-12  -6*L     12  -6*L;
6*L  2*L^2  -6*L  4*L^2];
k1 = E*I/L^3*[12    6*L;
6*L   4*L^2];
k2 = E*I/L^3*[-12   6*L;
-6*L  2*L^2];
k3 = E*I/L^3*[-12  -6*L;
6*L   2*L^2];
k4 = E*I/L^3*[12   -6*L;
-6*L  4*L^2];

%%

index = ones(1,2*(n+1));                        % flag to catch the zero displacement or rotation
F_ext(1,2*(n+1)) = 0;
for i=1:n+1
fprintf ('what is the boundary condition at node  %2.0f ?\n',i);
fprintf ('Nothing = 0   Fixed = 1    Simply supported = 2    Free = 3\n');
BC(i) = input(' ');
if BC(i)==1
index([2*i-1 2*i]) = 0;                 % zero displacement and rotation if FIXED
elseif BC(i)==2
index(2*i-1) = 0;                       % zero displacement if SIMPLY SUPPORTED
end

F_ext(2*i-1) = input('Enter the external force at this node     ');
F_ext(2*i) = input('Enter the external moment at this node     ');
end
ind0 = find(index==0);              % catch the zero displacements and rotations to cancel
% the corresponding rows and columns
ind1 = find(index~=0);                  % catch the non-zero displacements and rotations

%%
for i = 2:n+1
j = i-1;
A([1 2],[1 2]) = 0;
A([2*i-1 2*i],[2*i-1 2*i]) = k4;
A([2*i-1 2*i],[2*j-1 2*j]) = k3;
end

for i = 1:n
j = i+1;
B([2*i-1 2*i],[2*i-1 2*i]) = k1;
B([2*i-1 2*i],[2*j-1 2*j]) = k2;
B([2*(n+1)-1 2*(n+1)],[2*(n+1)-1 2*(n+1)]) = 0;
end
K = A+B;
Stiffness = K;              % copy the stiffness matrix
F_ext = F_ext';             % external forces' column
Force_Mat = F_ext;          % copy the force matrix

Force_Mat(ind0)=[];           % delete the zero displacement and rotation columns in the external forces matrix
Stiffness(ind0,:) = [];       % delete the zero displacement and rotation rows in the stiffness matrix
Stiffness(:,ind0) = [];       % delete the zero displacement and rotation columns in the stiffness matrix
Displacement(1:2*(n+1)) = 0;    % define a displacement vector of zero's initially
D = Stiffness\Force_Mat;        % find the local nodal displacements
Displacement(ind1) = D;         % dump the nodal displacements into the overall displacement vector
% in the non-zero displacement locations

Global_Nodal = K * Displacement';   % find the global nodal forces
Force = Global_Nodal(1:2:end);      % the odd elements are the forces
Moment = Global_Nodal(2:2:end);     % the even elements are the elements
y_disp = Displacement(1:2:end);
rotation = Displacement(2:2:end);
Reaction_Forces = [1:n+1;Force';Moment';y_disp;rotation];
fprintf('    Node     Force            Moment          Y_disp             Rotation\n')
fprintf('    %2.0f    %+9.2f        %+9.2f        %+4.4e        %+4.4e   \n',Reaction_Forces)```