Error using symengine Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for

Error in Assifinal_Vruddhi (line 177)

KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA'; % KAA

My code is following

axes
clear
clc
close all
format short g
ndofn = 3; % number of degrees of freedom per node
% each number corresponds to a force for the corresponding DOF
Force = [0 0 0  0 0 0  0 0 0  0 -1000 0  0 -1000 0  0 -1000 0  0-1000 0  0 0 0  0 0 0  0 0 0]';
% Angles between local and global reference systems
% alphaA contains angles at the first node A
alphaA = [pi/2  pi/2  0  0  0  0  0  0    -pi/2];
% alphaB contains angles at the first node B
alphaB = [pi/2  0     0  0  0  0  0 -pi/2 -pi/2];
% distances along X and Y directions between origins of local reference systems
% in A and B; to be used in the tranlsation matrix
dX = [240 240 240 240 240 240 240 240 240]; %in
dY = [0  -240 0   0   0   0   0  -240 0]; %in
% Connectivity matrix: each column i shows the nodes that element i is connected to;
% For example, element 1 is connected to the structures nodes 1 and 2;
% In other words A==1  B==2 for the first element which is represented by
% the first column
t = [1 2 3 4 5 7 8 9;
    2 3 4 5 6 8 9 10];
%for straight
syms L x A E G Ay Iz real
%translation matrix for 2D element -> simplification of 3D case
HxB = [1 0 0;
    0 1 0;
    0 240 1];
