ガウス消去法

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら　　※コードを貼り付ける場合は [code] と [/code] で囲って下さい。詳しくはこちら
返信
ami

ガウス消去法

#1

投稿記事 by ami » 5年前

大学の課題でガウス消去法のプログラムと、ピボット操作を行うガウス消去法による逆行列のプログラムを作る問題があります。

講義内に配布されたスライドの疑似言語を使って作ってみたのですが、上手く答えを出すことができません。(EasyIDECを使用して実行しました)
また、解を求める過程というのも出さなければならないのですが、どうすれば良いのでしょうか？
自分ではどこがまちがっているのかわからないため、質問させていただきました。
よろしくお願いします！

ガウス消去法：

コード: 

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

#define N	3
#define EPS .001

int main( int argc, char ** argv ){
	double a[N][N+1] = {
		{ 2, -4, 6, 24 },
		{ -1, 3, -5, -20 },
		{ 1, 2, 3, 4 }
	};
	double pivot, ratio;
	int i, j, k, l;
	
	//前進
	for( i=0; i<N; i++ ){
		pivot = a[i][i];
		if(fabs( pivot ) < EPS ){
			printf( "ピボットが許容範囲誤差以下\n" );
			return 1;
		}
		
		for( k=i+1; k<N+1; k++ ){
			ratio = a[k][i] / pivot;
			for(j=i; j<N+1; j++){
			a[k][j] = a[k][j] - a[i][j] * ratio;
			}
		}
	}
	
	//後退
		for( i=N-1; i>=0; i-- ){
			for( j=N-1; j>=i+1; j-- ){
				a[i][N] = a[i][N] - a[j][N] * a[i][j];
				a[i][j] = 0.0;
			}
			a[i][N] = a[i][N] / a[i][i];
			a[i][i] = 1.0;
		}

			
		
			printf( "X = %6.2lf\n", a[0][N] );
			printf( "Y = %6.2lf\n", a[1][N] );
			printf( "Z = %6.2lf\n", a[2][N] );

	return a[N][N+1];
	}
ピボット操作を行うガウス消去法(逆行列)

コード: 

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

#define N	3
#define EPS .00001

void print_hline(void);
void print_array(double b[N][N+1]);

static void irekae(double b[N][N+1],int i)
{
	int m,j,k;
	double work;
	
	m=i;
	for(k=i+1;k<N;k++)
		if(fabs (b[m][i])<fabs(b[k][i])) m=k;

	for(j=0; j<N+1;j++){
		work= b[m][j];b[m][j]=b[i][j];b[i][j]=work;
	}
}

int main( int argc, char ** argv ){
	double a[N][N] = {
		{ 0.0, 2.0, 1.0 },
		{ 4.0, 4.0, 1.0 },
		{ -2.0, 1.0, -2.0 }
	};
	double b[N][2*N];
	double pivot, ratio;
	int x, y, i, j, k;
		for( y=0; y<N; y++ ){
		for(x=0; x<N; x++){
			b[ y ][ x     ] = a[ y ][ x ];
			b[ y ][ x + N ] = 0.0;
		}
		b[ y ][ y + N ] = 1.0;
	}
	
	//前進
	for( i=0; i<N; i++ ){
		pivot = b[i][i];
		if(fabs( pivot ) < EPS ){
			printf( "ピボットが許容範囲誤差以下\n" );
			return 1;
		}
		
		for( k=i+1; k<N+1; k++ ){
			ratio = b[k][i] / pivot;
			for(j=i; j<N*2; j++){
			b[k][j] = b[k][j] - b[i][j] * ratio;
			}
		}
	}
	
	//後退
		for( i=N-1; i>=0; i-- ){
			for( j=N-1; j>=i+1; j-- ){
				for(k=N; k<N*2; k++){
				b[i][k] = b[i][k] - b[j][k] * b[i][j];
				}
			b[i][j] = 0.0;
		}
			for(k=N; k<N*2; k++){
			b[i][k] = b[i][k] / b[i][i];
		}
		b[i][i] = 1.0;
		}

			
		
		for( y=0; y<N; y++ ){
			for( x=N; x<2*N; x++ )
			printf( "%7.3lf", b[ y ][ x ] );
			printf( "\n" );
		}

	return b[N][2*N];
	}
ページトップ
かずま

Re: ガウス消去法

#2

投稿記事 by かずま » 5年前

ami さんが書きました: ↑
5年前
 自分ではどこがまちがっているのかわからないため、質問させていただきました。
間違っているところは、前進の

