java ガウスの消去法

Question

Javaで、ガウスの消去法を行うプログラムを書いているのですが、答えが違ってしまいます、、、どこが間違っているのでしょうか？

augMatrix.matrixはArrayList<ArrayList<Double>>型で、拡大係数行列で、

2.0 4.0 -2.0 -4.0 -10.0
-6.0 -12.0 12.0 24.0 60.0
4.0 2.0 2.0 -4.0 8.0
2.0 -4.0 -2.0 4.0 6.0

の行列が入っています

ans.vectorはArrayList<Double>型で、解を格納するものです。

今のプログラムでは、解が

10.0 -5.0 -5.0 2.5

となってしまうのですが、正しい解は

2.0 -1.0 3.0 1.0

です。

わかりずらいプログラムだと思いますが、よろしくお願いします。

Vector solve() {
		
		Vector ans = new Vector();
		ArrayList<ArrayList<Double>> b = new ArrayList<ArrayList<Double>>();
		//行列の行数
		int N=augMatrix.matrix.size();
		int pivot;
        double p,q,m;
        
        //初期化
        for(int i=0;i<1;i++) {
        	ArrayList<Double> t = new ArrayList<Double>();
        	for(int j=0;j<N+1;j++) {
        		t.add(j,0.0);
        	}
        	b.add(t);
        }
        
        //初期化
        for(int i=0;i<N;i++) {
        	ans.vector.add(0.0);
        }

        for(int i=0;i<N;i++){
        	m=0;
            pivot=i;

            for(int l=i;l<N;l++) {
            	if(Math.abs(augMatrix.matrix.get(l).get(i))>m) { 
            		//i列の中で一番値が大きい行を選ぶ 
            		m=Math.abs(augMatrix.matrix.get(l).get(i));
                    pivot=l;
            	}
            }

            if(pivot!=i) {
            	
            	//pivotがiと違えば、行の入れ替え             
            	for(int j=0;j<N+1;j++) {
                       
            		 b.get(0).set(j,augMatrix.matrix.get(i).get(j));                        	
                     augMatrix.matrix.get(i).set(j,augMatrix.matrix.get(pivot).get(j));
                     augMatrix.matrix.get(pivot).set(j,b.get(0).get(j));
            	 }
             }
        }

        for(int k=0;k<N;k++) {
        	//対角要素を保存
        	p=augMatrix.matrix.get(k).get(k);
        	//対角要素は１になることがわかっているから
        	augMatrix.matrix.get(k).set(k,1.0);
     
        	for(int j=k+1;j<N+1;j++) {
        		augMatrix.matrix.get(k).set(j,(augMatrix.matrix.get(k).get(j))/p);
            }

            for(int i=k+1;i<N;i++) {		
            	q=augMatrix.matrix.get(i).get(k);
                for(int j=k+1;j<N+1;j++) {     	
                	augMatrix.matrix.get(i).set(j,-q*(augMatrix.matrix.get(k).get(j)));
                }       
                	//0となることがわかっているところ
                    augMatrix.matrix.get(i).set(k,0.0);
                }
        }

        //解の計算
        for(int i=N-1;i>=0;i--) {	
        	ans.vector.set(i,(augMatrix.matrix.get(i).get(N)));
            for(int j=N-1;j>i;j--) {
            	ans.vector.set(i,-1*(augMatrix.matrix.get(i).get(j))*ans.vector.get(j));
            }
        }

        //最終的な行列の表示
        result1(augMatrix.matrix);
        
        //解の出力
        result2(ans.vector);
     
		return ans;
	}

void result1(ArrayList<ArrayList<Double>> matrix) {
		for(int i=0;i<matrix.size();i++) {
			for(int j=0;j<matrix.get(0).size();j++) {
				System.out.print(matrix.get(i).get(j)+" ");
			}
		System.out.println("");
		}
}

void result2(ArrayList<Double> vector) {
		System.out.println("解は");
        for(int i=0;i<vector.size();i++) {
        	System.out.println(vector.get(i));
        }
}

Comments

[jimbe]

どのような値を入れてどのような答えになったのでしょうか。
teratail はデバッグを依頼する場所では無いと思います。
公式に公開されている各クラスやメソッドが仕様通りに動作していないのならともかく、ご自身で書いたプログラムのデバッグをまずご自身で行っていなければ、依頼と思われます。

プログラムの進行とともに変化する各変数の値の遷移をどこまで調査し、何がどう分からずにご質問に至ったのでしょうか。

