本業だけじゃ生活できないので、ときどき他学科のプログラミング演習の課題を受託開発をすることで、見返りに唐揚げ定食とかうまかっちゃんとかよくおごってもらっている。
なぜかこの時期に「行列式の値」と「逆行列」を作成するプログラムを毎年書いている気がして、その都度作るのが非常に面倒臭いので、今後も困らないようにアップしておく。
行列式の値の求め方
・2 次の場合
たすきがけで簡単に解ける。
double a[2][2]={{1,2},{4,-1}};
double det=0.0;
det=a[0][0]*a[1][1];
det-=a[0][1]*a[1][0];
printf("%f\n",det); // -> -9.000000
・3 次の場合
サラスの方法で簡単に解ける。
double a[3][3]={{1,8,9},{-3,2,1},{4,1,5}};
double det=0.0;
det=a[0][0]*a[1][1]*a[2][2];
det+=a[1][0]*a[2][1]*a[0][2];
det+=a[2][0]*a[0][1]*a[1][2];
det-=a[2][0]*a[1][1]*a[0][2];
det-=a[1][0]*a[0][1]*a[2][2];
det-=a[0][0]*a[2][1]*a[1][2];
printf("%f\n",det); // -> 62.000000
・4 次以上の場合
通常の計算では 3×3 以下の余因子行列を作成して、符号を掛けた上で足し合わせるやり方を使うと思う。
これをプログラムで実装するには、再帰または多段ループが必要になり計算量も多い。
特に 5 次以降になると面倒。
でも、上三角行列を作って、対角部分 (左上から右下) を掛け合わせれば、簡単に行列式の値を求めることができる。
このやり方だと非常に簡単。
double a[4][4]={{2,-2,4,2},{2,-1,6,3},{3,-2,12,12},{-1,3,-4,4}};
double det=1.0,buf;
int n=4; //配列の次数
int i,j,k;
//三角行列を作成
for(i=0;i<n;i++){
for(j=0;j<n;j++){
if(i<j){
buf=a[j][i]/a[i][i];
for(k=0;k<n;k++){
a[j][k]-=a[i][k]*buf;
}
}
}
}
//対角部分の積
for(i=0;i<n;i++){
det*=a[i][i];
}
printf("%f\n",det); // -> 120.000000
逆行列の求め方
何といってもやはり掃き出し法が簡単。
2 次から n 次まで対応。
double a[4][4]={{1,2,0,-1},{-1,1,2,0},{2,0,1,1},{1,-2,-1,1}}; //入力用の配列
double inv_a[4][4]; //ここに逆行列が入る
double buf; //一時的なデータを蓄える
int i,j,k; //カウンタ
int n=4; //配列の次数
//単位行列を作る
for(i=0;i<n;i++){
for(j=0;j<n;j++){
inv_a[i][j]=(i==j)?1.0:0.0;
}
}
//掃き出し法
for(i=0;i<n;i++){
buf=1/a[i][i];
for(j=0;j<n;j++){
a[i][j]*=buf;
inv_a[i][j]*=buf;
}
for(j=0;j<n;j++){
if(i!=j){
buf=a[j][i];
for(k=0;k<n;k++){
a[j][k]-=a[i][k]*buf;
inv_a[j][k]-=inv_a[i][k]*buf;
}
}
}
}
//逆行列を出力
for(i=0;i<n;i++){
for(j=0;j<n;j++){
printf(" %f",inv_a[i][j]);
}
printf("\n");
}
/* 出力
2.000000 2.000000 -1.000000 3.000000
-4.000000 -5.000000 3.000000 -7.000000
3.000000 4.000000 -2.000000 5.000000
-7.000000 -8.000000 5.000000 -11.000000
*/
プログラムは Microsoft Visual Studio 2008 および Fedora 9 の gcc で実行確認済み。
手元に Linux のマシンがなかったので、稼働中の Web サーバを使って動作確認した。
このブログのフィードを取得
