行列
いつもながらありがたい青木先生の電卓シリーズ(行列式、逆行列)。ネットにつながればいつでも計算できる。
こちらはeclipseを用いてJavaサンプルソースを丁寧な注釈とともに公開しているページ。そちらの逆行列用クラスを流用したソースが以下
import java.text.*;public class GyakuGyouretsu {
//こちらのページの公開ソースの行列表現部分のみを演算対象部分のみの表現に修正したものです。
//http://www.geocities.jp/java_sample_program/index.htm
//pivotは、消去演算を行う前に、対象となる行を基準とし、それ以降の
//行の中から枢軸要素の絶対値が最大となる行を見つけ出し、対象の行と
//その行とを入れ替えることを行う関数である。
void pivot(int k,double a[][],int N){
double max,copy;
//ipは絶対値最大となるk列の要素の存在する行を示す変数で、
//とりあえずk行とする
int ip=k;
//k列の要素のうち絶対値最大のものを示す変数maxの値をとりあえず
//max=|a[k][k]|とする
max=Math.abs(a[k][k]);
//k+1行以降、最後の行まで、|a[i][k]|の最大値とそれが存在する行を
//調べる
for(int i=k+1; i<N; i++){
if(max<Math.abs(a[i][k])){
ip=i;
max=Math.abs(a[i][k]);
}
}
if(ip!=k){
for(int j=0; j<2*N; j++){
//入れ替え作業
copy =a[ip][j];
a[ip][j]=a[k][j];
a[k][j] =copy;
}
}
}//ガウス・ジョルダン法により、消去演算を行う
void sweep(int k,double a[][],int N){
double piv,mmm;
//枢軸要素をpivとおく
piv=a[k][k];
//k行の要素をすべてpivで割る a[k][k]=1となる
for(int j=0; j<2*N; j++)
a[k][j]=a[k][j]/piv;
//
for(int i=0; i<N; i++){
mmm=a[i][k];
//a[k][k]=1で、それ以外のk列要素は0となる
//k行以外
if(i!=k){
//i行において、k列から2N-1列まで行う
for(int j=k; j<2*N; j++)
a[i][j]=a[i][j]-mmm*a[k][j];
}
}
}
//逆行列を求める演算
void gyaku(double a[][],int N){
for(int k=0; k<N; k++){
pivot(k,a,N);
sweep(k,a,N);
}
}
public static void main(String[] args) {
//DecimalFormatクラスにより小数点以下5桁まで表示
DecimalFormat ft=new DecimalFormat("##0.00000");
double b[][] ={{1,-1,-1,1,-1,1,1,-1},
{1,1,-1,-1,-1,-1,1,1},
{1,-1,1,-1,-1,1,-1,1},
{1,-1,-1,1,1,-1,-1,1},
{1,1,1,1,-1,-1,-1,-1},
{1,1,-1,-1,1,1,-1,-1},
{1,-1,1,-1,1,-1,1,-1},
{1,1,1,1,1,1,1,1}};
double[][] invb = inv(b);
String out ="元行列\n";
for(int i=0;i<b.length;i++){
for(int j=0;j<b[i].length;j++){
out += b[i][j] + "\t";
}
out += "\n";
}
System.out.println(out);
out="";
for(int i=0;i<invb.length;i++){
for(int j=0;j<invb[i].length;j++){
out += invb[i][j] + "\t";
}
out += "\n";
}
System.out.println(out);
}
public static double[][] matPlusEmat(double[][] m){
double[][] ret = new double[m.length][m[0].length*2];
for(int i=0;i<ret.length;i++){
for(int j=0;j<m[i].length;j++){
ret[i][j]=m[i][j];
if(i==j){
ret[i][m[i].length+j]=1;
}else{
ret[i][m[i].length+j]=0;
}
}
}
return ret;
}
public static double[][] inv(double[][] b){
double[][] ret = new double[b.length][b[0].length];
double[][] a = matPlusEmat(b);
int N=a.length;
//Gyakugyouretuクラスのインスタンスの生成
GyakuGyouretsu g=new GyakuGyouretsu();
//gyakuの呼び出し
g.gyaku(a,N);
for(int i=0; i<N; i++){
for (int j=N; j<2*N; j++){
ret[i][j-N]=a[i][j];
}
}
return ret;
}
}
