matlab程序错误

[skz2005]

麻烦高手帮我看看，这个程序怎么改
function SECOND
%多层结构在顶部受到简谐振动荷载作用下的响应

M=[1,2,3,4];
M=diag(M);
K=[800,-800,0,0;-800,2400,-2400,0;0,-2400,4800,-4800;0,0,-4800,8000];
C=0.05*M+0.02*K;
beta=1/4;
dt=0.01;
K1=K+C/(2*beta*dt)+M/(beta*dt^2);
t=0;
u=[0;0;0;0];
u1=[0;0;0;0];
P=[cos(6.647*t);0;0;0];
u2=inv(K)*(P-C*u-K*u1);
A=[0,0,0,0];
B=[0,0,0,0];
D=[0,0,0,0];
for i=1:1000
t=t+dt;
PN=[cos(6.647*t);0;0;0];
dP=PN-P;
P=PN;
dP1=dP+M*(u1/(beta*dt)+u2/(2*beta))+C*(u1/(2*beta)+u2*dt*(1/(4*beta)-1));
du=inv(K1)*dP1;
du1=du/(2*beta*dt)-u1/(2*beta)-u2*dt*(0.25/beta-1);
du2=du/(beta*dt^2)-u1/(beta*dt)-u2/(2*beta);
u=du+u;
u1=du1+u1;
u2=du2+u2;
m=u';
n=u1';
q=u2';
A=cat(1,A,m);
B=cat(1,B,n);
D=cat(1,D,q);
end;
fid=fopen('output u.dat','w');
fprintf(fid,'%12.4f\r\n',A);
fid=fopen('output u1.dat','w');
fprintf(fid,'%12.4f\r\n',B);
fid=fopen('output u2.dat','w');
fprintf(fid,'%12.4f\r\n',D);
end

出来的数据结果不对，不知道哪儿错了。。。


题目：
编写用Newmark β法（或者Wilsonθ法）计算结构振动响应的程序

谢谢。

[zsy312]

A,B,C的维数开始为1X4，第二次为1X8，此时u1，u2维数还是4X1，矩阵相乘不能衔接，不出错就怪了。

[zsy312]

再者，第34行你的t怎么突然变成行向量了，加上第23行，PN的行列不对应，第一行是4列，其他行都是一列0，又是错误。