chnyyong
2010-05-12, 21:29
clear;
close all;
clc;
syms m1 m2 Iz1 Iz2 a1 a2 b1 b2 d C k1 k2 k3 s V
Vcount=320000:20000:400000;
A=[V*(m1+m2)*s+k1+k2+k3,-m2*d*s+(k1*a1-k2*b1-k3*d)/V+V*m1,-a2*m2*s-k3*(b2+a2)/V+V*m2,V*m2*s+k3;m1*V*d*s+k1*(a1+d)-k2*(b1-d),Iz1*s+[k1*a1*(a1+d)+k2*b1*(b1-d)]/V+m1*V*d,0,C*s;m1*a2*V*s+(k1+k2)*a2-k3*b2,[k3*b2*d+(k1*a1-k2*b1)*a2]/V+m1*V*a2,Iz2*s+k3*(b2+a2)*b2/V,-C*s-k3*b2;0,-1,1,s]
%上面的A就是传递函数的分母
H2=subs(A,{ m1 m2 Iz1 Iz2 a1 a2 b1 b2 d C k1 k3 V },{ 14170 11830 123290 55620 2.635 4.443 2.965 0.317 2.565 50000 292000 500000 20 });%把参数值（见罗论文中的附表）代入A中
j=length(Vcount)
ADET=det(H2)%求矩阵的行列式，此处即是求特征方程
for i=1:j
K=320000+20000*(i-1)
h=subs(ADET,k2 ,K);%把k2的值代入特征方程
N=sym2poly(h);%求特征方程的系数向量
E=roots(N)%解特征方程的根，即要求的特征根
t=find(real(E)==max(real(E))&imag(E)>=0); %特征根是成对出现且与x轴对称，因此只找出虚部大于等于零这个根即可
R=E(t);
All(i)=R;
end
RealP=real(All);%求特征根的实部
ImagP=imag(All);%求特征根的虚部
%subplot(2,1,2)
plot(RealP,ImagP,'b-*','LineWidth',1.5)%画特征根图，即根轨迹图
set(gca,'XDir','reverse')
xlabel('实轴/rad?s^{-1}');
ylabel('虚轴/rad?s^{-1}');
%title('k2变化，k1、k3不变化时的根轨迹图');
%axis([-8,0,0,7])
box off
for i=1:j
K=320000+20000*(i-1)
%NUM=1:length(RealP);
text(RealP(i),ImagP(i),num2str(K))%在数据点处标出速度的值
end
这是我毕业设计的一个程序，是求k1的变化曲线，现在想把k2的变化曲线跟k1的放在同一个坐标中，应该怎么样写程序，请高手帮下忙，谢谢！
close all;
clc;
syms m1 m2 Iz1 Iz2 a1 a2 b1 b2 d C k1 k2 k3 s V
Vcount=320000:20000:400000;
A=[V*(m1+m2)*s+k1+k2+k3,-m2*d*s+(k1*a1-k2*b1-k3*d)/V+V*m1,-a2*m2*s-k3*(b2+a2)/V+V*m2,V*m2*s+k3;m1*V*d*s+k1*(a1+d)-k2*(b1-d),Iz1*s+[k1*a1*(a1+d)+k2*b1*(b1-d)]/V+m1*V*d,0,C*s;m1*a2*V*s+(k1+k2)*a2-k3*b2,[k3*b2*d+(k1*a1-k2*b1)*a2]/V+m1*V*a2,Iz2*s+k3*(b2+a2)*b2/V,-C*s-k3*b2;0,-1,1,s]
%上面的A就是传递函数的分母
H2=subs(A,{ m1 m2 Iz1 Iz2 a1 a2 b1 b2 d C k1 k3 V },{ 14170 11830 123290 55620 2.635 4.443 2.965 0.317 2.565 50000 292000 500000 20 });%把参数值（见罗论文中的附表）代入A中
j=length(Vcount)
ADET=det(H2)%求矩阵的行列式，此处即是求特征方程
for i=1:j
K=320000+20000*(i-1)
h=subs(ADET,k2 ,K);%把k2的值代入特征方程
N=sym2poly(h);%求特征方程的系数向量
E=roots(N)%解特征方程的根，即要求的特征根
t=find(real(E)==max(real(E))&imag(E)>=0); %特征根是成对出现且与x轴对称，因此只找出虚部大于等于零这个根即可
R=E(t);
All(i)=R;
end
RealP=real(All);%求特征根的实部
ImagP=imag(All);%求特征根的虚部
%subplot(2,1,2)
plot(RealP,ImagP,'b-*','LineWidth',1.5)%画特征根图，即根轨迹图
set(gca,'XDir','reverse')
xlabel('实轴/rad?s^{-1}');
ylabel('虚轴/rad?s^{-1}');
%title('k2变化，k1、k3不变化时的根轨迹图');
%axis([-8,0,0,7])
box off
for i=1:j
K=320000+20000*(i-1)
%NUM=1:length(RealP);
text(RealP(i),ImagP(i),num2str(K))%在数据点处标出速度的值
end
这是我毕业设计的一个程序，是求k1的变化曲线，现在想把k2的变化曲线跟k1的放在同一个坐标中，应该怎么样写程序，请高手帮下忙，谢谢！