「わかりずらいプログラムだと思」われるのでしたら、分かり易くなるよう修正されては如何でしょう。
インデントを整えるだけでもかなり見栄えは違うと思いますが。
(タブとスペースの混在でインデントを付けていると teratail では崩れやすいですので、どちらかに統一されると良いかと思います。)
また、ローカル変数をメソッドの最初で定義する書式は、ローカル変数のスコープを最大にし、(一定以上長いメソッドの場合は特に)変数値の追跡をむやみに難しくさせます。
augMatrix.matrix のような広域な変数の使用を見直すのと同様、処理単位毎にメソッドを分ける等して各変数のスコープをなるべく狭くすることも、分かり易いプログラムにするコツのように思います。

[hanao_]

すみませんでした、、、自分なりにコードを編集してみました。

変数のスコープについては今まで意識したことがなかったので、これから気を付けたいと思います。
勉強になります。

Answer 1

+       dump
        for(int k=0;k<N;k++) {
            //対角要素を保存
            p=augMatrix.matrix.get(k).get(k);
            //対角要素は１になることがわかっているから
            augMatrix.matrix.get(k).set(k,1.0);

            for(int j=k+1;j<N+1;j++) {
                augMatrix.matrix.get(k).set(j,(augMatrix.matrix.get(k).get(j))/p);
            }
+           dump

            for(int i=k+1;i<N;i++) {        
                q=augMatrix.matrix.get(i).get(k);
                for(int j=k+1;j<N+1;j++) {         
                    augMatrix.matrix.get(i).set(j,-q*(augMatrix.matrix.get(k).get(j)));
                }       
                    //0となることがわかっているところ
                    augMatrix.matrix.get(i).set(k,0.0);
                }
+               dump
        }

上記3か所にダンプ用コードを入れて実行し、配列の中身を出力した結果か以下になります。

---dump data start---
[-6.0, -12.0, 12.0, 24.0, 60.0]
[2.0, 4.0, -2.0, -4.0, -10.0]
[4.0, 2.0, 2.0, -4.0, 8.0]
[2.0, -4.0, -2.0, 4.0, 6.0]
---dump data end---
---dump data start---
[1.0, 2.0, -2.0, -4.0, -10.0]
[2.0, 4.0, -2.0, -4.0, -10.0]
[4.0, 2.0, 2.0, -4.0, 8.0]
[2.0, -4.0, -2.0, 4.0, 6.0]
---dump data end---
---dump data start---
[1.0, 2.0, -2.0, -4.0, -10.0]
[0.0, -4.0, 4.0, 8.0, 20.0]
[0.0, -8.0, 8.0, 16.0, 40.0]
[0.0, -4.0, 4.0, 8.0, 20.0]
---dump data end---
---dump data start---
[1.0, 2.0, -2.0, -4.0, -10.0]
[0.0, 1.0, -1.0, -2.0, -5.0]
[0.0, -8.0, 8.0, 16.0, 40.0]
[0.0, -4.0, 4.0, 8.0, 20.0]
---dump data end---
---dump data start---
[1.0, 2.0, -2.0, -4.0, -10.0]
[0.0, 1.0, -1.0, -2.0, -5.0]
[0.0, 0.0, -8.0, -16.0, -40.0]
[0.0, 0.0, -4.0, -8.0, -20.0]
---dump data end---
---dump data start---
[1.0, 2.0, -2.0, -4.0, -10.0]
[0.0, 1.0, -1.0, -2.0, -5.0]
[0.0, 0.0, 1.0, 2.0, 5.0]
[0.0, 0.0, -4.0, -8.0, -20.0]
---dump data end---
---dump data start---
[1.0, 2.0, -2.0, -4.0, -10.0]
[0.0, 1.0, -1.0, -2.0, -5.0]
[0.0, 0.0, 1.0, 2.0, 5.0]
[0.0, 0.0, 0.0, 8.0, 20.0]
---dump data end---
---dump data start---
[1.0, 2.0, -2.0, -4.0, -10.0]
[0.0, 1.0, -1.0, -2.0, -5.0]
[0.0, 0.0, 1.0, 2.0, 5.0]
[0.0, 0.0, 0.0, 1.0, 2.5]
---dump data end---
---dump data start---
[1.0, 2.0, -2.0, -4.0, -10.0]
[0.0, 1.0, -1.0, -2.0, -5.0]
[0.0, 0.0, 1.0, 2.0, 5.0]
[0.0, 0.0, 0.0, 1.0, 2.5]
---dump data end---
---dump data start---
[1.0, 2.0, -2.0, -4.0, -10.0]
[0.0, 1.0, -1.0, -2.0, -5.0]
[0.0, 0.0, 1.0, 2.0, 5.0]
[0.0, 0.0, 0.0, 1.0, 2.5]
---dump data end---

