I used this code for another truss and it worked and now I am getting an error

This is the error I keep getting:

Index in position 1 exceeds array bounds (must not exceed 4).

Error in CEE210_CP4_Franke_IDK (line 289)

SupportReaction(j,:))

Here is my code: (Any help would be greatly appreciated!)

%% **********************************************************************
% Start of Input Section
% Input data for Benchmark Truss........................

% Units used: 
	LUnit = 'ft';	% Length = ft or m
	FUnit = 'lb';	% Force = lb or kip or N or kN

%.. Input problem data
    nNodes = 5;             % Number of nodes (pin joints) in the truss
    nMembers = 9;          % Number of members
    nReactions = 4;         % Number of force vector reactions
    nLoads = 2;             % Number of externally applied loads
    nDim = 3;               % Dimensions 2=2D, 3=3D 
% *** Once your code is modified and tested, nDim must always equal 3 and
% the variable nDim cannot be used in your code.
    nDOF = nDim * nNodes;   % Number of degrees of freedom
    tol = 1.e-8;            % Remainder tolerance on norm of residual force
 
% Input coordinates of each node
    Node(1,:) = [0, 0, 12];		% A
    Node(2,:) = [5, -5, 0];		% B
    Node(3,:) = [5, 5, 0];		% C
    Node(4,:) = [-5, 5, 0];		% D
    Node(5,:) = [-5, -5, 0];	% E
  
	
% Input member connections
%.. First input is the start node; Second input is the end node
    Member(1,:) = [1,2];	% AB
    Member(2,:) = [1,3];	% AC
    Member(3,:) = [1,4];	% AD
    Member(4,:) = [1,5];	% AE
    Member(5,:) = [2,3];	% BC
    Member(6,:) = [3,4];	% CD
    Member(7,:) = [4,5];	% DE
    Member(8,:) = [5,2];	% DG
    Member(9,:) = [2,4];	% EB
   
    
% Input support reactions    
%.. First input is the node number 
%.. Second input is resistance in the e1 direction
%.. Third input is resistance in the e2 direction
%.. Fourth input is resistance in the e3 direction (zero for 2-D trusses)
%*** 1 for resistance in that direction, 0 for no resistance ***
    Reaction(1,:) = [2, 1,1,1];		% Pin at Node 1
    Reaction(2,:) = [3, 0,0,1];		% Roller at Node 5
    Reaction(3,:) = [4, 1,0,0];
    Reaction(4,:) = [5, 0,0,0];
% Input externallly applied loads
%.. First input is the node number, remaining inputs are the magnitude of
%.. the force in the e1, e2 and e3 directions, respectively.
    Load(1,:) = [1, 0, 0, -2400];		% 1000 lb load at Node 2
    Load(2,:) = [4, 0, -1200, 0];			% 1000 lb load at Node 4
	
% End of Input Section for Benchmark Truss......
% ************************************************************************


%%  Compute the unit vector (direction) for every member

PosVec = zeros(nMembers, 3);		%Initialize Position Vector matrix
UnVec = zeros(nMembers, 3);		%Initialize Unit Vector matrix
OppUnVec = zeros(nMembers, 3);	%Initialize opposite direction matrix
Length = zeros(nMembers, 3);  %Initialize Length matrix

for i = 1 : nMembers
    SN = Member(i,1); 
    EN = Member(i,2); 

    
     PosVec(i,:) = Node(EN,:)- Node(SN,:);
     Length(i,:) = norm(PosVec(i,:));
     UnitVec(i,:) = PosVec(i,:)/Length(i,:);
     negUnitVec(i,:) = - UnitVec(i,:);
	OppUnVec(i,:) = -PosVec(i,:) / norm(PosVec(i,:));
end

%% The coefficient matrix, C, contains the unit vectors (direction vectors)
% for each member force associated with each node.  The unit vectors are
% then separated into e1, e2, e3 components and placed in separate rows. 
% Therefore, each node has 2 rows (for 2-D) or 3 rows (for 3-D) associated
% with it; one for each of the e1, e2, e3 components.
% This corresponds with the Degrees of Freedom, nDOF, which is 2*nNodes
% for 2-D or 3*nNodes for 3-D.
% So the coefficient matrix has 2 or 3 rows for each node and one column
% for each member force.

C = zeros(3*nNodes, nMembers);	%Initialize coefficient matrix

% Loop through all nodes, create the vector force equation for that node
% and store it in the proper row of the coefficient matrix, C.

for i = 1 : nNodes
	for j = 1 : nMembers
		SN = Member(j,1); 
		EN = Member(j,2);
		if (SN == i)  % If member j starts at node i
			C(3*i-2, j) = UnVec(j,1);
			C(3*i-1, j) = UnVec(j,2);
            C(3*i, j) = UnVec(j,3);
		elseif (EN == i)  % If member j ends at node i
			C(3*i-2, j) = OppUnVec(j,1);
			C(3*i-1, j) = OppUnVec(j,2);
            C(3*i, j) = OppUnVec(j,3);
		end  % if SN
	end  % for nMembers
