频响函数计算程序

shock

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圆
麻烦各位看看程序有什么问题，如何改进？