正しく消去法ができていると思いますか？

Answer 2

argMatrix.matrix がコード中にたくさんあって見にくいいので、
solve() の先頭で、変数 m にしています。

import java.util.ArrayList;
import java.util.Arrays;

class Matrix {
	ArrayList<ArrayList<Double>> matrix;

	Matrix(Double[][] d) {
		matrix = new ArrayList<>();
		for (var e : d)
			matrix.add(new ArrayList<Double>(Arrays.asList(e)));
	}
}
	
class Gauss {
	Double[][] d = {
		{  2.0,   4.0, -2.0, -4.0, -10.0 },
		{ -6.0, -12.0, 12.0, 24.0,  60.0 },
		{  4.0,   2.0,  2.0, -4.0,   8.0 },
		{  2.0,  -4.0, -2.0,  4.0,   6.0 },
	};
	Matrix augMatrix = new Matrix(d);
	int N = augMatrix.matrix.size();

	ArrayList<Double> solve() {
		ArrayList<ArrayList<Double>> m = augMatrix.matrix;
		for (int k = 0; k < N; k++) {
			double t = Math.abs(m.get(k).get(k));
			int p = k;
			for (int i = k+1; i < N; i++)
				if (m.get(i).get(k) > t) {
					t = Math.abs(m.get(i).get(k));
					p = i;
				}
			if (p != k)
				for (int j = 0; j <= N; j++) {
					t = m.get(k).get(j);
					m.get(k).set(j, m.get(p).get(j));
					m.get(p).set(j, t);
				}
			for (int i = k+1; i < N; i++) {        
				t = m.get(i).get(k) / m.get(k).get(k);
				for (int j = k; j <= N; j++)        
					m.get(i).set(j, m.get(i).get(j) - t * m.get(k).get(j));      
			}
		}

		for (int k = N; --k > 0; )    
			for (int i = k; --i >= 0; ) {
				double t = m.get(i).get(k) / m.get(k).get(k);
				for (int j = 0; j <= N; j++)        
					m.get(i).set(j, m.get(i).get(j) - t * m.get(k).get(j));
			}

		for (int k = 0; k < N; k++) {
			m.get(k).set(N, m.get(k).get(N) / m.get(k).get(k));
			m.get(k).set(k, 1.0);
		}

		result1(m);

		ArrayList<Double> vector = new ArrayList<>();
		for (int i = 0; i < N; i++)
			vector.add(m.get(i).get(N));
		return vector;
	}

	void result1(ArrayList<ArrayList<Double>> matrix) {
		for (int i = 0; i < matrix.size(); i++) {
			for (int j = 0; j < matrix.get(0).size(); j++)
				System.out.printf("%8.3f", matrix.get(i).get(j));
			System.out.println();
		}
	}
}

class Main {
	public static void main(String[] args) {
		ArrayList<Double> ans = new Gauss().solve();
		System.out.println();
		for (var e : ans) System.out.printf("%8.3f", e);
		System.out.println();
	}
}

実行結果

   1.000   0.000   0.000   0.000   2.000
   0.000   1.000   0.000   0.000  -1.000
   0.000   0.000   1.000   0.000   3.000
   0.000   0.000   0.000   1.000   1.000

   2.000  -1.000   3.000   1.000

追記
上記のコードに間違いがありました。

