ない-A-数MATLAB

-1

clear 

% input problem mesh 
node = [ 0 0; 
     20 0; 
     20 15; 
      0 15]*12; 

conn = [1 2; 
     1 3; 
     2 4; 
     2 3; 
     4 3; 
     1 4]; 

% properties 
A = 1.0; 
E = 10e6; 

% boundary conditions 
P = 1000; 
f=zeros(8,1); 
f(6)=P; 

isol = [3 4 5 6 ]; %unconstrained dof 

%-----------------END OF INPUT---------------------- 

nn = size(node,1); %number of nodes 
ne = size(conn,1); %number of elements 
ndof = 2*nn; % number of dof 

K = zeros(ndof,ndof); %global stiffness matrix 
d = zeros(ndof, 1); %displacement vector 

for e=1:ne %loop through all the elements 

     n1=conn(e,1); 
     n2=conn(e,2); 

     x1=node(n1,1); y1=node(n1,2); % x and y coordinate for the first node 
     x2=node(n2,1); y2=node(n2,2); %x and y coordinate for the second node 

     L =sqrt( (x2-x1)^2 + (y2-y2)^2 ); %Length of the element 

     C = (x2-x1)/L; %Cosine 
     S = (y2-y1)/L; %Sine 
     C2 = C*C; % Cosine square 
     S2 = S*S; % Sine square 
     CS = C*S; 

     ke =(A*E/L)*[ C2 CS -C2 -CS; 
        CS S2 -CS -S2; 
        -C2 -CS C2 CS; 
        -CS -S2 CS S2 ]; %Stiffness matrix for element e 

     sctr = [2*n1-1 2*n1 2*n2-1 2*n2]; %Location where ke is to scatter 
              %in the global stiffness matrix 
              %scatter array 


     K(sctr,sctr) = K(sctr,sctr) + ke; 
end 

%solve for the nodal displacement, d 
d(isol) = K(isol,isol)\f(isol); 
fprintf('\n N o d a l D i s p l a c e m e n t\n') 
fprintf(' NID  X-Force Y-Force\n') 
for i=1:nn 
     fprintf('%5d %8.3e  %8.3e\n',i,d(2*i-1),d(2*i)) 
end 

%compute the reaction forces 
f = K*d; 

fprintf('\n E x t e r n a l F o r c e s\n') 
fprintf(' NID  X-FORCE Y-FORCE\n') 
for i=1:nn 
     fprintf('%5d %8.3e  %8.3e\n',i,f(2*i-1),f(2*i)) 
end 

%compute the element stress 
fprintf('\n E l e m e n t S t r e s s\n') 
fprintf('EID  STRAIN  STRESS\n') 
for e = 1:ne 

     n1 = conn(e,1); 
     n2 = conn(e,2); 

     x1 = node(n1,1); y1 = node(n1,2); 
     x2 = node(n2,1); y2 = node(n2,2); 

     L = sqrt((x2-x1)^2 + (y2-y1)^2); 
     C =(x2-x1)/L; 
     S =(y2-y1)/L; 

     B = 1/L*[-C -S C S]; %the element b matrix 

     sctr = [2*n1-1 2*n1 2*n2-1 2*n2 ]; 
     strain = B*d(sctr); 
     stress = E*strain; 

     fprintf('% 5d  %8.3e %8.3e\n',i,strain,stress) 

end

ここCの値が0であるとSの値はCS = C * Sが出力としてNaNを与えているが、1 C2 = C * C、S2 = S * S、です。 NaNの値がS2,C2,CSの算術演算で得られた実際の値になるのを防ぐにはどうすればいいですか？これは、MATLABで行われる有限要素解析コードです。ない-A-数MATLAB