C言語 Linux

Question

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

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＝１、x=ー２、x3=３が得られません。

また、

３x1+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＝１、x=２、x3=５が得られません。

“♯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;
}

”

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

Comments

[KSwordOfHaste]

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

Answer 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
と出力されているように、ガウス＝ザイデル法の収束条件（対角優位）を満たしていないため収束しません。

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

% ./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

Comments

[uv-]

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

[MasashiKimura]

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

[uv-]

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

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

[MasashiKimura]

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

Answer 2

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

1. 前回の質問で正しいコードを間違ったものに直してます。

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

1. 収束条件

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

Answer 3

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

2番目の連立方程式は、解x1＝１、x=２、x3=５が得られました。
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; 
}