% Cantilevered beam, three nodes, cant3.m
% W8x31 beam, 2 conc forces
clear                      % variables
format short e
f2 = -1e3                  % force at center

ex = 30e6                  % steel elastic modulus
ei = ex*110; s = 27.5      % E x I and section modulus

Len = 100
L2 = 2*ei/Len; L4 = 4*ei/Len
L6 = 6*ei/Len^2; L12 = 12*ei/Len^3
beta = L12 * 1e7

stiff1 = [ L12    L6   -L12     L6     % stiff matrix
          L6     L4   -L6      L2
         -L12   -L6    L12    -L6
          L6     L2   -L6      L4 ]

stiff = zeros(6); stiffa = zeros(6);
stiff(1:4,1:4)= stiff1;
stiffa(3:6,3:6)= stiff1;
stiff = stiff + stiffa
df =   [ f2; 0; f2; 0; 0*beta; 0*beta ] % combined d, F
zienk = stiff;
zienk(5,5)=zienk(5,5)+beta;
zienk(6,6)=zienk(6,6)+beta

format short e
displ = inv(zienk) * df;     % displacements
displt = displ'
force = (stiff * displ)'     % force vector

% print the results
fprintf('\n    Cantilevered beam, three nodes, cant3 ')
fprintf('\n      Your name, MENG 421, ')
disp(date)
fprintf('\n   Node   Deflec    Rotation     Force    Moment\n')
for i = 1:2:length(df-2)
  disp(sprintf('    %2.0f   %6.2f   %11.2e   %6.0f   %7.2e',...
    (i+1)/2,displ(i),displ(i+1),force(i),force(i+1)))
end

%     Cantilevered beam, three nodes, cant3 
%       Your name, MENG 421, 23-Apr-2004

%    Node   Deflec    Rotation     Force    Moment
%      1    -1.06     7.58e-003    -1000   -5.82e-011
%      2    -0.35     6.06e-003    -1000   -9.32e-012
%      3    -0.00     7.58e-007     2000   -3.00e+005