end  % for nNodes


%% The External Force matrix, F, contains the magnitude of all externally
% applied loads, (input as e1, e2, e3 components), stored in the
% proper rows to corespond with the node it is applied at.
% Therefore, the F matrix has 2 or 3 rows for each node and one column.

F = zeros(3*nNodes, 1);		%Initialize External Force matrix

% Loop through all externally applied loads and place each component of the
% load (e1, e2, e3) in the proper rows in the external load matrix, F.

for i = 1 : nLoads
    j = Load(i,1);		% Which node, j is the load i on
	F(3*j-2) = Load(i,2);
	F(3*j-1) = Load(i,3);
    F(3*j) = Load(i,4);
end  % for i,nLoads


%% The Restraint matrix, Crest, contains the unit vectors (directions)
% for each member force associated with those nodes that are restrained by
% external supports.
% The Crest matrix has one row for each reaction component and one column
% for each member.

%.... Loop through all nodes and determine if there is a reaction at that
%.... node and seperate the reaction node equations from the force load
%.... equations.

nReaction = 1;       %Used to count the number of reaction equations

%.... The nReactionEquation vector determines if a row should be put in the
%.... Crest matrix or the Cfree matrix
nReactionEquation = zeros(1, 3*nNodes);

for (j = 1 : nReactions)
	for (i = 1 : nNodes)
		if (Reaction(j,1) == i)
			if (Reaction(j,2) == 1)
				Crest(nReaction,:) = C(3*i-2,:);
				nReactionEquation(3*i-2) = 1;       
				% If nReactionEquation is 1 it should be in the Crest
				% matrix and not the Cfree matrix.
				nReaction = nReaction+1;
			end
			if (Reaction(j,3) == 1)
				Crest(nReaction,:) = C(3*i-1,:);
				nReactionEquation(3*i-1) = 1;
				% If nReactionEquation is 1 it should be in the Crest
				% matrix and not the Cfree matrix.
				nReaction = nReaction+1;
            end
            if (Reaction(j,4) == 1)
				Crest(nReaction,:) = C(3*i,:);
				nReactionEquation(3*i) = 1;
				% If nReactionEquation is 1 it should be in the Crest
				% matrix and not the Cfree matrix.
				nReaction = nReaction+1;  
            end
           
        end
        
	end
end

nForce = 1;             %Used to count the number of Cfree equations

for (i = 1 : 3*nNodes)
	if (nReactionEquation(i) == 1)
	else
		Cfree(nForce,:) = C(i,:);
		Ffree(nForce) = F(i);
		nForce = nForce + 1;
	end
end


%% Solve first for the member forces and then solve for the support
% reactions.
% All of the values in the Cfree and Ffree matrices are known so you can
% now solve for the T matrix, which is the tension force in each member.

