% Indetermin beam, three nodes, indet1.m
% W8x31 beam
clear                      % variables
format short e
f2 = -8e3                  % 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 =   [ 0*beta; 0*beta; f2; 0; 0*beta; 0 ]   % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
zienk(2,2)=zienk(2,2)+beta;
zienk(5,5)=zienk(5,5)+beta

displ = inv(zienk) * df;     % displacements
force = stiff * displ        % force vector

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

%  Indeterminate beam, three nodes, indet1
%   Your name, MENG 421, 20-Feb-2004

%  Node      Deflec     Rotation     Force    Moment
%     1     -0.0000    -7.58e-010    5500     3.00e+005
%     2     -0.1768    -7.58e-004   -8000     0.00e+000
%     3     -0.0000     3.03e-003    2500     0.00e+000