求各位帮忙看下求频响函数的程序！

niluohe

m=220;k=100000;
af=0.5;bt=0.02;         %比例阻尼系数
M=[m 0 0;0 m 0;0 0 m];         %质量矩阵
K=[2*k -k 0;-k 2*k -k;0 -k k];           %刚度矩阵
[v,w]=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 编辑 ]

ch_j1985

报错
??? Complex values cannot be converted to logicals.

ChaChing

没报错!? 那一行?

To 楼主 : 明早有空再看看!

[ 本帖最后由 ChaChing 于 2008-12-19 21:38 编辑 ]

ch_j1985

第一次试的时候报错，这一行H1(i)=Y(i,1).*(v(1,1).^2)+Y(i,2).*(v(1,2).^2)+Y(i,3).*(v(1,3).^2);
不过现在又试了一下，确实没有报错！

ChaChing

今天看了下楼主的程序, 不知该如何说呢?
没大略说明求解过程及公式, 真是考记忆! 有点懒得帮忙改! 仅说明与个人习惯不同处, 参考参考!
1.不需要点乘/除就不要用
2.式子尽量不要太复杂, 愈简化愈好

ChaChing

楼主看过5楼的回覆了吗? 可有醒思? 若估计没错, 楼主matlab初学者吧!
以下先假设楼主的式子没有用错, 我根据我的习惯重写,  楼主参考参考!  楼主可以自行检查是与你所列一样结果!

m=220; k=100000; af=0.5; bt=0.02;
M=[m 0 0;0 m 0;0 0 m]; K=[2*k -k 0;-k 2*k -k;0 -k k]; C=af*M+bt*K;
[v,w]=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 编辑 ]

ChaChing

为何楼上假设楼主的式子没有用错? 除了与matlab使用无关外, 尚因为楼主并未给出求解过程及使用公式, 所以不确定对错! 但总觉得楼主好像误用了!
1.频响函数式子中之频率比率应该是rr=w/wn!?
2.频响函数式子中之系数应该是1/Kr并非1/Mr!?
3.频响函数式子1/(a+j*b)之实虚部分R & I计算有误!?
4.求H1/H2/H3的式子好像也不对!?
楼主检查看看吧! 希望楼主早日完成! 不会太过狠拍吧!