T = -(Cfree)\(Ffree');      % finish this equation

% Since all values in the Crest and T matrices are now known, you can solve
% for the support reactions.

Reactions =  -Crest*T;      % finish this equation

% Calculate the residual to find the amount of error and ensure that
% equilibrium was satisfied.
Residual1 = (Cfree * T + Ffree');
Residual2 = (Crest * T + Reactions);
Res = norm(Residual1) + norm(Residual2);

% Create the support reaction matrix.
nReaction = 1;
for (i = 1 : nNodes)
	if (nReactionEquation(3*i-2) == 1)
		SupportReaction(i,1) = i;
		SupportReaction(i,2) = Reactions(nReaction);
		nReaction = nReaction+1;   
	end
	if (nReactionEquation(3*i-1) == 1)
		SupportReaction(i,1) = i;
		SupportReaction(i,3) = Reactions(nReaction);
		nReaction = nReaction+1;   
    end
    if (nReactionEquation(3*i) == 1)
		SupportReaction(i,1) = i;
		SupportReaction(i,4) = Reactions(nReaction);
		nReaction = nReaction+1; 
    end
end


%% Create output for command window 
fprintf('%s\n'     ,  '----------------------------------------')
fprintf('%s\n'     ,  '-------------    Truss    --------------')
fprintf('%s\n'     ,  '----------------------------------------')
fprintf('%s',' Node Coordinates (', LUnit, ')')
fprintf ('\n')
for i = 1 : nNodes
	fprintf('%s%4i%8.3f%8.3f%8.3f\n','   Node: ', i, Node(i,:)')
end % for i, nNodes
	fprintf ('\n')
	fprintf('%s',' Member Forces (', FUnit, ')')
	fprintf ('\n')
for i = 1 : nMembers
	fprintf('%s%4i%12.3f\n','   Member Force: ',i, T(i)) 
end % for i, nMembers
	fprintf ('\n')
for i = 1 : nReactions
	for j = 1 : nNodes
		if (Reaction(i,1) == j)
			fprintf('%s%s%s%4i%12.3f%12.3f%12.3f\n',...
				'   Support Reaction (',FUnit,') at Node ',...
				SupportReaction(j,:)) 
			fprintf ('\n')
		end % if Reaction
	end % for j, nNodes
end % for i, nReactions
fprintf('\n%s%8.3f\n','   Residual Error: ', Res)
fprintf('\n\n\n')
      
      
%% Create the plot............................................
% Once your code is modified the plot will be 3-D and so you may need the
% "Rotate 3D" command to view the truss figure properly.

%.. Plot members with color indicating tension or compression 
Marker = ceil(25 / sqrt(nNodes));
fig1 = figure(1); clf; grid on; axis equal; 
hold on;
xlabel(cat(2, 'X  (',LUnit,')' )); 
ylabel(cat(2, 'Y  (',LUnit,')' ));
zlabel(cat(2, 'Z  (',LUnit,')' ));
title('Truss Analysis');
small = max(T)*tol;
for m = 1 : nMembers
	if (T(m) < -small)
		MColor = 'b';       % Color compression members blue
	elseif (T(m) > small)
		MColor = 'r';       % Color tension members red
	else
		MColor = 'k';       % Color "zero force members" black
	end % if T
	SN = Member(m,1); 
	EN = Member(m,2); 
	X = [Node(SN,1); Node(EN,1)];
	Y = [Node(SN,2); Node(EN,2)];
    Z = [Node(SN,3); Node(EN,3)];
	p = plot3(X,Y,Z);
	set(p, 'Color', MColor, 'LineWidth', 6/nDim);
end % for m, nMembers
		view(3);
%.. Plot external reaction forces
Rlength = 0.1 * max(max(Node));  % establish size of line
e = [ 1, 0, 0 ; 0, 1, 0 ; 0, 0, 1];
for i = 1 : nReactions
	for j = 1 : nNodes
		if (Reaction(i,1) == j)
			for k = 1 : 2
				if(Reaction(i, k+1) == 1)
					xR(1,:) = Node(j,:);
					xR(2,:) = Node(j,:) - Rlength * e(k, 1:nDim);
					xR(3,:) = Node(j,:) - Rlength * e(k, 1:nDim);
                    p = plot3(xR(:,1),xR(:,2),xR(:,3));
					set(p,'Color','g','LineWidth',8/nDim);
				end % if Reaction == 1
			end % for k
		end % if Reaction == j
	end % for j, nNodes
end % for i, nReactions

%.. Plot loads
Llength = 0.1 * max(max(Node));  % establish size of line
Fmax = max(max(abs(Load)));     % establish maximum force level
e = [ 1, 0, 0 ; 0, 1, 0 ; 0, 0, 1];      
for i = 1 : nLoads
	for j = 1 : nNodes
		if (Load(i,1) == j)
			for k = 1 : 2
				F = Load(i, k+1) / Fmax;
				xL(1,:) = Node(j,:);
				xL(2,:) = Node(j,:) + Llength * F * e(k, 1:nDim); 
				xL(3,:) = Node(j,:) + Llength * F * e(k, 1:nDim); 
                p = plot3(xL(:,1), xL(:,2), xL(:,3));
			set(p,'Color','c','LineWidth',8/nDim);
			end % for k
		end % if Load == j
	end % for j, nNodes
end % for i, nLoads
      
%.. Plot nodes      
plot3(Node(:,1), Node(:,2),Node(:,3), 'o', 'LineWidth', 2,...
	'MarkerFaceColor','w',...
	'MarkerEdgeColor','k',...
	'MarkerSize',Marker);
view(3);
		
% End of Program

Answer 1

[Walter Roberson]

for (i = 1 : nNodes)
	if (nReactionEquation(3*i-2) == 1)
		SupportReaction(i,1) = i;
		SupportReaction(i,2) = Reactions(nReaction);
		nReaction = nReaction+1;   
	end
	if (nReactionEquation(3*i-1) == 1)
		SupportReaction(i,1) = i;
		SupportReaction(i,3) = Reactions(nReaction);
		nReaction = nReaction+1;   
    end
    if (nReactionEquation(3*i) == 1)
		SupportReaction(i,1) = i;
		SupportReaction(i,4) = Reactions(nReaction);
		nReaction = nReaction+1; 
    end
end

That code only writes into SupportReaction(i,:) if at least one of those three conditions is true, and you are not pre-allocating SupportReaction. If there is an i for which none of the three tests is true, then SupportReaction(i,:) will not be written that time. However if there is a later i for which one of the tests is true, then writing into that later SupportReaction(i,:) would implicitly extend SupportReaction to have that many rows.

Therefore, SupportReaction will have as many rows as the last i for which one of the conditions hold. And that implies that SupportReaction will have fewer rows than nNodes if there are trailing i for which none of the conditions match.

K>> reshape(nReactionEquation,3,[])

ans =

     0     1     0     1     0
     0     1     0     0     0
     0     1     1     0     0

We can see from this at a glance that none of the conditions match for the final i in your data.

You could get around this problem by simply pre-allocating

   SupportReaction = zeros(3, nNodes);