質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.12%

C言語 Linux

解決済

回答 3

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 2,496

score 26

以前にも同じ質問をさせていただいたのですが、解決できなかったので再送します。

Gnome端末を使っています。

x1+2x2+3x3=6
3x1+10x2-4x3=-29
-2x1-4x2+x3=9
という連立方程式をGauss-seidel法で求める問題です。

下記のようにプログラミングすると、

ILL CONDITION  i=0
Iteration    x1    x2    x3
1  1.0000  1.1034  1.7126
2  -0.2241  0.7406  1.2793
3  0.1135  0.8353  1.3964
4  0.0233  0.8098  1.3651
5  0.0475  0.8166  1.3735
6  0.0410  0.8148  1.3713
7  0.0428  0.8153  1.3719
8  0.0423  0.8152  1.3717
9  0.0424  0.8152  1.3717
10  0.0424  0.8152  1.3717
11  0.0424  0.8152  1.3717

となり、x1=1、x=ー2、x3=3が得られません。

また、

3x1+x2+x3=10
x1+5x2+2x3=21
x1+2x2+5x3=30

を同じプログラミングで解こうとすると、

Iteration    x1    x2    x3
1  1.0000  0.9524  0.9032
2  0.8144  0.8752  0.9145
3  0.8210  0.8738  0.9144
4  0.8212  0.8738  0.9144
5  0.8212  0.8738  0.9144

のようになり、x1=1、x=2、x3=5が得られません。

“♯include<stdio.h>
♯include<math.h>
♯define eps 1.e-5
int main(void)
{

double a[3][3],b[3],x[3];
double s,w,anorm,xnorm;
int i,j,k,n;
scanf("%d",&n);
for(i=0;i<n;i++){
for(j=0;j<n;j++){
scanf("%lf",&a[i][j]);
}
scanf("%lf",&b[i]);}

for(i=0;i<n;i++){
s=0.0;
for(j=0;j<n;j++){
if(i==j) continue;
s+=fabs(a[i][j]);
}

if(fabs(a[i][i])<=s){
printf("ILL CONDITION  i=%d\n",i);
break;
}
}

for(i=0;i<n;i++){
x[i]=0.0;
}

printf("Iteration    x1    x2    x3\n");
for(k=1;;k++){
anorm=0.0;
xnorm=0.0;
for(i=0;i<n;i++){
w=b[i];
for(j=0;j<n;j++){
if(j==i) continue;
w -=a[i][j]*x[j];
}

w /=b[i];
anorm+=fabs(x[i]-w);
xnorm+=fabs(w);
x[i]=w;
}

printf("%6d",k);
for(i=0;i<n;i++){
printf("  %.4f",x[i]);
}
printf("\n");
if(anorm/xnorm<eps)
break;
}
return 0;
}

