如何求解有阻尼振动方程

david

对于质量MM3、刚度KK3、阻尼矩阵CC3均为对称矩阵的情况下，如何求解有阻尼振动方程的固有频率？在钟一谔的 《转子动力学》一书中，他介绍了通过用质量矩阵的逆矩阵、质量矩阵、刚度矩阵、以及单位矩阵、零矩阵而组合成一个新的矩阵MCKD（2n*2n）然后再对这个矩阵用QR算法来求解？进而求出特征值，但我对这个方法不了解，请各位大哥给予指点，谢谢!qq：363420100
看到有人这么运算，但不知道为什么，给讲解一下
MCKD=[-inv(MM3)*CC3 -inv(MM3)*KK3;eye(33) zeros(33,33)]
[QQ,RR]=qr(MCKD);%矩阵正交三角化分解，其中QQ为酉矩阵RR为上三角阵
[AD3,WD3]=eig(RR);%求特征值和特征向量

[ 本帖最后由 david 于 2008-4-17 10:18 编辑 ]

gh688

如果你考虑阻尼的话最好把方程划为状态空间的形式求解，如果你想省时间可以直接用eig，当然也可以自己编程。你的3个矩阵都对称所以应该比较好解(指的是解得精确性)。另外在计算中尽量不要出现求逆inv，你的QR分解也没有必要。

david

