/*! * ガウス・ザイデル法による連立方程式Ax = bの計算 * @param[in] a n×n正値対称行列 * @param[in] b n次元定数ベクトル * @param[inout] x 引数として初期値を入れておき、ルーチン内で解を代入 * @param[inout] max_iter 最大反復回数 * @retval true 解を求めるのに成功 * @retval false 解なし */ static bool GaussSeidel(vector< vector<double> > &A, vector<double> &b, vector<double> &x, const int n, int max_iter = 1000) { vector<double> x_old; x_old.resize(n); double denom, dev, max_dev; int i, j, iter = 0; // 行列AとベクトルBの正規化 for(i = 0; i < n; i++){ denom = A[i][i]; if(denom < RX_FEQ_EPS) return false; b[i] /= denom; for(j = 0; j < n; j++) A[i][j] /= denom; } // ガウス・ザイデル反復の実行 while(true){ for(i = 0; i < n; i++){ x_old[i] = x[i]; x[i] = 0; for(j = 0; j < n; j++){ if(j != i){ x[i] -= A[i][j]*x[j]; } } x[i] += b[i]; } // 収束判定 max_dev = fabs(x_old[0]-x[0])/x[0]; for(i = 1; i < n; i++){ dev = fabs(x_old[i]-x[i])/x[i]; max_dev = (dev > max_dev) ? dev : max_dev; } if(max_dev <= RX_FEQ_EPS){ return true; } else{ iter++; if(iter > max_iter){ return false; } } } } |