newmark－beta法计算多自由度系统响应

alexxc

编了一个用newmark－beta法计算荷载作用下的多自由度系统响应的函数，可是运行时出现以下错误：
>> vtb5(0.1,100);
??? Error: File: e:\MATLAB7\work\vtb6.m Line: 16 Column: 9
The expression to the left of the equals sign is not a valid target for an assignment.
不知哪里出了问题？请教各位高手帮忙

function vtb5 (deltt,tf)
%用newmark－beta法计算荷载作用下的多自由度系统的响应
fid1=fopen('disp','wt');%建立一个位移文件disp.dat
m=2*[1 0 0;0 1 0;0 0 1];
c=1.5*[2 -1 0;-1 2 -1;0 -1 2];
k=50*[2 -1 0;-1 2 -1;0 -1 2];
x0=[1 1 1]';
v0=[1 1 1]';
beta=1/6;
md=inv(m+deltt/2*c+beta*deltt^2*k);
m1=inv(m);
[E,F]=eig(m1*k);
diag(sqrt(F));%计算固有频率
for t=0:deltt:tf;
    f=[2*sin(3.754*t) -2*cos(2.2*t) 1.0*sin(2.8*t)]';
    if t=0; %(line 16)        xdd0=m1*(f-k*x0-c*v0);%计算初始加速度
        
    else
        xdd=md*(f-c*(v0+deltt/2*xdd0)-k*(x0+deltt*v0+(1/2-beta)*deltt^2*xdd0));%计算加速度
        xd=v0+deltt/2*(xdd0+xdd);%计算速度
        x=x0+deltt*v0+deltt^2/2*xdd0+beta*deltt^3*(xdd-xdd0)/deltt;%计算位移
        xdd0=xdd;v0=xd;x0=x;
        fprintf(fid1,'%10.4f',x0);%向文件中写数据
        t %显示计算时间步长
    end
end   
    fid2=fopen('disp','rt');%打开disp.dat文件
    n=tf/deltt;%disp.dat文件中位移的个数
    x=fscanf(fid2,'%f',[3,n]);%将disp.dat文件中的位移写成矩阵
    t=1:n;
    figure('numbertitle','off','name','point1 displacement','pos',[450 180 400 420]);
    plot(t,x(1,t)),grid,xlabel('time*0.1s'),title('点1的位移与时间的关系');
    figure('numbertitle','off','name','point2 displacement','pos',[350 160 400 420]);
    plot(t,x(2,t)),grid,xlabel('time*0.1s'),title('点2的位移与时间的关系');
    figure('numbertitle','off','name','point3 displacement','pos',[250 140 400 420]);
    plot(t,x(3,t)),grid,xlabel('time*0.1s'),title('点3的位移与时间的关系');

xjzuo

%%%逻辑运算符号不对, 可改为:
%%%%%%
if t==0;
%%%%%%

alexxc

谢谢楼上这位兄弟的回复，犯了个低级错误。。

only8810205028

请问若是更多自由度，且质量矩阵不为对角阵
我将楼主的程序稍微变了下，可是得出了这样的结果

1. k=sysK;% K ----- 调用刚度矩阵
   
2. m=sysM; % M ----- 调用质量矩阵
   
3. c=zeros(51,51);
   
4. x0=zeros(51,1);
   
5. x0(1:51,1)=1;
   
6. v0=zeros(51,1);
   
7. v0(1:51,1)=1;
   
8. f=zeros(51,1);
   
9. 
   
10. 
    
11. m1=inv(m);
    
12. [E,F]=eig(m1*k);
    
13. fn=diag(sqrt(F));%计算固有频率
    
14. 
    
15. beta=1/6;
    
16. tf=20;
    
17. deltt=0.01;
    
18. md=inv(m+deltt/2*c+beta*deltt^2*k);
    
19. for t=0:deltt:tf;
    
20.     f(1:51,1)=2*sin(3.754*t);
    
21.     %    f=[2*sin(3.754*t) -2*cos(2.2*t) 1.0*sin(2.8*t)]';
    
22.     if t==0; %(line 16)      
    
23.       xdd0=m1*(f-k*x0-c*v0);%计算初始加速度   
    
24.     else
    
25.         xdd=md*(f-c*(v0+deltt/2*xdd0)-k*(x0+deltt*v0+(1/2-beta)*deltt^2*xdd0));%计算加速度
    
26.         xd=v0+deltt/2*(xdd0+xdd);%计算速度
    
27.         x=x0+deltt*v0+deltt^2/2*xdd0+beta*deltt^3*(xdd-xdd0)/deltt;%计算位移
    
28.         xdd0=xdd;v0=xd;x0=x;
    
29.         fprintf(fid1,'%10.4f',x0);%向文件中写数据
    
30.         t %显示计算时间步长
    
31.     end
    
32. end   
    
33. 
    
34. 
    
35.     fid2=fopen('disp','rt');%打开disp.dat文件
    
36.     n=tf/deltt;%disp.dat文件中位移的个数
    
37.     x=fscanf(fid2,'%f',[3,n]);%将disp.dat文件中的位移写成矩阵
    
38.     t=1:n;
    
39.     figure('numbertitle','off','name','point1 displacement','pos',[450 180 400 420]);
    
40.     plot(t,x(1,t)),grid,xlabel('time*0.1s'),title('点1的位移与时间的关系');
    
41.     figure('numbertitle','off','name','point2 displacement','pos',[350 160 400 420]);
    
42.     plot(t,x(2,t)),grid,xlabel('time*0.1s'),title('点2的位移与时间的关系');
    
43.     figure('numbertitle','off','name','point3 displacement','pos',[250 140 400 420]);
    
44.     plot(t,x(3,t)),grid,xlabel('time*0.1s'),title('点3的位移与时间的关系');
    

复制代码

望高手解答！不胜感激！