[V,D]=eig(K,M);
[DD,I]=sort(diag(sqrt(D))/(2*pi));%固有频率
DDD=DD*60;  %临界转速
我看到有人用上面这段程序来求解，但不明白为什么除以(2*pi），这是其一；其二是固有频率*60就是临界转速吗？

gh688

原帖由 david 于 2008-4-17 17:18 发表

                               
登录/注册后可看大图

[V,D]=eig(K,M);
[DD,I]=sort(diag(sqrt(D))/(2*pi));%固有频率
DDD=DD*60;  %临界转速
我看到有人用上面这段程序来求解，但不明白为什么除以(2*pi），这是其一；其二是固有频率*60就是临界转速吗？

1 　如果不考虑阻尼的话，你可以用上面的公式，但是eig(K,M)缺少负号，应该是eig(-K,M)
2　eig求出的特征向量的单位是rad/s，除2*pi变成频率单位HZ，频率*60变成转速单位rpm
3　eig求出的只是固有频率，这个量和转速有关，确切的说不是准确临界转速，这是两个概念．确切的临界转速你要参考campbell图形．

david

1、我想求考虑阻尼情况下的固有频率，如何求？
2、campbell图形是什么？如何获得？刚学matlab，请多多指导，谢谢

咕噜噜

有阻尼那就不叫固有频率了，那叫振动频率
很多振动书上都有，可以找找看
campbell图形？没见过

gh688

原帖由 咕噜噜 于 2008-4-17 18:00 发表

                               
登录/注册后可看大图

有阻尼那就不叫固有频率了，那叫振动频率
很多振动书上都有，可以找找看
campbell图形？没见过

我刚才误以为他就是求转子系统的临界转速，所以把转子动力学的坎贝图给搬弄出来了，呵呵

david

是旋转系统，水轮发电机主轴系统，

david

function Drawing_diagram=Drawing(M,K,C,numb)
w=75;          %主轴的转速
wmr=w*2*pi/60;            %转轮的角频率 角速度
M0=1047*10^3;             %转子质量
Ms0=286*10^3;         %水轮机转轮的质量
%-------------the force on water seal-------%
Qr=580;            %额定流量
Rt=4;              %转轮半径
gw=1.002e-3;            %水的粘度
ros=1e3;              %水的密度
gg=0.5;             %局部阻力系数
gl=0.1;             %沿程阻力系数
ls=0.47;            %密封长度
b=0.00286;               %水封平均间隙宽度
g0=0.92;             %流量损失系数
s0=2*10^(-3);             %偏心距
Q1=(1-g0)*Qr;             %有效流量
C1=Q1/(2*pi*Rt*b);             %密封中水流的速度
hou1=0.063079;            %进水叶片平均厚度
hou2=0.037872;           %出水叶片平均厚度
r1=3.691;            %进水边曲线AB重心距轴线的距离，该值是约值
r2=1.920;             %出水边曲线A1B1重心距轴线的距离，该值是约值
B1=3.251;             %进水边曲线AB的长度，该值是约值
B2=5.328;             %出水边曲线A1B1的长度，该值是约值
c1=30;          %导水叶片出口角范围在20-40之间
c2=15;          %叶片出口安放角，见水力机械原理（印度）P125
Jw=1021149.781;              %水流等效转动惯量
Z=13;             %转轮叶片总数
%--------------------%
uu=1/sqrt(1+gg+(gl*ls/(2*b)));
Kw=(ros*(1+gg+gl*ls/(2*b))*C1^2/2)*Rt*(ls^2)*pi*gl*uu^4*(s0/b)/(2*b);            %水封处水力不平衡力鲁式3-3
K(13,13)=K(13,13)+Kw;
K(27,27)=K(27,27)+Kw;
%-----------hydraulic resistance-------%
cfa=6*pi*gw*(1+Rt/sqrt(2*gw/(wmr*ros)));%？？？？？？？？？？？？？？？
mfa=6*pi*Rt*sqrt(gw*ros/(2*wmr))*(1+Rt/(3*sqrt(2*gw/wmr*ros)));%？？？？？？？？？？？？
M(13,13)=M(13,13)+mfa;
M(27,27)=M(27,27)+mfa;
C(13,13)=C(13,13)+cfa;
C(27,27)=C(27,27)+cfa;
e0=1.2e-2;               %转子质量偏心距
es0=2e-2;             %转轮质量偏心距

%---------Newmark----%
delt=0.008;             % time step size
r=0.5;                  % parameters
alfa=0.25*(0.5+r)^2;              %parameters
a0=1/(alfa*(delt)^2);           %a0-a7积分常数
a1=r/(alfa*delt);
a2=1/(alfa*delt);
a3=1/(2*alfa)-1;
a4=r/alfa-1;
a5=delt*(r/alfa-2)/2;
a6=delt*(1-r);
a7=r*delt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------------------------calaculate temproal response---------------%
u0=zeros(33,1);            %initial displacement
u0(3)=1e-3;               % initial displacement of rotor in y-axis赋值
u0(13)=1e-3;             % initial displacement of shuilunji in y-axis赋值
ud0=zeros(33,1);            %initial velocity
ud0(33)=75*2*pi/60;            %initial rotational velocity
ud0(17)=4.5e-3;                 %initial velocity of rotor in z-axis赋值
ud0(27)=3.8e-3;                %initial velocity of shuilunji in z-axis赋值
Mw=ros*Q1*((Q1*cos(c1)*r1/(2*r1*pi*B1-hou1*B1*Z))-r2*(pi*w*r2/30-cot(c2)*Q1/(2*r2*pi*B2-hou2*B2*Z)))+Jw*ud0(33);%水力矩
f=zeros(33,1);
udd0=inv(M)*(f-K*u0-C*ud0);           %initial acceleration起始加速度
udd0=zeros(33,1);
Kef=K+a0*M+a1*C;             %effective stiffness matrix有效刚度矩阵
i=1;j=1;                  % parameter for controling loop
delt=0.002
for t=0:delt:1  
     f(3,1)=M0*e0*wmr^2*(cos(wmr*t)-sin(wmr*t));%转子处的机械不平衡力
     f(17,1)=M0*e0*wmr^2*(sin(wmr*t)+cos(wmr*t));%转子处的机械不平衡力
     f(13,1)=Ms0*es0*wmr^2*(cos(wmr*t)-sin(wmr*t));%转轮处的机械不平衡力
     f(27,1)=Ms0*es0*wmr^2*(sin(wmr*t)+cos(wmr*t));%转子处的机械不平衡力
     f(33,1)=f(33,1)+Mw*sin(wmr*t);
     f1=f+M*(a0*u0+a2*ud0+a3*udd0)+C*(a1*u0+a4*ud0+a5*udd0);          %effective loads有效载荷
    %%-----------------%
    u=Kef\f1;          %一个时间间隔的位移
    ddu=a0*(u-u0)-a2*ud0-a3*udd0;           %一个时间间隔的加速度
    du=ud0+a6*udd0+a7*ddu;           %一个时间间隔的速度
    u0=u;
    ud0=du;
    udd0=ddu;
    if (t==0)|(t<=1)   
        temp_ttt(j)=t;
        temp_f1_n2(j)=sqrt(f1(3)^2+f1(17)^2);
        temp_f1_n7(j)=sqrt(f1(13)^2+f1(27)^2);
        
        %temp_u1y(j)=u(1)+rand(1)*0.000009*sin(rand(1)*100);temp_du1y(j)=du(1);temp_ddu1y(j)=ddu(1);
        %temp_u1z(j)=u(15)+rand(1)*0.000009*sin(rand(1)*100);temp_du1z(j)=du(15);temp_ddu1z(j)=ddu(15);
        %temp_u1y(j)=u(1)+unifrnd(-0.00009,0.00009,[i,1]);temp_du1y(j)=du(1);temp_ddu1y(j)=ddu(1);
        %temp_u1z(j)=u(15)+unifrnd(-0.00009,0.00009,[i,1]);temp_du1z(j)=du(15);temp_ddu1z(j)=ddu(15);
        temp_u1y(j)=u(1);temp_du1y(j)=du(1);temp_ddu1y(j)=ddu(1);
        temp_u1z(j)=u(15);temp_du1z(j)=du(15);temp_ddu1z(j)=ddu(15);     
        
        temp_u2y(j)=u(3);temp_du2y(j)=du(3);temp_ddu2y(j)=ddu(3);
        temp_u2z(j)=u(17);temp_du2z(j)=du(17);temp_ddu2z(j)=ddu(17);
        
        temp_u3y(j)=u(5);temp_du3y(j)=du(5);temp_ddu3y(j)=ddu(5);
        temp_u3z(j)=u(19);temp_du3z(j)=du(19);temp_ddu3z(j)=ddu(19);
        temp_u3x(j)=u(29);
        
        temp_u4y(j)=u(7);temp_du4y(j)=du(7);temp_ddu4y(j)=ddu(7);
        temp_u4z(j)=u(21);temp_du4z(j)=du(21);temp_ddu4z(j)=ddu(21);
        temp_u4x(j)=u(30);
        
        temp_u5y(j)=u(9);temp_du5y(j)=du(9);temp_ddu5y(j)=ddu(9);
        temp_u5z(j)=u(23);temp_du5z(j)=du(23);temp_ddu5z(j)=ddu(23);
        temp_u5x(j)=u(31);
        
        temp_u6y(j)=u(11);temp_du6y(j)=du(11);temp_ddu6y(j)=ddu(11);
        temp_u6z(j)=u(25);temp_du6z(j)=du(25);temp_ddu6z(j)=ddu(25);
        temp_u6x(j)=u(32);
        
        temp_u7y(j)=u(13);temp_du7y(j)=du(13);temp_ddu7y(j)=ddu(13);
        temp_u7z(j)=u(27);temp_du7z(j)=du(27);temp_ddu7z(j)=ddu(27);
        temp_u7x(j)=u(33);
        
        j=j+1;
    end
end
temp_f1_n2=temp_f1_n2';
temp_f1_n7=temp_f1_n7';
temp_ttt=temp_ttt';
temp_u1y=temp_u1y';
temp_u1z=temp_u1z';
temp_u2y=temp_u2y';
temp_u2z=temp_u2z';
temp_u3y=temp_u3y';
temp_u3z=temp_u3z';
temp_u3x=temp_u3x';
temp_u4y=temp_u4y';
temp_u4z=temp_u4z';
temp_u4x=temp_u4x';
temp_u5y=temp_u5y';
temp_u5z=temp_u5z';
temp_u5x=temp_u5x';
temp_u6y=temp_u6y';
temp_u6z=temp_u6z';
temp_u6x=temp_u6x';
temp_u7y=temp_u7y';
temp_u7z=temp_u7z';
temp_u7x=temp_u7x';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----------calaculate dynamic response--%
u0=zeros(33,1); %initial displacement
ud0=zeros(33,1); %initial velocity
f=zeros(33,1);
udd0=inv(M)*(f-K*u0-C*ud0); %initial acceleration
Kef=K+a0*M+a1*C;  %effective stiffness matrix
ii=1;j=1; % parameter for controling loop
for t=0:delt:2
%for t=0:delt:10  
    %     if (t==5)
    %         MMs0=Ms0-1e3;
    %         Ms0=MMs0;
    %     end
       %------------------centrifugal force--%
     f(3,1)=M0*e0*wmr^2*(cos(wmr*t)-sin(wmr*t));%转子处的机械不平衡力
     f(17,1)=M0*e0*wmr^2*(sin(wmr*t)+cos(wmr*t));%转子处的机械不平衡力
     f(13,1)=Ms0*es0*wmr^2*(cos(wmr*t)-sin(wmr*t));%转轮处的机械不平衡力
     f(27,1)=Ms0*es0*wmr^2*(sin(wmr*t)+cos(wmr*t));%转子处的机械不平衡力
     f(33,1)=f(33,1)+Mw;

    %---------assuming there is an exciting force acting on the runner suddenly----------%
    if (t>0)&(t<=2)
   
        f(13,1)=f(13,1)+Ms0*0.01;
        f(27,1)=f(27,1)+Ms0*0.02;
        f(33,1)=f(33,1)+Mw*0.3;
    end
    f2=f+M*(a0*u0+a2*ud0+a3*udd0)+C*(a1*u0+a4*ud0+a5*udd0); %effective loads
    u=Kef\f2;
    ddu=a0*(u-u0)-a2*ud0-a3*udd0;
    du=ud0+a6*udd0+a7*ddu;
    u0=u;
    ud0=du;
    udd0=ddu;
    if (t==0)|(t<=2)
        
        ttt(j)=t;
        %f2_n2(j)=sqrt(f1(3)^2+f1(17)^2);
        %f2_n7(j)=sqrt(f1(13)^2+f1(27)^2);
        f2_n2(j)=sqrt(f2(3)^2+f2(17)^2);
        f2_n7(j)=sqrt(f2(13)^2+f2(27)^2);
        u1y(j)=u(1);du1y(j)=du(1);ddu1y(j)=ddu(1);
        u1z(j)=u(15);du1z(j)=du(15);ddu1z(j)=ddu(15);
      
        u2y(j)=u(3);du2y(j)=du(3);ddu2y(j)=ddu(3);
        u2z(j)=u(17);du2z(j)=du(17);ddu2z(j)=ddu(17);
        
        u3y(j)=u(5);du3y(j)=du(5);ddu3y(j)=ddu(5);
        u3z(j)=u(19);du3z(j)=du(19);ddu3z(j)=ddu(19);
         u3x(j)=u(29);
      
        u4y(j)=u(7);du4y(j)=du(7);ddu4y(j)=ddu(7);
        u4z(j)=u(21);du4z(j)=du(21);ddu4z(j)=ddu(21);
         u4x(j)=u(30);
        
        u5y(j)=u(9);du5y(j)=du(9);ddu5y(j)=ddu(9);
        u5z(j)=u(23);du5z(j)=du(23);ddu5z(j)=ddu(23);
         u5x(j)=u(31);
      
        u6y(j)=u(11);du6y(j)=du(11);ddu6y(j)=ddu(11);
        u6z(j)=u(25);du6z(j)=du(25);ddu6z(j)=ddu(25);
         u6x(j)=u(32);
   
        u7y(j)=u(13);du7y(j)=du(13);ddu7y(j)=ddu(13);
        u7z(j)=u(27);du7z(j)=du(27);ddu7z(j)=ddu(27);
   
         u7x(j)=u(33);
        j=j+1;
    end
end
ttt=ttt';
f2_n2=f2_n2';
f2_n7=f2_n7';
u1y=u1y';
u1z=u1z';
u2y=u2y';
u2z=u2z';
u3y=u3y';
u3z=u3z';
u3x=u3x';
u4y=u4y';
u4z=u4z';
u4x=u4x';
u5y=u5y';
u5z=u5z';
u5x=u5x';
u6y=u6y';
u6z=u6z';
u6x=u6x';
u7y=u7y';
u7z=u7z';
u7x=u7x';
size_ttt=size(ttt);
size_temp_ttt=size(temp_ttt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------------------各个节点的瞬态响应和稳态不平衡响应图--------------------%
position=1;
%--------------

figure(numb);%创建图形窗口
subplot(1,2,1);%创建子图图区
plot(temp_u1y*1000,temp_u1z*1000);%创建二维线形图,以temp_u1y的元素为横坐标，temp_u1z中同维的元素为纵坐标
hold on;
plot(temp_u1y(1,1)*1000,temp_u1z(1,1)*1000,'rd');
plot(temp_u1y(size_temp_ttt(1,1),1)*1000,temp_u1z(size_temp_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点1瞬态响应');
axis square;
grid off;%网格显示控制
box on;%坐标轴外框显示

subplot(1,2,2);
plot(temp_u2y*1000,temp_u2z*1000);
hold on;
plot(temp_u2y(1,1)*1000,temp_u2z(1,1)*1000,'rd');
plot(temp_u2y(size_temp_ttt(1,1),1)*1000,temp_u2z(size_temp_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点2瞬态响应');
axis square;
grid off;box on;

figure(numb+1);
subplot(1,2,1);
plot(temp_u3y*1000,temp_u3z*1000);
hold on;
plot(temp_u3y(1,1)*1000,temp_u3z(1,1)*1000,'rd');
plot(temp_u3y(size_temp_ttt(1,1),1)*1000,temp_u3z(size_temp_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点3瞬态响应');
axis square;
grid off;box on;

subplot(1,2,2);
plot(temp_u4y*1000,temp_u4z*1000);
hold on;
plot(temp_u4y(1,1)*1000,temp_u4z(1,1)*1000,'rd');
plot(temp_u4y(size_temp_ttt(1,1),1)*1000,temp_u4z(size_temp_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点4瞬态响应');
axis square;
grid off;box on;

figure(numb+2);
subplot(1,2,1);
plot(temp_u5y*1000,temp_u5z*1000);
hold on;
plot(temp_u5y(1,1)*1000,temp_u5z(1,1)*1000,'rd');
plot(temp_u5y(size_temp_ttt(1,1),1)*1000,temp_u5z(size_temp_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点5瞬态响应');
axis square;
grid off;box on;

subplot(1,2,2);
plot(temp_u6y*1000,temp_u6z*1000);
hold on;
plot(temp_u6y(1,1)*1000,temp_u6z(1,1)*1000,'rd');
plot(temp_u6y(size_temp_ttt(1,1),1)*1000,temp_u6z(size_temp_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点6瞬态响应')
axis square;
grid off;box on;

figure(numb+3);
subplot(1,2,1);
plot(temp_u7y*1000,temp_u7z*1000);
hold on;
plot(temp_u7y(1,1)*1000,temp_u7z(1,1)*1000,'rd');
plot(temp_u7y(size_temp_ttt(1,1),1)*1000,temp_u7z(size_temp_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点7瞬态响应');
axis square;
grid off;box on;

figure(numb+4);
subplot(1,2,1);
plot(u1y*1000,u1z*1000);
hold on;
plot(u1y(1,1),u1z(1,1),'rd');
plot(u1y(size_ttt(1,1),1)*1000,u1z(size_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点1动态响应')
axis square;
grid off;box on;

subplot(1,2,2);
plot(u2y*1000,u2z*1000);
hold on;
plot(u2y(1,1)*1000,u2z(1,1)*1000,'rd');
plot(u2y(size_ttt(1,1),1)*1000,u2z(size_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点2动态响应');
axis square;
grid off;box on;

figure(numb+5);
subplot(1,2,1);
plot(u3y*1000,u3z*1000);
hold on;
plot(u3y(1,1)*1000,u3z(1,1)*1000,'rd');
plot(u3y(size_ttt(1,1),1)*1000,u3z(size_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点3动态响应');
axis square;
grid off;box on;

subplot(1,2,2);
plot(u4y*1000,u4z*1000);
hold on;
plot(u4y(1,1)*1000,u4z(1,1)*1000,'rd');
plot(u4y(size_ttt(1,1),1)*1000,u4z(size_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点4动态响应');
axis square;
grid off;box on;

figure(numb+6);
subplot(1,2,1);
plot(u5y*1000,u5z*1000);
hold on;
plot(u5y(1,1)*1000,u5z(1,1)*1000,'rd');
plot(u5y(size_ttt(1,1),1)*1000,u5z(size_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点5动态响应');
axis square;
grid off;box on;

subplot(1,2,2);
plot(u6y*1000,u6z*1000);
hold on;
plot(u6y(1,1)*1000,u6z(1,1)*1000,'rd');
plot(u6y(size_ttt(1,1),1)*1000,u6z(size_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点6动态响应');
axis square;
grid off;
box on;

figure(numb+7);
subplot(1,2,1);
plot(u7y*1000,u7z*1000);
hold on;
plot(u7y(1,1)*1000,u7z(1,1)*1000,'rd');
plot(u7y(size_ttt(1,1),1)*1000,u7z(size_ttt(1,1),1)*1000,'r*');
xlabel('y方向位移 /mm');
ylabel('z方向位移 /mm');
title('节点7动态响应');
axis square;
grid off;box on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  %----------------------------各个节点的时域图-----------------------------%
  %-------------节点1------------%
figure(numb+8);
%subplot(1,2,1);
plot(temp_ttt,sqrt(temp_u1y.^2+temp_u1z.^2));
%plot(temp_ttt,sqrt(temp_u1y.^2+temp_u1z.^2).*sin(wmr.*temp_ttt));
xlabel('时间 /s');
ylabel('幅值 /m');
title('节点1的时域图');
axis square;
grid off;  
box on;
%----------------节点1的频域图---------------%
%subplot(1,2,2);
%fs=80;
%nq=100;
%y=fft(sqrt(temp_u1y.^2+temp_u1z.^2),nq);%进行fft变换
%mag=abs(y);%求幅值
%f=(0:length(y)-1)'*fs/length(y);%进行对应的频率转换
%plot(f,mag);%做频谱图
%axis([0,1,0,1]);
%xlabel('频率(Hz)');
%ylabel('幅值');
%title('正弦信号y=2*pi*10t幅频谱图N=128');
%grid;

%-------------节点2------------%
figure(numb+9);
plot(temp_ttt,sqrt(temp_u2y.^2+temp_u2z.^2));
%plot(temp_ttt,sqrt(temp_u2y.^2+temp_u2z.^2).*sin(wmr.*temp_ttt));
xlabel('时间 /s');
ylabel('幅值 /m');
title('节点2的时域图');
grid off;  
box on;
%-------------节点3------------%
figure(numb+10);
plot(temp_ttt,sqrt(temp_u3y.^2+temp_u3z.^2));
%plot(temp_ttt,sqrt(temp_u3y.^2+temp_u3z.^2).*sin(wmr.*temp_ttt));
xlabel('时间 /s');
ylabel('幅值 /m');
title('节点3的时域图');
grid off;  
box on;
  %-------------节点4------------%
figure(numb+11);
plot(temp_ttt,sqrt(temp_u4y.^2+temp_u4z.^2));
%plot(temp_ttt,sqrt(temp_u4y.^2+temp_u4z.^2).*sin(wmr.*temp_ttt));
xlabel('时间 /s');
ylabel('幅值 /m');
title('节点4的时域图');
grid off;  
box on;
  %-------------节点5------------%
figure(numb+12);
plot(temp_ttt,sqrt(temp_u5y.^2+temp_u5z.^2));
%plot(temp_ttt,sqrt(temp_u5y.^2+temp_u5z.^2).*sin(wmr.*temp_ttt));
xlabel('时间 /s');
ylabel('幅值 /m');
title('节点5的时域图');
grid off;  
box on;
  %-------------节点6------------%
figure(numb+13);
plot(temp_ttt,sqrt(temp_u6y.^2+temp_u6z.^2));
%plot(temp_ttt,sqrt(temp_u6y.^2+temp_u6z.^2).*sin(wmr.*temp_ttt));
xlabel('时间 /s');
ylabel('幅值 /m');
title('节点6的时域图');
grid off;  
box on;
  %-------------节点7------------%
figure(numb+14);

plot(temp_ttt,sqrt(temp_u7y.^2+temp_u7z.^2));
%plot(temp_ttt,sqrt(temp_u7y.^2+temp_u7z.^2).*sin(wmr.*temp_ttt));
xlabel('时间 /s');
ylabel('幅值 /m');
title('节点7的时域图');
grid off;  
box on;

gh688

你有问题最好精练的概括一下，让大家给你出个注意．你这个程序包含了好多部分，从响应到时域频域分析太多了，也不知道你想讨论的东西．

david

我仿真出来的时域图，好像是错误的，也想请诸位帮我分析一下程序，待会我把时域图贴上来大家给看看

[ 本帖最后由 david 于 2008-4-17 20:59 编辑 ]

咕噜噜

哈，是坎贝尔图啊，没反映过来

vibrater

原帖由 gh688 于 2008-4-17 16:56 发表

                               
登录/注册后可看大图

如果你考虑阻尼的话最好把方程划为状态空间的形式求解，如果你想省时间可以直接用eig，当然也可以自己编程。你的3个矩阵都对称所以应该比较好解(指的是解得精确性)。另外在计算中尽量不要出现求逆inv，你的QR分解也没 ...

为什莫尽量不要inv？

gh688

比较病态的矩阵求逆的误差比较大，可以尽量采取别的形式啊，你可以找几个矩阵实验以下看看，举个例子：如果解方程组Ax=b，直接x=inv(A)*b效果一般不如用高斯消去法好，具体的你看一下求解的步骤就明白了