桁架矩阵位移法matlab计算程序

%% Yjj 211908214006
clear all;
clc;
jd=input('请输入节点数量');
free=input('请输入自由度');
gj=input('请输入杆件数量');
P=zeros(free,1);
P=input('请输入外荷载矩阵（[x;y;z]）');
d=zeros(free,1);
member=zeros(1,4,gj);
K=zeros(4,4,gj);
for i=1:1:gj
angle(i)=input(sprintf('请输入第%d杆件角度（角度制）',i));
E(i)=input(sprintf('请输入第%d杆件弹性模量（N/mm2）',i));
A(i)=input(sprintf('请输入第%d杆件截面面积（mm2）',i));
L(i)=input(sprintf('请输入第%d杆件长度（mm）',i));
member(:,:,i)=input(sprintf('请输入第%d杆件定位向量([1 2 3 4])',i));
end
for i=1:1:gj
K(:,:,i)=[cosd(angle(i))*cosd(angle(i)) cosd(angle(i))*sind(angle(i)) -cosd(angle(i))*cosd(angle(i)) -cosd(angle(i))*sind(angle(i));
cosd(angle(i))*sind(angle(i)) sind(angle(i))*sind(angle(i)) cosd(angle(i))*sind(angle(i)) -sind(angle(i))*sind(angle(i));
-cosd(angle(i))*cosd(angle(i)) -cosd(angle(i))*sind(angle(i)) cosd(angle(i))*cosd(angle(i)) cosd(angle(i))*sind(angle(i));
-cosd(angle(i))*sind(angle(i)) -sind(angle(i))*sind(angle(i)) cosd(angle(i))*sind(angle(i)) sind(angle(i))*sind(angle(i))]*E(i)*A(i)/L(i);
end
Ss=zeros(free);
S=zeros(free);
for i=1:1:gj
for n=1:1:4
for j=1:1:4
if (member(1,n,i)<=free && member(1,j,i)<=free)
Ss(member(1,n,i),member(1,j,i))=K(n,j,i);
end
end
end
S=S+Ss ;
end
d=S\P
v=zeros(4,1,gj);
for i=1:1:gj
c=1;
for j=1:1:4
if(member(1,j,i)<free+1)
v(j,1,i)=d(c,1);
c=c+1;
end
end
end
T=zeros(4,4,gj);
for i=1:1:gj
T(:,:,i)=[cosd(angle(i)) sind(angle(i)) 0 0;
-sind(angle(i)) cosd(angle(i)) 0 0;
0 0 cosd(angle(i)) sind(angle(i));
0 0 -sind(angle(i)) cosd(angle(i))]; end
u=zeros(4,1,gj)
for i=1:1:gj
u(:,:,i)=T(:,:,i)*v(:,:,i);
end
k=zeros(4,4,gj);
for i=1:1:gj
k(:,:,i)=[1 0 -1 0;
0 0 0 0;
-1 0 1 0;
0 0 0 0]*E(i)*A(i)/L(i); end
Q=zeros(4,1,gj);
for i=1:1:gj
Q(:,:,i)=k(:,:,i)*u(:,:,i)
end
F=zeros(4,1,gj);
for i=1:1:gj
F(:,:,i)=T(:,:,i)'*Q(:,:,i);
end
zfl=jd*2;
r=zeros(zfl,1);
for i=1:1:gj
for j=1:1:4
r(member(1,j,i),1)=F(j,1,i);
end
end
zfl=zfl-free;
R=r(free+1:end)。