| Nums of M | 質量の数（自由度の数）（１〜１０の整数） |
| Nums of S | ばねの数（１〜１０の整数） |
| Nums of D | 減衰（速度に比例）の数（ばねと並列に接続し１〜１０の整数） |
| Mass (M0,M1,...M9) | 質量（kg） |
| Force | 質量に作用する周期的な外力を Fsin(2πft) とする、振幅 F(N) |
| Stiff.(K0,K1,...K9) | ばね定数（N/m） |
| Damp. (D0,D1,...D9) | 速度に比例する減衰（N*sec/m） |
| Connect to | ばねが接続される質量の番号（０ 〜 質量の数−１ の整数）
固定端の場合は −１とする |
| Display M# | 計算結果を表示する質量の番号（０ 〜 質量の数−１ の整数） |
for (i=0; i<nd; i++){
e1=(int)(stiff[i][1]);/* connected
end to mass #e1 */
e2=(int)(stiff[i][2]);/* connected
end to mass #e2 */
if (e1==-1)
c0[e2][e2]=c0[e2][e2]+damp[i][0];
else if (e2==-1)
c0[e1][e1]=c0[e1][e1]+damp[i][0];
else{
c0[e1][e2]=c0[e1][e2]-damp[i][0];
c0[e2][e1]=c0[e2][e1]-damp[i][0];
c0[e1][e1]=c0[e1][e1]+damp[i][0];
c0[e2][e2]=c0[e2][e2]+damp[i][0];
}
連立１次方程式を [ A + jB ]{ x }={ P }
ここで、
[ A ] ばねによる係数マトリックス
[ B ] 速度に比例する減衰による係数マトリックス
{ x } 質量の変位ベクトル
{ P } 質量に働く外力ベクトル
とすると変位ベクトルは、
{ x }=[ C + jD ]{ P }=[ E + jF ]
[ C + jD ]を求めるために、マトリックス[ A ]及び[ B ]を以下のように配置
し、逆マトリックスを計算
[A - B I 0 ]
[B A 0 I ]
ここで、
[ I ]は単位マトリックス
以下のプログラムでは、m1[][]に[ A ]、[ B ]及び[ I ]を配置し、
逆マトリックスをガウスの消去法により計算
for (j=0; j<nm; j++){
for (i=0; i<2*nm; i++){
i0[ i ][ j ]=0;
i0[ j ][ j ]=1;
}
}
for (i=0; i<nm; i++){
for (j=0; j<nm;
j++){
if (i==j) //
w1=2πf
m1[i][i]=k0[i][i]-w1*w1*mass[i];//
system matrix [ A ]
else{
m1[i][j]=k0[i][j];
m1[i][i]=k0[i][i]-w1*w1*mass[i];//
system matrix [ A ]
}
}
for (k=0; k<nm; k++){
m1[nm+i][k]=w1*c0[i][k];// system
matrix [ B ]
}
}
for (i=0; i<nm; i++){
for (j=0; j<nm; j++){
m1[i][nm+j]= -m1[nm+i][j];
m1[nm+i][nm+j]= m1[i][j];
}
}
for (i=0; i<2*nm; i++){
nj= 2*nm;
for (j=nj; j<4*nm; j++){
m1[i][j]=0;
}
m1[i][2*nm+i]=1;
}
for (k=0; k<2*nm; k++){
pivi=m1[k][k];
k1=k+1;
for (j=k1; j<4*nm; j++){
m1[k][j]=(m1[k][j])/pivi;
}
for (i=0; i<2*nm; i++){
if(i !=k){
w=m1[i][k];
for (j=k1; j<4*nm; j++){
m1[i][j]=m1[i][j]-w*m1[k][j];
}
}
}
}
// end of inversion
以下省略