御指摘いただきたいです。

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • KSwordOfHaste

    2017/01/17 16:13

    コードを載せるときは```で囲み字下げしましょう。

    キャンセル

回答 3

checkベストアンサー

+1

なんとなく質問者さんが理解されていないように感じるため、修正したコードの全文を掲載します。

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define eps 1.e-5
int main(void)
{
  double a[3][3],b[3],x[3];
  double s,w,anorm,xnorm;
  int i,j,k,n;
  scanf("%d", &n);
  for(i=0; i<n; i++){
    for(j=0; j<n; j++){
      scanf("%lf", &a[i][j]);
    }
    scanf("%lf", &b[i]);
  }

  for(i=0;i<n;i++){
    s=0.0;
    for(j=0;j<n;j++){
      if(i==j) continue;
      s+=fabs(a[i][j]);
    }

    if(fabs(a[i][i])<=s){
      printf("ILL CONDITION  i=%d\n",i);
      exit(1);
    }
  }

  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      printf("%lfx%d ", a[i][j], j+1);
    }
    printf("= %lf\n", b[i]);
  }

  for(i=0; i < n; i++){
    x[i] = 0.0;
  }

  printf("Iteration    x1    x2    x3\n");
  for(k=1; ; k++){
    anorm = 0.0;
    xnorm = 0.0;
    for(i = 0; i < n; i++){
      w=b[i];
      for(j=0; j<n; j++){
        if(j==i) continue;
        w -= a[i][j]*x[j];
      }

//    w /= b[i]; // ←ここが間違っている
      w /= a[i][i];
      anorm += fabs(x[i]-w);
      xnorm += fabs(w);
      x[i] = w;
    }

    printf("%6d", k);
    for(i = 0; i < n; i++){
      printf("  %.4f", x[i]);
    }
    printf("\n");
    if(anorm/xnorm < eps) break;
  }
  return 0;
}

ガウス=ザイデル法の更新式をよく見なおしてください。

なお、ひとつ目の例は、
ILL CONDITION  i=0 
と出力されているように、ガウス=ザイデル法の収束条件(対角優位)を満たしていないため収束しません。

以上のように修正すれば、2つ目の例は収束します。

% ./a.out    
3
3 1 1 10
1 5 2 21
1 2 5 30
Iteration    x1    x2    x3
     1  3.3333  3.5333  3.9200
     2  0.8489  2.4622  4.8453
     3  0.8975  2.0824  4.9876
     4  0.9767  2.0096  5.0008
     5  0.9965  2.0004  5.0005
     6  0.9997  1.9998  5.0001
     7  1.0000  1.9999  5.0000
     8  1.0000  2.0000  5.0000

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2017/01/17 15:51

    お返事ありがとうございます。
    二つ目の例の場合でも、収束はしていますが、得たい解答(x1=1、x=2、x3=5)が得られていないです。どのようにすれば、正答が得られるか教えていただけると嬉しいです。

    キャンセル

  • 2017/01/17 15:53

    修正後
    Iteration x1 x2 x3
    ....
    8 1.0000 2.0000 5.0000
    となっていますよ。あっていると思うのですが。

    キャンセル

  • 2017/01/17 16:20


    申し訳ございません。
    あっていました。

    おこがましいですが、再び質問させていただきます。
    一つ目の例で収束しないということは、ガウス=ザイデル法では一つ目の例は解けないということなのでしょうか。

    キャンセル

  • 2017/01/17 16:24

    そのままでは解けません。
    これを解くために、
    Ax=b ではなく、A^T Ax =A^T b を解くテクニックがあるようです。

    キャンセル

0

インデントを整えたり、コメントを入れて、プログラムをわかりやすくしましょう。
バグが発見しやすくなります。

2番目の連立方程式は、解x1=1、x=2、x3=5が得られました。
1番目の連立方程式は、解がGauss-Seidel 法では発散するようです。

#include<stdio.h> 
#include<math.h> 


/* 計算誤差の許容値 */
#define eps (1.e-5) 


int main(void) 
{

    double a[3][3], b[3], x[3]; 
    double s, anorm, xnorm, sum; 
    int i, j, k;
    int n; /* 連立方程式の大きさ */
    double new_x; /* 新しい近似解 */


    /* データ入力 */

    scanf("%d", &n);

    for (i = 0; i < n; i++) { 
        for (j = 0; j < n; j++) { 
            scanf("%lf", &a[i][j]); 
        } 
        scanf("%lf", &b[i]);
    }


    for (i = 0; i < n; i++) { 
        s = 0.0;

        for (j = 0; j < n; j++) { 
            if (i == j) 
                continue; 

            s += fabs(a[i][j]); 
        }

        if (fabs(a[i][i]) <= s) { 
            printf("ILL CONDITION  i=%d\n", i);
            /* この条件のとき、近似解は発散? */
            break; 
        } 
    }

    /* 近似解の初期値代入 */
    for (i = 0; i < n; i++) { 
        x[i] = 0.0; 
    }


    printf("Iteration    x1    x2    x3\n"); 

    for (k = 1; ; k++) { 
        anorm = 0.0; 
        xnorm = 0.0; 

        for (i = 0; i < n; i++) { 
            sum = 0.0;

            for (j = 0; j < n; j++) { 
                if (j == i) 
                    continue;

                sum += a[i][j] * x[j]; 
            }

            /* 反復計算後の近似解 */
            new_x = 1.0 / a[i][i] * (b[i] - sum);
            /* 近似解の変化量を加算 */
            anorm += fabs(x[i] - new_x);
            /* 近似解の総和計算 */
            xnorm += fabs(new_x);
            /* 新しい近似解を代入 */
            x[i] = new_x; 
        }

        printf("%6d", k);

        for (i = 0; i < n; i++) { 
            printf("  %.4f", x[i]); 
        } 

        printf("\n");

        /* 計算終了条件 */
        if (anorm / xnorm < eps) 
            break; 
    } 
    return 0; 
}

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

0

ちょこっとwikiでGauss-seidel法をみただけですがGauss-seidel法のアルゴリズムとコードを見比べておかしいと思える点をコメントします。

  1. 前回の質問で正しいコードを間違ったものに直してます。
    w /= a[i][i]が正しいですが、これをw /= b[i]に修正してしまってます。コメントした方は「Gauss-seidel法については知らないがコードをみた感じ間違いではないか?」とおっしゃっているだけですので、その指摘が正しいかどうかは質問者さんご自身で判断すべきだったと思います。

  2. 収束条件
    Gauss-seidelで解が収束する充分条件は(必要条件かはわかりません)「係数行列の各行で非対角要素の絶対値の和が対角要素の絶対値よりも小さい場合」とwikiにあります。プログラムコードでもそれをチェックしてはいますが単にメッセージを表示するだけで実行は継続していますね?最初の方程式は1行目がこの条件を満たしていませんので"IL CONDITION"とメッセージが出ます。実行を継続してみるとおっしゃるとおり収束しないように見えます。エラーメッセージを出すだけでプログラムを中断しない論理となっているのはどういう意図でしょう?「充分条件は満たしていなくてもとにかくやってみる」ということならわかりますが・・・当然ながら本当に収束できないケースでは無限ループに陥ってしまいます。収束条件を満たさないケースで処理をあきらめたいなら「収束条件を満たさないことが判明した時点で」returnなりexitで実行を中断すべきでしょう。
    なお、収束条件を満たすように式を変形するような配慮をプログラム上でしてやれば解が求まるとは思いますが、そのような論理を入れると最早Gauss-seidel法ではなくGaussの消去法になってしまうのかも知れません。

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.12%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る