请高手帮忙看一下MATLAB程序错在哪儿|地震工程『抗震』 - 领先的岩土技术社区，岩土领域的媒体、社区与应用平台！

离线ewarren

副总工

发帖: 556

土币: 2883

威望: 4403

原创币: 0

只看楼主倒序阅读使用道具楼主发表于: 2010-05-15

程序如下，运行时提示：“Matrix mast be square”，上网查了一下据说是要把“*”改为“.*”就可以了，可是我试了一下，貌似还不行，请高手帮帮忙啊

复制代码

for i=1:1:1000
for k=0:1:1000
j=2*pi; %θ=2π
w=(2*pi)/(i/100); %自振圆频率ω
t=k/100; %时间
u(1)=0; %初始位移为0
u1(1)=0; %初始速度为0
c11=exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01+(g/sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01;
c12=exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01;
d11=(1/w^2)*(1+3*g/w*0.01+(4*g^2-1)/w^2*0.01^2)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01+(g+3*g/w*0.01+g(4*g^2-1)/w^2*0.01^2-3/2*w*0.01-2*g/w^2*0.01^2)*(1/w^2*sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01-(1/w^2)*(g/w*0.01+(4*g^2-1)/w^2*0.01^2);
d12=-(2/w^3*0.01)*(2*g+(4*g^2-1)/w*0.01)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01-(2*g^2+(4*g^2-1)*g/w*0.01-1-2*g/w*0.01)*(2/0.01*w^3*sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01+(1/w^2*0.01)*(2*(4*g^2-1)/w^2*0.01-0.01);
d13=(1/w^3*0.01)*(g+(4*g^2-1)/w^2*0.01)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01+(g^2+(4*g^2-1)*g/w*0.01-1/2-2*g/w*0.01)*(1/0.01*w^3*sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01-((4*g^2-1)/w^2*0.01-g/w)*(1/w^2*0.01);
c21=-(w/sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01;
c22=-(g/sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01+exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01;
d21=-(3+4*g/w*0.01)*(1/2*w^2*0.01)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01-(1+3*g/2*w*0.01+(2*g^2-1)/w^2*0.01^2)*(1/w*sqrt(1-g^2));
d22=(1+2*g/w*0.01)*(2/w^2*0.01)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01+(g+(2*g^2-1)/w*0.01)*(2/0.01*w^2*sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01-4*g/w^3*0.01^2;
d23=-(1+4*g/w*0.01)*(1/2*w^2*0.01)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01-(g/2+(2*g^2-1)/w*0.01)*(1/0.01*w^2*sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01+(1/(2*w^2*0.01))*(4*g/w*0.01-1);
u(k+2)=c11*u(k)+c12*u1(k)+d11*cos(j*(k+1))+d12*cos(j*(k+2))+d13*cos(j*(k+3));
u1(k+2)=c21*u(k)+c22*u1(k)+d21*cos(j*(k+1))+d22*cos(j*(k+2))+d23*cos(j*(k+3));
u11(k+2)=-2*g*w*u1(k+1)-w^2*u(k+1);
end
wy=max(u);
sd=max(u1);
jsd=max(u11);
end
%-------输出结果-------
fid=fopen('抛物线内插法位移解.txt','w');
fprintf(fid,'%15.5f\r\n',wy);
fclose(fid);
fid=fopen('抛物线内插法速度解.txt','w');
fprintf(fid,'%15.5f\r\n',sd);
fclose(fid);
fid=fopen('抛物线内插法加速度解.txt','w');
fprintf(fid,'%15.5f\r\n',jsd);
fclose(fid);
%---------画图-------
figure;
T=0.01:0.01:10;
plot(T,wy);
title('y(t)');
xlabel('T');
ylabel(g);
grid on;
%--------画图-------
figure;
T=0.01:0.01:10;
plot(T,sd);
title('y′(t)');
xlabel('T');
ylabel(g);
grid on;
%--------画图--------
figure;
T=0.01:0.01:10;
plot(T,jsd);
title('y′′(t)');
xlabel('T');
ylabel(g);
grid on;