1、此文是基于matlab的受轴弯剪的梁单元的结构求解，代码如下：

% 基于matlab的受轴弯剪的梁单元的结构求解；
% 输入量为：节点数目，坐标； 单元节点编号，自由度编号; 已知力，位移； 未知位移；
clc      
clear all 

num_nod=3;  %节点数目；
ele_nod=[1 2;2 3]; %单元节点编号；
nod_coor=[0 0;12 0; 24.73 12.73] ; %节点坐标，向右为x轴，向下为y轴；
ele_dof=[1 2 3 4 5 6;4 5 6 7 8 9]; %单元自由度编号，对应于子块搬家对应的家；
force=zeros(3*num_nod,1);
displacement=zeros(3*num_nod,1);
stiffness=zeros(3*num_nod);
ele_num=size(ele_nod,1);
force(5) = 1000000000.0; %已知节点力；
displacement(1)=0.0;  %约束条件；
displacement(2)=0.0;
displacement(3)=0.0;
displacement(7)=0.0;
displacement(8)=0.0;
displacement(9)=0.0;
unknown_f_a=[4;5;6]; %未知节点位移；

% 截面基本参数
b=0.5; 
h=0.5; 
E=2.1e11;


%初使化所要的矩阵，F K 位移矩阵；

A=b*h;
I=b*h^3/12; 

for e = 1 : ele_num
    L(e) = sqrt(  (nod_coor(ele_nod(e,1),1)-nod_coor(ele_nod(e,2),1))^2+ ...
                  (nod_coor(ele_nod(e,1),2)-nod_coor(ele_nod(e,2),2))^2   );
      C=(nod_coor(ele_nod(e,2),1)-nod_coor(ele_nod(e,1),1))/L(e);     
      S=(nod_coor(ele_nod(e,2),2)-nod_coor(ele_nod(e,1),2))/L(e);      
      t=[  C S 0 0 0 0; 
          -S C 0 0 0 0; 
           0 0 1 0 0 0; 
           0 0 0 C S 0; 
           0 0 0 -S C 0; 
           0 0 0 0 0 1];      %坐标转换矩阵    
      k11=A*E/L(e);
      k22=12*E*I/L(e)^3;
      k23=6*E*I/L(e)^2;
      k33=4*E*I/L(e);
      k36=2*E*I/L(e);
      k_local = [k11 0 0 -k11 0 0;0 k22 k23 0 -k22 k23;0 k23 k33 0 -k23 k36;-k11 0 0 k11 0 0;0 -k22 -k23 0 k22 -k23; 0 k23 k36 0 -k23 k33];
      k=t'*k_local*t; 
      
      ele_dof_a = ele_dof(e,:);      
      for i = 1 : 6
          for j = 1 : 6
          
              stiffness(ele_dof_a(1,i),ele_dof_a(1,j))=...
              stiffness(ele_dof_a(1,i),ele_dof_a(1,j))+k(i,j);
          end
      end       
              
end

%因为此时的K是奇异的，那么要消除对应的值；

for i = 1:size(unknown_f_a)
   dis_new(i,1)=displacement(unknown_f_a(i,1),1);
   force_new(i,1)=force(unknown_f_a(i,1),1);
end

for i = 1:size(unknown_f_a)
    for j =1:size(unknown_f_a)
        stiff_new(i,j) = stiffness(unknown_f_a(i,1),unknown_f_a(j,1));
    end
end

%消除K的奇异性后，求解所对应自由度处的位移；
dis_new = inv(stiff_new)*force_new;
for i = 1 : size(unknown_f_a)
fprintf('对应的第%d个自由度的未知位移大小为 %d；\n',unknown_f_a(i),dis_new(i))
end

for i = 1 :size(unknown_f_a)

    displacement(unknown_f_a(i,1),1) = dis_new(i,1);
end


for e = 1 :ele_num
    
    x = [nod_coor(ele_nod(e,1),1) nod_coor(ele_nod(e,2),1)];
     y = [nod_coor(ele_nod(e,1),2) nod_coor(ele_nod(e,2),2)];
     plot(x,-y,'b')
     hold on     
end

for i=1:num_nod

    nod_coor_def(i,1)=nod_coor(i,1)+displacement(3*i-2,1);
    nod_coor_def(i,2)=nod_coor(i,2)+displacement(3*i-1,1);
end

for e=1:ele_num
    x=[nod_coor_def(ele_nod(e,1),1) nod_coor_def(ele_nod(e,2),1)];
    y=[nod_coor_def(ele_nod(e,1),2) nod_coor_def(ele_nod(e,2),2)];
plot(x,-y,'r')
hold on
end