求各位帮忙看下求频响函数的程序！
m=220;k=100000;af=0.5;bt=0.02; %比例阻尼系数
M=; %质量矩阵
K=; %刚度矩阵
=eig(inv(M)*K); %求特征向量和特征值
w=sqrt(w);
fn=diag(w)/(2*pi); %无阻尼固有频率
C=af*M+bt*K; %阻尼矩阵
zeta=(v'*M*v)\(v'*C*v)/2/w; %阻尼比
zeta=diag(zeta); w=diag(w);
wd=sqrt(diag(eye(3))-zeta.*zeta).*w; %有阻尼固有频率
fnd=wd/(2*pi);
v=v./v(1,1); %按v(1,1)归一化
Mr=diag(v'*M*v); %模态质量
Kr=diag(v'*K*v); %模态刚度
Cr=diag(v'*C*v); %模态阻尼
i=1;
for k=1:1:3
for wx=0.1:0.1:100
R(i,k)=(1-(w(k)./wx).^2)/(Mr(k).*((1-(w(k)./wx).^2).^2+(2*(zeta(k).*w(k))./wx).^2));
I(i,k)=(2*(zeta(k).*w(k))./wx)./(Mr(k).*((1-(w(k)./wx).^2).^2+(2*(zeta(k).*w(k))./wx).^2));
Y(i,k)=R(i,k)+j.*I(i,k);
i=i+1;
end
i=1;
end
i=1:1:1000;
H1(i)=Y(i,1).*(v(1,1).^2)+Y(i,2).*(v(1,2).^2)+Y(i,3).*(v(1,3).^2);
H2(i)=Y(i,1).*(v(2,1).*v(1,1))+Y(i,2).*(v(2,2).*v(1,2))+Y(i,3).*(v(2,3).*v(1,3));
H3(i)=Y(i,1).*(v(3,1).*v(1,1))+Y(i,2).*(v(3,2).*v(1,2))+Y(i,3).*(v(3,3).*v(1,3));
wx=i/10;
subplot(211); plot(wx,abs(H1)); title('H11幅频特性曲线')
subplot(212); plot(wx,180*angle(H1)/pi); title('H11相频特性曲线')
figure,subplot(211); plot(wx,real(H1)); title('H11实频特性曲线')
subplot(212); plot(wx,imag(H1)); title('H11虚频特性曲线')
figure,plot(real(H1),imag(H1)); title('H11导纳圆Nyquist Circle')
要求画出第一列频响函数的加速度幅频、相频、实频、虚频、Nyquist圆
麻烦各位看看程序有什么问题，如何改进？
本人菜鸟，各位狠拍！
[ 本帖最后由 ChaChing 于 2009-10-16 21:39 编辑 ]
回复 楼主 niluohe 的帖子
报错??? Complex values cannot be converted to logicals.
回复 沙发 ch_j1985 的帖子
没报错!? 那一行?To 楼主 : 明早有空再看看!
[ 本帖最后由 ChaChing 于 2008-12-19 21:38 编辑 ]
回复 板凳 ChaChing 的帖子
第一次试的时候报错，这一行H1(i)=Y(i,1).*(v(1,1).^2)+Y(i,2).*(v(1,2).^2)+Y(i,3).*(v(1,3).^2);不过现在又试了一下，确实没有报错！ 今天看了下楼主的程序, 不知该如何说呢?
没大略说明求解过程及公式, 真是考记忆! 有点懒得帮忙改! 仅说明与个人习惯不同处, 参考参考!
1.不需要点乘/除就不要用
2.式子尽量不要太复杂, 愈简化愈好 楼主看过5楼的回覆了吗? 可有醒思? 若估计没错, 楼主matlab初学者吧!
以下先假设楼主的式子没有用错, 我根据我的习惯重写,楼主参考参考!楼主可以自行检查是与你所列一样结果!
m=220; k=100000; af=0.5; bt=0.02;
M=; K=; C=af*M+bt*K;
=eig(M\K); v=v./v(1,1);Mr=diag(v'*M*v); Kr=diag(v'*K*v); Cr=diag(v'*C*v);
wn=sqrt(diag(w)); fn=wn/(2*pi); zeta=Cr./(2*Mr.*wn); wd=sqrt(1-zeta.*zeta).*wn; fd=wd/(2*pi);
wx=0.1:0.1:100;
for k=1:3, for i=1:length(wx),
rr=wn(k)/wx(i); a=1-rr^2; b=2*zeta(k)*rr;
R(i,k)=1/Mr(k)*a/(a^2+b^2); I(i,k)=1/Mr(k)*b/(a^2+b^2); Y(i,k)=R(i,k)+j*I(i,k);
end; end
HH=Y*(v.*repmat(v(1,:),3,1))'; H1=HH(:,1);
figure,subplot(211); plot(wx,abs(H1)); title('H11 Mag'); subplot(212); plot(wx,180*angle(H1)/pi); title('H11 Phase')
figure,subplot(211); plot(wx,real(H1)); title('H11 Real'); subplot(212); plot(wx,imag(H1)); title('H11 Imag')
figure; plot(real(H1),imag(H1)); title('H11 Nyquist Plot')
[ 本帖最后由 ChaChing 于 2008-12-21 14:39 编辑 ] 为何楼上假设楼主的式子没有用错? 除了与matlab使用无关外, 尚因为楼主并未给出求解过程及使用公式, 所以不确定对错! 但总觉得楼主好像误用了!
1.频响函数式子中之频率比率应该是rr=w/wn!?
2.频响函数式子中之系数应该是1/Kr并非1/Mr!?
3.频响函数式子1/(a+j*b)之实虚部分R & I计算有误!?
4.求H1/H2/H3的式子好像也不对!?
楼主检查看看吧! 希望楼主早日完成! 不会太过狠拍吧!
页:
[1]