#code(C){{
/*!
 * ガウス・ザイデル法による連立方程式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;
			}
		}
	}
}
}}

トップ   編集 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS