MATLAB矩阵位移法

for i=1n wk=qk(E,A,I,l(i)) wl=ql(a(i)); kk2=wkwl; kk1(,6i-56i)=kk2; kk=wl'kk2 k(,6i-56i)=kk; end qqq=max(max(r)); K=zeros(m,m); KR=zeros(qqq-m,m); for p=1n aa=r(p,); for i=16 for j=16 if aa(i)=m & aa(j)=m kkk=k(,6p-56p); K(aa(i),aa(j))=K(aa(i),aa(j))+kkk(i,j); end if aa(i)=m+1 & aa(j)=m kkk=k(,6p-56p); KR(aa(i)-m,aa(j))=KR(aa(i)-m,aa(j))+kkk(i,j); end end end end K KR Fell=zeros(1,qqq); for p=1n bb=r(p,); tt=ql(a(p)); ll=tt'Fg(,p); lll=ll'; for i=16 Fell(1,bb(i))=Fell(1,bb(i))+lll(1,i); end end Fel=-Fell'; Fej=Fed+Fel; Fejj=Fej'; Fe111=Fejj(1m); Fe1=Fe111' Fe222=Fejj(m+1qqq); Fe2=Fe222' d=KFe1 d1=d'; d1(m+1qqq)=0; for p=1n cc=r(p,); for i=16 dd(1,i)=d1(1,cc(i)); end p F=kk1(,6p-56p)dd'+Fg(,p) end FR=KRd-Fe2 定义函数qk function y=qk(E,A,I,l) y=[EAl 0 0 -EAl 0 0;0 12EIlll -6EIll 0 -12EIlll -6EIll;0 -6EIll 4EIl 0 6EIll 2EIl;-EAl 0 0 EAl 0 0;0 -12EIlll 6EIll 0 12EIlll 6EIll;0 -6EIll 2EIl 0 6EIll 4EIl]; 定义函数ql function yy=ql(a) yy=[cos(a) sin(a) 0 0 0 0;-sin(a) cos(a) 0 0 0 0;0 0 1 0 0 0;0 0 0 cos(a) sin(a) 0;0 0 0 -sin(a) cos(a) 0;0 0 0 0 0 1];