1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
/*!
 * ガウス・ザイデル法による連立方程式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
Last-modified: 2011-10-27 (木) 15:09:18 (3859d)