clear
k=[ 120 -120    0    0    0;        #剛性マトリクス
   -120  270 -150    0    0;
      0 -150  330 -180    0;
      0    0 -180  390 -210;
      0    0    0 -210  420];

m=[0.055     0     0     0     0;   #質量マトリクス
       0 0.045     0     0     0;
       0     0 0.045     0     0;
       0     0     0 0.045     0;
       0     0     0     0 0.050];

r=inv(m)*k;                              # [k]{u]=p2[m][u]
[u,p2]=eig(r);                           # u−固有ベクトル　p2−固有値
for i=1:5; w2(i)=p2(i,i); endfor         # p2行列をw2ベクトルに代入
for i=1:5; tt(i)=2*pi/sqrt(p2(i,i)); endfor
[s,i]=sort(w2);                          # w2のソート
u=u(:,i);                                # iインデックスにより固有ベクトルの並べ換え
beta=inv(u'*m*u)*u'*m*[1 1 1 1 1]';      # beta−刺激係数
for i=1:5;b(i,:)=beta'.*u(i,:); endfor   # [βｕ] 刺激関数ベクトル
b
bb=zeros(6,5);
for i=1:5;bb(i+1,:)=b(6-i,:); endfor     # 各モード形の表示
title("各次モード［βＵ］"); ylabel("(階)"); grid "on"
plot(bb(:,1),0:5, bb(:,2),0:5, bb(:,3),0:5, bb(:,4),0:5, bb(:,5),0:5)