if (m.get(i).get(k) > t) { を
if (Math.abs(m.get(i).get(k)) > t) { に修正してください。

今回のデータでは、たまたま結果が正しくなっていました。

追記2

            if (p != k)
                for (int j = 0; j <= N; j++) {
                    t = m.get(k).get(j);
                    m.get(k).set(j, m.get(p).get(j));
                    m.get(p).set(j, t);
                }

これは、if (p != k) java.util.Collections.swap(m, k, p); と書き換え可能です。

Answer 3

二重コレクション(List<List<Double>>)は扱い難いので、クラスとしてまとめる形で配列に戻し(?)てみました。
division 等でコメントにしてある println や printMatrix を有効にすれば動作が確認し易いかと思いますので、ご自身のコードでも同じように表示してみると違う個所が分かるかと思います。

package teratail_java.q371994;

import java.io.PrintStream;

public class Q371994 {

  public static void main(String[] args) {
/*
    GaussianEliminationMatrix matrix = new GaussianEliminationMatrix(
        new double[] {
            2, 4,-2, 8, //2[x1]+4[x2]-2[x3]=8
            1, 2, 1, 6, // [x1]+2[x2]+ [x3]=6
            1, 3, 2, 9, // [x1]+3[x2]+2[x3]=9
        });
*/
    GaussianEliminationMatrix matrix = new GaussianEliminationMatrix(
        new double[] {
             2.0,  4.0,-2.0,-4.0, -10.0,
            -6.0,-12.0,12.0,24.0,  60.0,
             4.0,  2.0, 2.0,-4.0,   8.0,
             2.0, -4.0,-2.0, 4.0,   6.0,
        });
    matrix.printMatrix("start ---");

    matrix.forwardElimination();
    matrix.printMatrix("forward ---");

    double[] x = matrix.backwardSubstitution();

    System.out.println("解は");
    for(int i=0; i<x.length; i++) {
      System.out.println(x[i]);
    }

    //Vector<Double> ans = new Vector<Double>(x.length);
    //for(Double b : x) ans.add(b);
    //return ans;
  }
}

final class GaussianEliminationMatrix {
  private int size;
  private double[][] a_b;

  GaussianEliminationMatrix(double[] matrix) {
    size = (int)(Math.sqrt(matrix.length));
    if(matrix.length != size*size+size) throw new IllegalArgumentException();

    int w = size + 1;
    a_b = new double[size][w];
    for(int i=0; i<matrix.length; i++) {
      a_b[i/w][i%w] = matrix[i];
    }
  }

  double getA(int i, int j) {
    if(j >= size) throw new ArrayIndexOutOfBoundsException("j");
    return a_b[i][j];
  }
  double getB(int i) { return a_b[i][size]; }

  //前進消去
  void forwardElimination() {
    for(int k=0; k<size-1; k++) {
      int p = pivoting(k);
      swap(k, p);
      division(k);
      for(int i=k+1; i<size; i++) {
        subtraction(k, i);
      }
    }
  }
  //k行k列から下の中で一番絶対値が大きい行
  private int pivoting(int k) {
    int pivot = 0;
    double m = 0;
    for(int l=k; l<size; l++) {
      if(m < Math.abs(getA(l,k))) {
        m = Math.abs(getA(l,k));
        pivot = l;
      }
    }
    return pivot;
  }
  //行の入れ替え
  private void swap(int k, int p) {
    double t;
    for(int i=0; i<size+1; i++) {
      t = a_b[k][i];
      a_b[k][i] = a_b[p][i];
      a_b[p][i] = t;
    }
    //System.out.println("swap: k="+k+", p="+p);
    //printMatrix("swap");
  }
  //行の徐算(k行k列を1にする)
  private void division(int k) {
    double q = getA(k,k);
    //System.out.println("division: k="+k+", q="+q);
    for(int l=k; l<size+1; l++) {
      a_b[k][l] /= q;
    }
    //printMatrix("division");
  }
  //行の減算(i行k列を0にする)
  private void subtraction(int k, int i) {
    double p = getA(i,k) / getA(k,k);
    //System.out.println("subtraction: k="+k+", i="+i+", p="+p);
    for(int l=k; l<size+1; l++) {
      a_b[i][l] -= p * a_b[k][l];
    }
    //printMatrix("subtraction");
  }

  //後退代入
  double[] backwardSubstitution() {
    double[] x = new double[size];

    x[size-1] = getB(size-1) / getA(size-1,size-1);
    for(int i=size-2; i>=0; i--) {
      double s=0;
      for(int j=i+1; j<size; j++) {
        s += getA(i,j) * x[j];
      }
      x[i] = (getB(i) - s) / getA(i,i);
    }

    return x;
  }

  //表示
  void printMatrix(String header) {
    System.out.println(header);
    printMatrix(System.out);
  }
  void printMatrix(PrintStream out) {
    for(int i=0; i<size; i++) {
      for(int j=0; j<size+1; j++) {
        System.out.print(a_b[i][j] + " ");
      }
      System.out.println();
    }
  }
}

実行結果

start ---
2.0 4.0 -2.0 -4.0 -10.0 
-6.0 -12.0 12.0 24.0 60.0 
4.0 2.0 2.0 -4.0 8.0 
2.0 -4.0 -2.0 4.0 6.0 
forward ---
1.0 2.0 -2.0 -4.0 -10.0 
0.0 1.0 -0.25 -1.5 -3.25 
0.0 0.0 1.0 0.35294117647058826 3.3529411764705883 
0.0 0.0 0.0 3.2941176470588234 3.2941176470588234 
解は
2.0
-1.0
3.0
1.0