コード: 

		for( k=i+1; k<N+1; k++ ){
これだと、k が N になってもループを抜けなくて、
a[N][ i] という配列の範囲外の要素を参照してしまいます。
k<N+1 ではなく、k < N にしましょう。
逆行列の方も同じ。

最後の return でも a[N][N+1] という範囲外を参照しています。
また double の値を返そうとするのは無意味です。
main の成功を意味する return 0; にしましょう。
逆行列の方も同じ。

逆行列の方は、pivot が 0 になっているので、

コード: 

		if(fabs( pivot ) < EPS ){
			irekae(b, i);
			pivot = b[i][i];
		}
のように入れ替えを実行しましょう。
b は main で、double b[N][2*N]; と宣言されているので、
irekae の引数でもそのように宣言し、for文で、
j<N+1 ではなく、j < 2*N にしましょう。
ami さんが書きました: ↑
5年前
 また、解を求める過程というのも出さなければならないのですが、どうすれば良いのでしょうか？
解を求める過程の表示は、おそらく print_array を使うのでしょうが、
これは逆行列の方だけでいいんでしょうか？
ページトップ
ami

Re: ガウス消去法

#3

投稿記事 by ami » 5年前

ありがとうございます。早速直してみようと思います。
解を求める過程ですが、両方とも求めたいです。
ページトップ
ami

Re: ガウス消去法

#4

投稿記事 by ami » 5年前

すいません、逆行列のほうなのですが、この解だとどうしても「ピボットが許容範囲誤差以下」になるのですが、(値を変えると答えは出てくる)これであっているのでしょうか？
ページトップ
かずま

Re: ガウス消去法

#5

投稿記事 by かずま » 5年前

ami さんが書きました: ↑
5年前
 すいません、逆行列のほうなのですが、この解だとどうしても「ピボットが許容範囲誤差以下」になるのですが、(値を変えると答えは出てくる)これであっているのでしょうか？
ピボットは、0.0, 4.0, -2.0 で、0.0 は許容範囲以下になっているので、
「入れ替えを実行しましょう」と言ったはずですが、
「早速直してみたプログラム」はどんなものでしょうか？

また、「値を変えると答えは出てくる」とは
どの値をどのように変えたのですか？
質問には、具体的な情報を付けてください。

返信が遅くなりましたので、解を求める過程を追加したプログラムを示します。

コード: 

#include <stdio.h>  // printf, puts, putchar
#include <math.h>   // fabs

#define N   3
#define EPS 0.00001

void print_hline(void)
{
	puts("------------------------------------------------");
}

void print_array(double b[N][N * 2])
{
	int i, j;
	for (i = 0; i < N; i++) {
		for (j = 0; j < N * 2; j++)
			printf(" %7.3f", b[i][j]);
		putchar('\n');
	}
	print_hline();
}

void irekae(double b[N][N * 2], int i)
{
	int m = i, j, k;
	double work;

	for (k = i + 1; k < N; k++)
		if (fabs(b[m][i]) < fabs(b[k][i])) m = k;

	for (j = 0; j < N * 2; j++)
		work = b[m][j], b[m][j] = b[i][j], b[i][j] = work;
}

int main(void)
{
	double a[N][N] = {
		{  0.0, 2.0,  1.0 },
		{  4.0, 4.0,  1.0 },
		{ -2.0, 1.0, -2.0 },
	};
	double b[N][N * 2];
	double pivot, ratio;
	int i, j, k;
	for (i = 0; i < N; i++) {
		for (j = 0; j < N; j++) {
			b[i][j] = a[i][j];
			b[i][j + N] = 0.0;
		}
		b[i][i + N] = 1.0;
	}

	puts("前進");
	for (i = 0; i < N; i++) {
		pivot = b[i][i];
		if (fabs(pivot) < EPS) {
			irekae(b, i);
			pivot = b[i][i];
		}
		print_array(b);
		for (k = i + 1; k < N; k++) {
			ratio = b[k][i] / pivot;
			for (j = i; j < N * 2; j++) {
				b[k][j] -= b[i][j] * ratio;
			}
		}
	}

	puts("後退");
	for (i = N - 1; i >= 0; i--) {
		for (j = N - 1; j >= i + 1; j--) {
			for (k = N; k < N * 2; k++) {
				b[i][k] -= b[j][k] * b[i][j];
			}
			b[i][j] = 0.0;
		}
		for (k = N; k < N * 2; k++) {
			b[i][k] /= b[i][i];
		}
		b[i][i] = 1.0;
		print_array(b);
	}

	for (i = 0; i < N; i++) {
		for (j = N; j < 2 * N; j++)
			printf(" %7.3f", b[i][j]);
		printf("\n");
	}
#if 0  // #if 1 で、逆行列との積が単位行列になることの確認
	print_hline();
	for (i = 0; i < N; i++) {
		for (j = N; j < N * 2; j++) {
			double c = 0;
			for (k = 0; k < N; k++)
				c += a[i][k] * b[k][j];
			printf(" %7.3f", c);
		}
		putchar('\n');
	}
#endif
	return 0;
}
実行結果

コード: 

前進
   4.000   4.000   1.000   0.000   1.000   0.000
   0.000   2.000   1.000   1.000   0.000   0.000
  -2.000   1.000  -2.000   0.000   0.000   1.000
------------------------------------------------
   4.000   4.000   1.000   0.000   1.000   0.000
   0.000   2.000   1.000   1.000   0.000   0.000
   0.000   3.000  -1.500   0.000   0.500   1.000
------------------------------------------------
   4.000   4.000   1.000   0.000   1.000   0.000
   0.000   2.000   1.000   1.000   0.000   0.000
   0.000   0.000  -3.000  -1.500   0.500   1.000
------------------------------------------------
後退
   4.000   4.000   1.000   0.000   1.000   0.000
   0.000   2.000   1.000   1.000   0.000   0.000
   0.000   0.000   1.000   0.500  -0.167  -0.333
------------------------------------------------
   4.000   4.000   1.000   0.000   1.000   0.000
   0.000   1.000   0.000   0.250   0.083   0.167
   0.000   0.000   1.000   0.500  -0.167  -0.333
------------------------------------------------
   1.000   0.000   0.000  -0.375   0.208  -0.083
   0.000   1.000   0.000   0.250   0.083   0.167
   0.000   0.000   1.000   0.500  -0.167  -0.333
------------------------------------------------
  -0.375   0.208  -0.083
   0.250   0.083   0.167
   0.500  -0.167  -0.333
最初の「ガウス消去法」で解を求める過程を表示するプログラムは
どうなりますか？
ページトップ
返信

“C言語何でも質問掲示板” へ戻る