HxB_T = HxB'; % transpose
% applies to linear portions of bridge
% fuq is cross sectional unit flexibility matrix
fuq_straight = [1/(A*E) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
fux_straight = HxB_T*fuq_straight*HxB;
fBA_straight = int(fux_straight, x, 0, L);
% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_straight = eval(inv(fBA_straight));
%for columns
syms L1 x A1 E G Ay1 Iz1 real
%translation matrix for 2D element -> simplification of 3D case
HxBc = [1 0 0;
    0 1 0;
    0 240 1];
HxB_Tc = HxBc'; % transpose
% applies to columns
% fuq is cross sectional unit flexibility matrix
fuq_column = [1/(A1*E) 0 0;
    0 1/(G*Ay1) 0;
    0 0 1/(E*Iz1)];
fux_column = HxB_Tc*fuq_column*HxBc;
fBA_column = int(fux_column, x, 0, L1);
% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_column = eval(inv(fBA_column));
%for curved
syms alpha theta r E G A Ay Iz real
c = cos(alpha);
s = sin(alpha);
% rotation operator
lambda_QB_T = [c s 0;
    -s c 0;
    0 0 1];
lambda_QB = lambda_QB_T';
% translation operator accounts also for the different (rotated) reference
% systems at different cross sections in the curved member
HQB =   [1 0 0;
    0 1 0;
    -240 240 1];
HQB_T = HQB';
% unitary flexibility matrix
fuq = [1/(E*A) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
TQB = lambda_QB_T*HQB; % Translation and rotation operator
TQB_T = TQB'; % Transpose
fux = TQB_T * fuq * TQB.*r;
fBA_curved = int(fux, alpha, 0, theta); %Flexibility matrix of curved member
theta = pi/2;
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
r = 240; %in
c1 = cos(theta);
s1 = sin(theta);
% Stiffness matrix of cantilever curved element
kB_curved = eval(inv(fBA_curved));
%global matrice
nnodes = length(t)+1;
% Initializing the global stiffness matrix of the entire structure
% Its dimension is square with nnodes*ndofn rows and columns
KKG = zeros(nnodes*ndofn,nnodes*ndofn);
% Initializing the global stiffness matrix of the one element
% Its dimension is square with 2*ndofn rows and columns
KG = zeros(2*ndofn,2*ndofn);
%Initialize Matrices
HAB = zeros(3,3);
lA = zeros(3,3);
lB = zeros(3,3);
for ii=1:length(t)
    HAB = [1 0 0;
        0 1 0;
        -dY(ii) dX(ii) 1];    % translation operator for straight members
    
    cA = cos(alphaA(ii));
    sA = sin(alphaA(ii));
    cB = cos(alphaB(ii));
    sB = sin(alphaB(ii));
    
    lA = [cA -sA 0;
        sA  cA 0;
        0   0  1];
    
    lB = [cB -sB 0;
        sB  cB 0;
        0   0  1];
    
    angle = pi/2;
    
    TAB = [  cos(angle)    sin(angle)    0;    % rotation/translation operator for curved members
        -sin(angle)    cos(angle)    0     % local ref systems in A and B are off by angle
        -r*(1-cos(angle)) r*sin(angle)  1];
    
    disp(['element= ',num2str(ii)])
    disp(['   '])
    
    if ii==1 || ii==6   % for straight components
        
        KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_straight*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_straight*lB';      % KAB
        KG(4:6,4:6) = lB*kB_straight*lB';      % KBB
        
    elseif ii==2 || ii==5     % for all columns
        
        KG(1:3,1:3) = lA*HAB*kB_column*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_column*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_column*lB';      % KAB
        KG(4:6,4:6) = lB*kB_column*lB';      % KBB
        
    else        % for curved members
        
        KG(1:3,1:3) = lA*TAB*kB_curved*TAB'*lA';    % KAA
        KG(4:6,1:3) = -lB*kB_curved*TAB'*lA';    % KBA
        KG(1:3,4:6) = -lA*TAB*kB_curved*lB';    % KAB
        KG(4:6,4:6) = lB*kB_curved*lB';    % KBB
        
    end
    
    node1 = t(1,ii);
    node2 = t(2,ii);
    
    AA = node1*3-2 : node1*3;  % AA contains the dofs for the node A of the member
    BB = node2*3-2 : node2*3;  % BB contains the dofs for the node B of the member
    
    % Adding the element stiffness matrix coefficients in the global
    % stiffness matrix of the entire structure
    
    KKG(AA,AA) = KKG(AA,AA) + KG(1:3,1:3);
    KKG(BB,AA) = KKG(BB,AA) + KG(4:6,1:3);
    KKG(AA,BB) = KKG(AA,BB) + KG(1:3,4:6);
    KKG(BB,BB) = KKG(BB,BB) + KG(4:6,4:6);
    
end

element= 1

Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Caused by:
Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.

% Eliminate rows and columns corresponding to constrained DOFs
Force([1:3,28:30],:) = [];
KK = KKG([1,3:6,10:12,16:19,21],[1,3:6,10:12,16:19,21]);
displ = KK\Force;
%Displacements of unconstrained nodes
ux = displ([4,7,10,13,16,19,22,25]);
uy = displ([11,14,17,20]);
phi = displ([6,9,12,15,18,21,24]);
% Table of nodal displacements for the structure
disp(['     Node       '   '  ux [in]       '   'uy [in]         '   'phi [rad]'])
disp([[2;4;6;7] , ux ,  [0;uy;0] ,    phi])
% Displacement components for all nodes
displALL_x = [ux(1),ux(2),0,ux(3),0,ux(4),ux(5)];
displALL_y = [0,uy(1),0,uy(2),0,uy(3),0];
n_coord = [000 120 120 195 270 270 390;...
    035 075 000 085 000 075 035].*12;
xcoord = n_coord(1,:);
ycoord = n_coord(2,:);
scale = 100;
figure(1), plot(xcoord,ycoord,'ks')
hold on
plot(xcoord + scale.*displALL_x , ycoord + scale.*displALL_y,'ro')
legend('undeformed nodes','deformed nodes')
xlabel('X axis [in]')
ylabel('Y axis [in]')
axis equal
grid on

Answer 1

[Walter Roberson]

axes
clear
clc
close all
format short g
ndofn = 3; % number of degrees of freedom per node
% each number corresponds to a force for the corresponding DOF
Force = [0 0 0  0 0 0  0 0 0  0 -1000 0  0 -1000 0  0 -1000 0  0-1000 0  0 0 0  0 0 0  0 0 0]';
% Angles between local and global reference systems
% alphaA contains angles at the first node A
alphaA = [pi/2  pi/2  0  0  0  0  0  0    -pi/2];
% alphaB contains angles at the first node B
alphaB = [pi/2  0     0  0  0  0  0 -pi/2 -pi/2];
% distances along X and Y directions between origins of local reference systems
% in A and B; to be used in the tranlsation matrix
dX = [240 240 240 240 240 240 240 240 240]; %in
dY = [0  -240 0   0   0   0   0  -240 0]; %in
% Connectivity matrix: each column i shows the nodes that element i is connected to;
% For example, element 1 is connected to the structures nodes 1 and 2;
% In other words A==1  B==2 for the first element which is represented by
% the first column
t = [1 2 3 4 5 7 8 9;
    2 3 4 5 6 8 9 10];
%for straight
syms L x A E G Ay Iz real
%translation matrix for 2D element -> simplification of 3D case
HxB = [1 0 0;
    0 1 0;
    0 240 1];
HxB_T = HxB'; % transpose
% applies to linear portions of bridge
% fuq is cross sectional unit flexibility matrix
fuq_straight = [1/(A*E) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
fux_straight = HxB_T*fuq_straight*HxB;
fBA_straight = int(fux_straight, x, 0, L);

fBA_straight is the integral to the unbound symbolic variable L, so fBA_straight is going to be symbolic.

% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_straight = eval(inv(fBA_straight))

kB_straight = 

You assigned to A1 and L1 and ay1 not to A and L and ay, so kB_straight remains symbolic after the eval() step.

%for columns
syms L1 x A1 E G Ay1 Iz1 real
%translation matrix for 2D element -> simplification of 3D case
HxBc = [1 0 0;
    0 1 0;
    0 240 1];
HxB_Tc = HxBc'; % transpose
% applies to columns
% fuq is cross sectional unit flexibility matrix
fuq_column = [1/(A1*E) 0 0;
    0 1/(G*Ay1) 0;
    0 0 1/(E*Iz1)];
fux_column = HxB_Tc*fuq_column*HxBc;
fBA_column = int(fux_column, x, 0, L1);
% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_column = eval(inv(fBA_column));
%for curved
syms alpha theta r E G A Ay Iz real
c = cos(alpha);
s = sin(alpha);
% rotation operator
lambda_QB_T = [c s 0;
    -s c 0;
    0 0 1];
lambda_QB = lambda_QB_T';
% translation operator accounts also for the different (rotated) reference
% systems at different cross sections in the curved member
HQB =   [1 0 0;
    0 1 0;
    -240 240 1];
HQB_T = HQB';
% unitary flexibility matrix
fuq = [1/(E*A) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
TQB = lambda_QB_T*HQB; % Translation and rotation operator
TQB_T = TQB'; % Transpose
fux = TQB_T * fuq * TQB.*r;
fBA_curved = int(fux, alpha, 0, theta); %Flexibility matrix of curved member
theta = pi/2;
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
r = 240; %in
c1 = cos(theta);
s1 = sin(theta);
% Stiffness matrix of cantilever curved element
kB_curved = eval(inv(fBA_curved));
%global matrice
nnodes = length(t)+1;
% Initializing the global stiffness matrix of the entire structure
% Its dimension is square with nnodes*ndofn rows and columns
KKG = zeros(nnodes*ndofn,nnodes*ndofn);
% Initializing the global stiffness matrix of the one element
% Its dimension is square with 2*ndofn rows and columns
KG = zeros(2*ndofn,2*ndofn);

You assign KG as numeric zeros.

%Initialize Matrices
HAB = zeros(3,3);
lA = zeros(3,3);
lB = zeros(3,3);
for ii=1:length(t)
    HAB = [1 0 0;
        0 1 0;
        -dY(ii) dX(ii) 1];    % translation operator for straight members
    
    cA = cos(alphaA(ii));
    sA = sin(alphaA(ii));
    cB = cos(alphaB(ii));
    sB = sin(alphaB(ii));
    
    lA = [cA -sA 0;
        sA  cA 0;
        0   0  1];
    
    lB = [cB -sB 0;
        sB  cB 0;
        0   0  1];
    
    angle = pi/2;
    
    TAB = [  cos(angle)    sin(angle)    0;    % rotation/translation operator for curved members
        -sin(angle)    cos(angle)    0     % local ref systems in A and B are off by angle
        -r*(1-cos(angle)) r*sin(angle)  1];
    
    disp(['element= ',num2str(ii)])
    disp(['   '])
    
    whos
    if ii==1 || ii==6   % for straight components
        
        KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA';      % KAA

You use the symbolic kB_straight on the right-hand side, so the right hand side comes out symbolic involving unbound variables. But the left-hand side is purely numeric, and numeric arrays cannot hold expressions involving unbound symbolic variables.

        KG(4:6,1:3) = -lB*kB_straight*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_straight*lB';      % KAB
        KG(4:6,4:6) = lB*kB_straight*lB';      % KBB
        
    elseif ii==2 || ii==5     % for all columns
        
        KG(1:3,1:3) = lA*HAB*kB_column*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_column*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_column*lB';      % KAB
        KG(4:6,4:6) = lB*kB_column*lB';      % KBB
        
    else        % for curved members
        
        KG(1:3,1:3) = lA*TAB*kB_curved*TAB'*lA';    % KAA
        KG(4:6,1:3) = -lB*kB_curved*TAB'*lA';    % KBA
        KG(1:3,4:6) = -lA*TAB*kB_curved*lB';    % KAB
        KG(4:6,4:6) = lB*kB_curved*lB';    % KBB
        
    end
    
    node1 = t(1,ii);
    node2 = t(2,ii);
    
    AA = node1*3-2 : node1*3;  % AA contains the dofs for the node A of the member
    BB = node2*3-2 : node2*3;  % BB contains the dofs for the node B of the member
    
    % Adding the element stiffness matrix coefficients in the global
    % stiffness matrix of the entire structure
    
    KKG(AA,AA) = KKG(AA,AA) + KG(1:3,1:3);
    KKG(BB,AA) = KKG(BB,AA) + KG(4:6,1:3);
    KKG(AA,BB) = KKG(AA,BB) + KG(1:3,4:6);
    KKG(BB,BB) = KKG(BB,BB) + KG(4:6,4:6);
    
end

element= 1

  Name               Size            Bytes  Class     Attributes

  A                  1x1                 8  double              
  A1                 1x1                 8  double              
  Ay                 1x1                 8  sym                 
  Ay1                1x1                 8  double              
  E                  1x1                 8  double              
  Force             29x1               232  double              
  G                  1x1                 8  double              
  HAB                3x3                72  double              
  HQB                3x3                72  double              
  HQB_T              3x3                72  double              
  HxB                3x3                72  double              
  HxB_T              3x3                72  double              
  HxB_Tc             3x3                72  double              
  HxBc               3x3                72  double              
  Iz                 1x1                 8  sym                 
  Iz1                1x1                 8  double              
  KG                 6x6               288  double              
  KKG               27x27             5832  double              
  L                  1x1                 8  sym                 
  L1                 1x1                 8  double              
  TAB                3x3                72  double              
  TQB                3x3                 8  sym                 
  TQB_T              3x3                 8  sym                 
  alpha              1x1                 8  sym                 
  alphaA             1x9                72  double              
  alphaB             1x9                72  double              
  angle              1x1                 8  double              
  c                  1x1                 8  sym                 
  c1                 1x1                 8  double              
  cA                 1x1                 8  double              
  cB                 1x1                 8  double              
  dX                 1x9                72  double              
  dY                 1x9                72  double              
  fBA_column         3x3                 8  sym                 
  fBA_curved         3x3                 8  sym                 
  fBA_straight       3x3                 8  sym                 
  fuq                3x3                 8  sym                 
  fuq_column         3x3                 8  sym                 
  fuq_straight       3x3                 8  sym                 
  fux                3x3                 8  sym                 
  fux_column         3x3                 8  sym                 
  fux_straight       3x3                 8  sym                 
  ii                 1x1                 8  double              
  kB_column          3x3                72  double              
  kB_curved          3x3                 8  sym                 
  kB_straight        3x3                 8  sym                 
  lA                 3x3                72  double              
  lB                 3x3                72  double              
  lambda_QB          3x3                 8  sym                 
  lambda_QB_T        3x3                 8  sym                 
  ndofn              1x1                 8  double              
  nnodes             1x1                 8  double              
  nu                 1x1                 8  double              
  r                  1x1                 8  double              
  s                  1x1                 8  sym                 
  s1                 1x1                 8  double              
  sA                 1x1                 8  double              
  sB                 1x1                 8  double              
  t                  2x8               128  double              
  theta              1x1                 8  double              
  x                  1x1                 8  sym                 

Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Caused by:
Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.

% Eliminate rows and columns corresponding to constrained DOFs
Force([1:3,28:30],:) = [];
KK = KKG([1,3:6,10:12,16:19,21],[1,3:6,10:12,16:19,21]);
displ = KK\Force;
%Displacements of unconstrained nodes
ux = displ([4,7,10,13,16,19,22,25]);
uy = displ([11,14,17,20]);
phi = displ([6,9,12,15,18,21,24]);
% Table of nodal displacements for the structure
disp(['     Node       '   '  ux [in]       '   'uy [in]         '   'phi [rad]'])
disp([[2;4;6;7] , ux ,  [0;uy;0] ,    phi])
% Displacement components for all nodes
displALL_x = [ux(1),ux(2),0,ux(3),0,ux(4),ux(5)];
displALL_y = [0,uy(1),0,uy(2),0,uy(3),0];
n_coord = [000 120 120 195 270 270 390;...
    035 075 000 085 000 075 035].*12;
xcoord = n_coord(1,:);
ycoord = n_coord(2,:);
scale = 100;
figure(1), plot(xcoord,ycoord,'ks')
hold on
plot(xcoord + scale.*displALL_x , ycoord + scale.*displALL_y,'ro')
legend('undeformed nodes','deformed nodes')
xlabel('X axis [in]')
ylabel('Y axis [in]')
axis equal
grid on

[Walter Roberson]

kB_straight = eval(inv(fBA_straight))

You should never eval() a symbolic expression. The meaning of eval() of a symbolic expression is not documented. In practice, it is the same as eval(char(inv(fBA_straight))) but char() of a symbolic expression returns something that is not compatible with MATLAB and is not compatible with the symbolic engine.

You should use subs() on symbolic expressions -- possibly followed by double()

[Vruddhi]

KG(4:6,1:3) = -lB*kB_straight*HAB'*lA';      % KBA
KG(1:3,4:6) = -lA*HAB*kB_straight*lB';      % KAB
KG(4:6,4:6) = lB*kB_straight*lB';      % KBB

elseif ii==2 || ii==5     % for all columns
    
    KG(1:3,1:3) = lA*HAB*kB_column*HAB'*lA';      % KAA
    KG(4:6,1:3) = -lB*kB_column*HAB'*lA';      % KBA
    KG(1:3,4:6) = -lA*HAB*kB_column*lB';      % KAB
    KG(4:6,4:6) = lB*kB_column*lB';      % KBB
    
    else        % for curved members
        
        KG(1:3,1:3) = lA*TAB*kB_curved*TAB'*lA';    % KAA
        KG(4:6,1:3) = -lB*kB_curved*TAB'*lA';    % KBA
        KG(1:3,4:6) = -lA*TAB*kB_curved*lB';    % KAB
        KG(4:6,4:6) = lB*kB_curved*lB';    % KBB
        
        end
        
        node1 = t(1,ii);
        node2 = t(2,ii);
        
        AA = node1*3-2 : node1*3;  % AA contains the dofs for the node A of the member
        BB = node2*3-2 : node2*3;  % BB contains the dofs for the node B of the member
        
        % Adding the element stiffness matrix coefficients in the global
        % stiffness matrix of the entire structure
        
        KKG(AA,AA) = KKG(AA,AA) + KG(1:3,1:3);
        KKG(BB,AA) = KKG(BB,AA) + KG(4:6,1:3);
        KKG(AA,BB) = KKG(AA,BB) + KG(1:3,4:6);
        KKG(BB,BB) = KKG(BB,BB) + KG(4:6,4:6);
        
        end
        % Eliminate rows and columns corresponding to constrained DOFs
        Force([1:3,28:30],:) = [];
        KK = KKG([1,3:6,10:12,16:19,21],[1,3:6,10:12,16:19,21]);
        displ = KK\Force;
        %Displacements of unconstrained nodes
        ux = displ([4,7,10,13,16,19,22,25]);
        uy = displ([11,14,17,20]);
        phi = displ([6,9,12,15,18,21,24]);
        % Table of nodal displacements for the structure
        disp(['     Node       '   '  ux [in]       '   'uy [in]         '   'phi [rad]'])
        disp([[2;4;6;7] , ux ,  [0;uy;0] ,    phi])
        % Displacement components for all nodes
        displALL_x = [ux(1),ux(2),0,ux(3),0,ux(4),ux(5)];
        displALL_y = [0,uy(1),0,uy(2),0,uy(3),0];
        n_coord = [000 120 120 195 270 270 390;...
            035 075 000 085 000 075 035].*12;
        xcoord = n_coord(1,:);
        ycoord = n_coord(2,:);
        scale = 100;
        figure(1), plot(xcoord,ycoord,'ks')
        hold on
        plot(xcoord + scale.*displALL_x , ycoord + scale.*displALL_y,'ro')
        legend('undeformed nodes','deformed nodes')
        xlabel('X axis [in]')
        ylabel('Y axis [in]')
        axis equal
        grid on

My current error

>> KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA'; % KAA

KG(4:6,1:3) = -lB*kB_straight*HAB'*lA'; % KBA

KG(1:3,4:6) = -lA*HAB*kB_straight*lB'; % KAB

KG(4:6,4:6) = lB*kB_straight*lB'

Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

[Walter Roberson]

kB_straight is an expression in A, Ay, Iz, and L.

The right hand side of the assignment involves kB_straight, so it is a symbolic expression involving A, Ay, Iz, and L.

A symbolic expression involving unresolved variables cannot be assigned to a numeric location.

[Walter Roberson]

Force and KK end up very different sizes; you cannot possibly \ between the two of them.

format short g
ndofn = 3; % number of degrees of freedom per node
% each number corresponds to a force for the corresponding DOF
Force = [0 0 0,  0 0 0,  0 0 0,  0 -1000 0,  0 -1000 0,  0 -1000 0,  0 -1000 0,  0 0 0,  0 0 0,  0 0 0]';
% Angles between local and global reference systems
% alphaA contains angles at the first node A
alphaA = [pi/2  pi/2  0  0  0  0  0  0    -pi/2];
% alphaB contains angles at the first node B
alphaB = [pi/2  0     0  0  0  0  0 -pi/2 -pi/2];
% distances along X and Y directions between origins of local reference systems
% in A and B; to be used in the tranlsation matrix
dX = [240 240 240 240 240 240 240 240 240]; %in
dY = [0  -240 0   0   0   0   0  -240 0]; %in
% Connectivity matrix: each column i shows the nodes that element i is connected to;
% For example, element 1 is connected to the structures nodes 1 and 2;
% In other words A==1  B==2 for the first element which is represented by
% the first column
t = [1 2 3 4 5 7 8 9;
     2 3 4 5 6 8 9 10];
%for straight
syms A Ay E G Iz L x real
%translation matrix for 2D element -> simplification of 3D case
HxB = [1 0 0;
    0 1 0;
    0 240 1];
HxB_T = HxB'; % transpose
% applies to linear portions of bridge
% fuq is cross sectional unit flexibility matrix
fuq_straight = [1/(A*E) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
fux_straight = HxB_T*fuq_straight*HxB;
fBA_straight = int(fux_straight, x, 0, L);
% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_straight = subs(inv(fBA_straight), {sym('G'), sym('E')}, {G, E})

kB_straight = 

%for columns
syms A1 Ay1 E G Iz1 L1 x real
%translation matrix for 2D element -> simplification of 3D case
HxBc = [1 0 0;
    0 1 0;
    0 240 1];
HxB_Tc = HxBc'; % transpose
% applies to columns
% fuq is cross sectional unit flexibility matrix
fuq_column = [1/(A1*E) 0 0;
    0 1/(G*Ay1) 0;
    0 0 1/(E*Iz1)];
fux_column = HxB_Tc*fuq_column*HxBc;
fBA_column = int(fux_column, x, 0, L1);
% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_column = double(subs(inv(fBA_column), {sym('A1'), sym('Ay1'), sym('G'), sym('E'), sym('Iz1'), sym('L1')}, {A1, Ay1, G, E, Iz1, L1}));

%for curved
syms alpha theta r E G A Ay Iz real
c = cos(alpha);
s = sin(alpha);
% rotation operator
lambda_QB_T = [c s 0;
    -s c 0;
    0 0 1];
lambda_QB = lambda_QB_T';
% translation operator accounts also for the different (rotated) reference
% systems at different cross sections in the curved member
HQB =   [1 0 0;
    0 1 0;
    -240 240 1];
HQB_T = HQB';
% unitary flexibility matrix
fuq = [1/(E*A) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
TQB = lambda_QB_T*HQB; % Translation and rotation operator
TQB_T = TQB'; % Transpose
fux = TQB_T * fuq * TQB.*r;
fBA_curved = int(fux, alpha, 0, theta); %Flexibility matrix of curved member
theta = pi/2;
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
r = 240; %in
c1 = cos(theta);
s1 = sin(theta);
% Stiffness matrix of cantilever curved element
kB_curved = eval(inv(fBA_curved));
%global matrice
%nnodes = length(t)+1;
nnodes = max(t(:));
% Initializing the global stiffness matrix of the entire structure
% Its dimension is square with nnodes*ndofn rows and columns
KKG = zeros(nnodes*ndofn, nnodes*ndofn, 'sym');
% Initializing the global stiffness matrix of the one element
% Its dimension is square with 2*ndofn rows and columns
KG = zeros(2*ndofn, 2*ndofn, 'sym');
%Initialize Matrices
HAB = zeros(3,3);
lA = zeros(3,3);
lB = zeros(3,3);
for ii=1:length(t)
    HAB = [1 0 0;
        0 1 0;
        -dY(ii) dX(ii) 1];    % translation operator for straight members
    
    cA = cos(alphaA(ii));
    sA = sin(alphaA(ii));
    cB = cos(alphaB(ii));
    sB = sin(alphaB(ii));
    
    lA = [cA -sA 0;
        sA  cA 0;
        0   0  1];
    
    lB = [cB -sB 0;
        sB  cB 0;
        0   0  1];
    
    angle = pi/2;
    
    TAB = [  cos(angle)    sin(angle)    0;    % rotation/translation operator for curved members
        -sin(angle)    cos(angle)    0     % local ref systems in A and B are off by angle
        -r*(1-cos(angle)) r*sin(angle)  1];
    
    disp(['element= ',num2str(ii)])
    disp(['   '])
    
    if ii==1 || ii==6   % for straight components
        
        KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_straight*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_straight*lB';      % KAB
        KG(4:6,4:6) = lB*kB_straight*lB';      % KBB
        
    elseif ii==2 || ii==5     % for all columns
        
        KG(1:3,1:3) = lA*HAB*kB_column*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_column*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_column*lB';      % KAB
        KG(4:6,4:6) = lB*kB_column*lB';      % KBB
        
    else        % for curved members
        
        KG(1:3,1:3) = lA*TAB*kB_curved*TAB'*lA';    % KAA
        KG(4:6,1:3) = -lB*kB_curved*TAB'*lA';    % KBA
        KG(1:3,4:6) = -lA*TAB*kB_curved*lB';    % KAB
        KG(4:6,4:6) = lB*kB_curved*lB';    % KBB
        
    end
    
    node1 = t(1,ii);
    node2 = t(2,ii);
    
    AA = node1*3-2 : node1*3;  % AA contains the dofs for the node A of the member
    BB = node2*3-2 : node2*3;  % BB contains the dofs for the node B of the member
    
    % Adding the element stiffness matrix coefficients in the global
    % stiffness matrix of the entire structure
    
    KKG(AA,AA) = KKG(AA,AA) + KG(1:3,1:3);
    KKG(BB,AA) = KKG(BB,AA) + KG(4:6,1:3);
    KKG(AA,BB) = KKG(AA,BB) + KG(1:3,4:6);
    KKG(BB,BB) = KKG(BB,BB) + KG(4:6,4:6);
    
end

element= 1

element= 2

element= 3

element= 4

element= 5

element= 6

element= 7

element= 8

% Eliminate rows and columns corresponding to constrained DOFs
Force([1:3,28:30],:) = [];
KK = KKG([1,3:6,10:12,16:19,21],[1,3:6,10:12,16:19,21]);
whos KK Force

  Name        Size            Bytes  Class     Attributes

  Force      24x1               192  double              
  KK         13x13                8  sym                 

displ = KK\Force;

Error using \
Invalid operands.

%Displacements of unconstrained nodes
ux = displ([4,7,10,13,16,19,22,25]);
uy = displ([11,14,17,20]);
phi = displ([6,9,12,15,18,21,24]);
% Table of nodal displacements for the structure
disp(['     Node       '   '  ux [in]       '   'uy [in]         '   'phi [rad]'])
disp([[2;4;6;7] , ux ,  [0;uy;0] ,    phi])
% Displacement components for all nodes
displALL_x = [ux(1),ux(2),0,ux(3),0,ux(4),ux(5)];
displALL_y = [0,uy(1),0,uy(2),0,uy(3),0];
n_coord = [000 120 120 195 270 270 390;...
    035 075 000 085 000 075 035].*12;
xcoord = n_coord(1,:);
ycoord = n_coord(2,:);
scale = 100;
figure(1), plot(xcoord,ycoord,'ks')
hold on
plot(xcoord + scale.*displALL_x , ycoord + scale.*displALL_y,'ro')
legend('undeformed nodes','deformed nodes')
xlabel('X axis [in]')
ylabel('Y axis [in]')
axis equal
grid on

[Vruddhi]

axes
clear
clc
close all
format short g
ndofn = 3; % number of degrees of freedom per node
% each number corresponds to a force for the corresponding DOF
Force = [0 0 0  0 0 0  0 0 0  0 -1 0  0 -1 0  0 -1 0  0 -1 0  0 0 0  0 0 0  0 0 0]';
% Angles between local and global reference systems
% alphaA contains angles at the first node A
alphaA = [pi/2  pi/2  0  0  0  0  0  0    -pi/2];
% alphaB contains angles at the first node B
alphaB = [pi/2  0     0  0  0  0  0 -pi/2 -pi/2];
% distances along X and Y directions between origins of local reference systems
% in A and B; to be used in the tranlsation matrix
dX = [240 240 240 240 240 240 240 240 240]; %in
dY = [0  -240 0   0   0   0   0  -240 0]; %in
% Connectivity matrix: each column i shows the nodes that element i is connected to;
% For example, element 1 is connected to the structures nodes 1 and 2;
% In other words A==1  B==2 for the first element which is represented by
% the first column
t = [1 2 3 4 5 7 8 9 10;
    2 3 4 5 6 8 9 10 11];
%for straight
syms L x A E G Ay Iz real
%translation matrix for 2D element -> simplification of 3D case
HxB = [1 0 0;
    0 1 0;
    0 240 1];
HxB_T = HxB'; % transpose
% applies to linear portions of bridge
% fuq is cross sectional unit flexibility matrix
fuq_straight = [1/(A*E) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
fux_straight = HxB_T*fuq_straight*HxB;
fBA_straight = int(fux_straight, x, 0, L);
% Numerical values given and from SAP
E = 30000; %ksi
nu = 0.29;
G = E/(2*(1+nu)); %ksi
A = 53.407; %in^2
Iz = 1936; %in^4
Ay = 53.407; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_straight = subs(inv(fBA_straight),{sym('G'), sym('E')}, {G, E});
%for columns
syms L1 x A E G Ay Iz real
%translation matrix for 2D element -> simplification of 3D case
HxBc = [1 0 0;
    0 1 0;
    0 240 1];
HxB_Tc = HxBc'; % transpose
% applies to columns
% fuq is cross sectional unit flexibility matrix
fuq_column = [1/(A*E) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
fux_column = HxB_Tc*fuq_column*HxBc;
fBA_column = int(fux_column, x, 0, L1);
% Numerical values given and from SAP
E = 30000; %ksi
nu = 0.29;
G = E/(2*(1+nu)); %ksi
A = 53.407; %in^2
Iz = 1936; %in^4
Ay = 53.407; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_column =double(subs(inv(fBA_column), {sym('A'), sym('Ay'), sym('G'), sym('E'), sym('Iz'), sym('L1')}, {A, Ay, G, E, Iz, L1}));
%for curved
syms alpha theta r E G A Ay Iz real
c = cos(alpha);
s = sin(alpha);
% rotation operator
lambda_QB_T = [c s 0
    -s c 0;
    0 0 1];
lambda_QB = lambda_QB_T';
% translation operator accounts also for the different (rotated) reference
% systems at different cross sections in the curved member
HQB =   [1 0 0;
    0 1 0;
    -240 240 1];
HQB_T = HQB';
% unitary flexibility matrix
fuq = [1/(E*A) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
TQB = lambda_QB_T*HQB; % Translation and rotation operator
TQB_T = TQB'; % Transpose
fux = TQB_T * fuq * TQB.*r;
fBA_curved = int(fux, alpha, 0, theta); %Flexibility matrix of curved member
theta = pi/2;
E = 30000; %ksi
nu = 0.29;
G = E/(2*(1+nu)); %ksi
A = 53.407; %in^2
Iz = 1936; %in^4
Ay = 53.407; %in^2
r = 240; %in
c1 = cos(theta);
s1 = sin(theta);
% Stiffness matrix of cantilever curved element
kB_curved = subs(inv(fBA_curved));
%global matrice
%nnodes = length(t)+1;
nnodes = max(t(:));
% Initializing the global stiffness matrix of the entire structure
% Its dimension is square with nnodes*ndofn rows and columns
KKG = zeros(nnodes*ndofn,nnodes*ndofn,'sym');
% Initializing the global stiffness matrix of the one element
% Its dimension is square with 2*ndofn rows and columns
KG = zeros(2*ndofn,2*ndofn,'sym');
%Initialize Matrices
HAB = zeros(3,3);
lA = zeros(3,3);
lB = zeros(3,3);
for ii=1:length(t)
    HAB = [1 0 0;
        0 1 0;
        -dY(ii) dX(ii) 1];    % translation operator for straight members
    
    cA = cos(alphaA(ii));
    sA = sin(alphaA(ii));
    cB = cos(alphaB(ii));
    sB = sin(alphaB(ii));
    
    lA = [cA sA 0;
        -sA cA 0;
        0   0  1];
    
    lB = [cB sB 0;
        -sB cB 0;
        0   0  1];
    
    angle = pi/2;
    
    TAB = [cos(angle)    sin(angle)    0;    % rotation/translation operator for curved members
        -sin(angle)    cos(angle)    0     % local ref systems in A and B are off by angle
        -r*(1-cos(angle)) r*sin(angle)  1];
    
    disp(['element= ',num2str(ii)])
    disp(['   '])
    
    if ii==3 || ii==4  || ii==5 || ii==6 || ii==7  % for straight components
        
        KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_straight*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_straight*lB';      % KAB
        KG(4:6,4:6) = lB*kB_straight*lB';      % KBB
        
    elseif ii==1 || ii==9     % for all columns
        
        KG(1:3,1:3) = lA*HAB*kB_column*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_column*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_column*lB';      % KAB
        KG(4:6,4:6) = lB*kB_column*lB';      % KBB
        
    else        % for curved members
        
        KG(1:3,1:3) = lA*TAB*kB_curved*TAB'*lA';    % KAA
        KG(4:6,1:3) = -lB*kB_curved*TAB'*lA';    % KBA
        KG(1:3,4:6) = -lA*TAB*kB_curved*lB';    % KAB
        KG(4:6,4:6) = lB*kB_curved*lB';    % KBB
        
    end
    
    node1 = t(1,ii);
    node2 = t(2,ii);
    
    AA = node1*3-2 : node1*3;  % AA contains the dofs for the node A of the member
    BB = node2*3-2 : node2*3;  % BB contains the dofs for the node B of the member
    
    % Adding the element stiffness matrix coefficients in the global
    % stiffness matrix of the entire structure
    
    KKG(AA,AA) = KKG(AA,AA) + KG(1:3,1:3);
    KKG(BB,AA) = KKG(BB,AA) + KG(4:6,1:3);
    KKG(AA,BB) = KKG(AA,BB) + KG(1:3,4:6);
    KKG(BB,BB) = KKG(BB,BB) + KG(4:6,4:6);
    
end

element= 1

element= 2

element= 3

element= 4

element= 5

element= 6

element= 7

element= 8

element= 9

% Eliminate rows and columns corresponding to constrained DOFs
Force([1:3,28:30],:) = [];
KK = KKG([4:27],[4:27]);
displ = KK\Force;
%Displacements of unconstrained nodes
ux = displ([2,3,8,9]);
uy = displ([2,3,4,5,6,7,8]);
phi = displ([2,3,4,5,6,7,8]);
% Table of nodal displacements for the structure
disp(['     Node       '   '  ux [in]       '   'uy [in]         '   'phi [rad]'])

     Node         ux [in]       uy [in]         phi [rad]

disp([[2;3;4;5;6;7;8] , ux ,  [0;uy;0] ,    phi])

Error using sym/cat>checkDimensions
CAT arguments dimensions not consistent.

Error in sym/cat>catMany (line 33)
[resz, ranges] = checkDimensions(sz,dim);

Error in sym/cat (line 25)
ySym = catMany(dim, args);

Error in sym/horzcat (line 19)
ySym = cat(2,args{:});

% Displacement components for all nodes
displALL_x = [ux(2),ux(3),uy(7),ux(8)];
displALL_y = [uy(4),uy(5),uy(6),uy(7)];
n_coord = [000 000 020 040 060 080 100 120 140 140;...
    000 020 040 040 040 040 040 040 020 000].*12;
xcoord = n_coord(1,:);
ycoord = n_coord(2,:);
scale = 100;
figure(1), plot(xcoord,ycoord,'ks')
hold on
plot(xcoord + scale.*displALL_x , ycoord + scale.*displALL_y,'ro')
legend('undeformed nodes','deformed nodes')
xlabel('X axis [in]')
ylabel('Y axis [in]')
axis equal
grid on

I tryed doing that, I am new to matlab , I need displacement values from this code . And want a plot . Which i am unable to run. kindly guide in urgent.

[Walter Roberson]

I can add arbitary padding to ux to get it to fit the table, but I can't do anything about the fact that you try to index element 7 of a vector of length 4.

format short g
ndofn = 3; % number of degrees of freedom per node
% each number corresponds to a force for the corresponding DOF
Force = [0 0 0  0 0 0  0 0 0  0 -1 0  0 -1 0  0 -1 0  0 -1 0  0 0 0  0 0 0  0 0 0]';
% Angles between local and global reference systems
% alphaA contains angles at the first node A
alphaA = [pi/2  pi/2  0  0  0  0  0  0    -pi/2];
% alphaB contains angles at the first node B
alphaB = [pi/2  0     0  0  0  0  0 -pi/2 -pi/2];
% distances along X and Y directions between origins of local reference systems
% in A and B; to be used in the tranlsation matrix
dX = [240 240 240 240 240 240 240 240 240]; %in
dY = [0  -240 0   0   0   0   0  -240 0]; %in
% Connectivity matrix: each column i shows the nodes that element i is connected to;
% For example, element 1 is connected to the structures nodes 1 and 2;
% In other words A==1  B==2 for the first element which is represented by
% the first column
t = [1 2 3 4 5 7 8 9 10;
    2 3 4 5 6 8 9 10 11];
%for straight
syms L x A E G Ay Iz real
%translation matrix for 2D element -> simplification of 3D case
HxB = [1 0 0;
    0 1 0;
    0 240 1];
HxB_T = HxB'; % transpose
% applies to linear portions of bridge
% fuq is cross sectional unit flexibility matrix
fuq_straight = [1/(A*E) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
fux_straight = HxB_T*fuq_straight*HxB;
fBA_straight = int(fux_straight, x, 0, L);
% Numerical values given and from SAP
E = 30000; %ksi
nu = 0.29;
G = E/(2*(1+nu)); %ksi
A = 53.407; %in^2
Iz = 1936; %in^4
Ay = 53.407; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_straight = subs(inv(fBA_straight),{sym('G'), sym('E')}, {G, E});
%for columns
syms L1 x A E G Ay Iz real
%translation matrix for 2D element -> simplification of 3D case
HxBc = [1 0 0;
    0 1 0;
    0 240 1];
HxB_Tc = HxBc'; % transpose
% applies to columns
% fuq is cross sectional unit flexibility matrix
fuq_column = [1/(A*E) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
fux_column = HxB_Tc*fuq_column*HxBc;
fBA_column = int(fux_column, x, 0, L1);
% Numerical values given and from SAP
E = 30000; %ksi
nu = 0.29;
G = E/(2*(1+nu)); %ksi
A = 53.407; %in^2
Iz = 1936; %in^4
Ay = 53.407; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_column =double(subs(inv(fBA_column), {sym('A'), sym('Ay'), sym('G'), sym('E'), sym('Iz'), sym('L1')}, {A, Ay, G, E, Iz, L1}));
%for curved
syms alpha theta r E G A Ay Iz real
c = cos(alpha);
s = sin(alpha);
% rotation operator
lambda_QB_T = [c s 0
    -s c 0;
    0 0 1];
lambda_QB = lambda_QB_T';
% translation operator accounts also for the different (rotated) reference
% systems at different cross sections in the curved member
HQB =   [1 0 0;
    0 1 0;
    -240 240 1];
HQB_T = HQB';
% unitary flexibility matrix
fuq = [1/(E*A) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];
TQB = lambda_QB_T*HQB; % Translation and rotation operator
TQB_T = TQB'; % Transpose
fux = TQB_T * fuq * TQB.*r;
fBA_curved = int(fux, alpha, 0, theta); %Flexibility matrix of curved member
theta = pi/2;
E = 30000; %ksi
nu = 0.29;
G = E/(2*(1+nu)); %ksi
A = 53.407; %in^2
Iz = 1936; %in^4
Ay = 53.407; %in^2
r = 240; %in
c1 = cos(theta);
s1 = sin(theta);
% Stiffness matrix of cantilever curved element
kB_curved = subs(inv(fBA_curved));
%global matrice
%nnodes = length(t)+1;
nnodes = max(t(:));
% Initializing the global stiffness matrix of the entire structure
% Its dimension is square with nnodes*ndofn rows and columns
KKG = zeros(nnodes*ndofn,nnodes*ndofn,'sym');
% Initializing the global stiffness matrix of the one element
% Its dimension is square with 2*ndofn rows and columns
KG = zeros(2*ndofn,2*ndofn,'sym');
%Initialize Matrices
HAB = zeros(3,3);
lA = zeros(3,3);
lB = zeros(3,3);
for ii=1:length(t)
    HAB = [1 0 0;
        0 1 0;
        -dY(ii) dX(ii) 1];    % translation operator for straight members
    
    cA = cos(alphaA(ii));
    sA = sin(alphaA(ii));
    cB = cos(alphaB(ii));
    sB = sin(alphaB(ii));
    
    lA = [cA sA 0;
        -sA cA 0;
        0   0  1];
    
    lB = [cB sB 0;
        -sB cB 0;
        0   0  1];
    
    angle = pi/2;
    
    TAB = [cos(angle)    sin(angle)    0;    % rotation/translation operator for curved members
        -sin(angle)    cos(angle)    0     % local ref systems in A and B are off by angle
        -r*(1-cos(angle)) r*sin(angle)  1];
    
    disp(['element= ',num2str(ii)])
    disp(['   '])
    
    if ii==3 || ii==4  || ii==5 || ii==6 || ii==7  % for straight components
        
        KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_straight*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_straight*lB';      % KAB
        KG(4:6,4:6) = lB*kB_straight*lB';      % KBB
        
    elseif ii==1 || ii==9     % for all columns
        
        KG(1:3,1:3) = lA*HAB*kB_column*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_column*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_column*lB';      % KAB
        KG(4:6,4:6) = lB*kB_column*lB';      % KBB
        
    else        % for curved members
        
        KG(1:3,1:3) = lA*TAB*kB_curved*TAB'*lA';    % KAA
        KG(4:6,1:3) = -lB*kB_curved*TAB'*lA';    % KBA
        KG(1:3,4:6) = -lA*TAB*kB_curved*lB';    % KAB
        KG(4:6,4:6) = lB*kB_curved*lB';    % KBB
        
    end
    
    node1 = t(1,ii);
    node2 = t(2,ii);
    
    AA = node1*3-2 : node1*3;  % AA contains the dofs for the node A of the member
    BB = node2*3-2 : node2*3;  % BB contains the dofs for the node B of the member
    
    % Adding the element stiffness matrix coefficients in the global
    % stiffness matrix of the entire structure
    
    KKG(AA,AA) = KKG(AA,AA) + KG(1:3,1:3);
    KKG(BB,AA) = KKG(BB,AA) + KG(4:6,1:3);
    KKG(AA,BB) = KKG(AA,BB) + KG(1:3,4:6);
    KKG(BB,BB) = KKG(BB,BB) + KG(4:6,4:6);
    
end

element= 1

element= 2

element= 3

element= 4

element= 5

element= 6

element= 7

element= 8

element= 9

% Eliminate rows and columns corresponding to constrained DOFs
Force([1:3,28:30],:) = [];
KK = KKG([4:27],[4:27]);
displ = KK\Force;
%Displacements of unconstrained nodes
ux = displ([2,3,8,9]);
uy = displ([2,3,4,5,6,7,8]);
phi = displ([2,3,4,5,6,7,8]);
whos ux uy phi

  Name      Size            Bytes  Class    Attributes

  phi       7x1                 8  sym                
  ux        4x1                 8  sym                
  uy        7x1                 8  sym                

% Table of nodal displacements for the structure
disp(['     Node       '   '  ux [in]       '   'uy [in]         '   'phi [rad]'])

     Node         ux [in]       uy [in]         phi [rad]

disp(vpa([[2;3;4;5;6;7;8] , [ux;0;0;0] ,  uy ,    phi],5))

% Displacement components for all nodes
displALL_x = [ux(2),ux(3),uy(7),ux(8)];

Index exceeds the number of array elements. Index must not exceed 4.

Error in indexing (line 968)
R_tilde = builtin('subsref',L_tilde,Idx);

displALL_y = [uy(4),uy(5),uy(6),uy(7)];
n_coord = [000 000 020 040 060 080 100 120 140 140;...
    000 020 040 040 040 040 040 040 020 000].*12;
xcoord = n_coord(1,:);
ycoord = n_coord(2,:);
scale = 100;
figure(1), plot(xcoord,ycoord,'ks')
hold on
plot(xcoord + scale.*displALL_x , ycoord + scale.*displALL_y,'ro')
legend('undeformed nodes','deformed nodes')
xlabel('X axis [in]')
ylabel('Y axis [in]')
axis equal
grid on

[Vruddhi]

How should i correct the above I tryed understand your comment but couldnt implement and make changes in my code. Please help

[Walter Roberson]

You have

ux = displ([2,3,8,9]);

which makes ux a vector of length 4.

You tried to do

disp([[2;3;4;5;6;7;8] , ux ,  [0;uy;0] ,    phi])

which would require that ux is a vector of length 7.

You tried to do

displALL_x = [ux(2),ux(3),uy(7),ux(8)];

which would have required that ux is a vector of length at least 8.

These definitions of ux are incompatible.

You also have

uy = displ([2,3,4,5,6,7,8]);

which makes uy have a vector of length 7. It is questionable whether it makes sense to extract the same elements of displ for both ux and uy.

Getting back to

disp([[2;3;4;5;6;7;8] , ux ,  [0;uy;0] ,    phi])

this requires that uy is a vector of length 5, which is incompatible with the way it is initialized.

I cannot tell you how to repair these problems; you need to re-think what ux and uy mean.

Question: in

displALL_x = [ux(2),ux(3),uy(7),ux(8)];

are you sure that you want uy(7) in a vector that is intended to be ALL_x ?