連立方程式を解くプログラムで、なぜベクトルで引き算を行うのかわかりません。

Question

#include <stdio.h>

#define N  3   //the dimension of equation

void vec_diff(float a[N], float b[N]){
    /* Calcurate the difference of two vectors. Be caution that b[N] changes.*/
    for (int i = 0; i < N; i++){
        b[i] -= a[i];
    }
}

int main(){
    float m[N][N] = {{5,-1,-1}, {2,1,-3},{1,1,1}};    // The matrix
    float b[N] = {0,-5,6};

    printf("The coefficient matrix is : \n");
    for (int i = 0; i < N; i++){
        for (int j = 0; j < N; j++){
            printf("%1.f ", m[i][j]);
            if (j == N-1){
                printf("\n");
            }
        }
    }

    printf("\nUse Gauss method to solve equations : \n");
    for (int i = 0; i < N; i++){
        for (int j = i+1; j < N; j++){
            float coef = m[j][i] / m[i][i];
            float del[N];

            for (int k = 0; k < N; k++){
                del[k] = m[i][k] * coef;
            }
            vec_diff(del, m[j]);
            b[j] -= b[i] * coef;
        }
    }

    for (int i = N -1; i >= 0; i--){
        float x = 1. / m[i][i];
        m[i][i] *= x;
        b[i] *= x;

        for (int j = 0; j < i; j++){
            b[j] -= b[i]*m[j][i];
            m[j][i] = 0;
        }
    }

    for (int i = 0; i < N; i++){
        for (int j = 0; j < N; j++){
            printf("%1.f ", m[i][j]);
            if (j == N - 1){
                printf("\n");
            }
        }
    }

    for (int i = 0; i < N; i++){
        printf("%f ", b[i]);
    }

    return 0;
}

のプログラムの

#define N  3   //the dimension of equation

void vec_diff(float a[N], float b[N]){
    /* Calcurate the difference of two vectors. Be caution that b[N] changes.*/
    for (int i = 0; i < N; i++){
        b[i] -= a[i];
    }
}

の部分に関してなのですが、
関数vec_diff(float a[N],float b[N])とb[i]-=a[i];よりベクトルの引き算をfor (int i = 0; i < N; i++)より３回行っていることはわかるのですが、「何のため何を３回」行っているのかわかりません。
プログラムを読んでみたのですが、意図した３点の座標を入れるための準備のようなものとしかわかりませんでした。

こちらのサイトを参考にしました。

＜編集７/６＞
あの後b[i] -= a[i];の a[i]とb[i]の計算を表示しよとしたのですが、うまくいきません。

#include <stdio.h>

#define N  3   //the dimension of equation

void vec_diff(float a[N], float b[N]) {
	/* Calcurate the difference of two vectors. Be caution that b[N] changes.*/
	for (int i = 0; i < N; i++) {
		b[i] -= a[i];
		printf("である",b[i] -= a[i]);//付け加えた部分、b[i] -= a[i];の計算を表示しようとしたのですが、これをいれると計算は表示されないことはおろか間違った答えが出てきます。
	}
}


int main() {
	float m[N][N] = { { 5,-1,-1 },{ 2,1,-3 },{ 1,1,1 } };    // The matrix
	float b[N] = { 0,-5,6 };

	printf("The coefficient matrix is : \n");
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N; j++) {
			printf("%1.f ", m[i][j]);
			if (j == N - 1) {
				printf("\n");
			}
		}
	}

	printf("\nUse Gauss method to solve equations : \n");
	for (int i = 0; i < N; i++) {
		for (int j = i + 1; j < N; j++) {
			float coef = m[j][i] / m[i][i];
			float del[N];

			for (int k = 0; k < N; k++) {
				del[k] = m[i][k] * coef;
			}
			vec_diff(del, m[j]);
			b[j] -= b[i] * coef;
		}
	}

	for (int i = N - 1; i >= 0; i--) {
		float x = 1. / m[i][i];
		m[i][i] *= x;
		b[i] *= x;

		for (int j = 0; j < i; j++) {
			b[j] -= b[i] * m[j][i];
			m[j][i] = 0;
		}
	}

	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N; j++) {
			printf("%1.f ", m[i][j]);
			if (j == N - 1) {
				printf("\n");
			}
		}
	}

	for (int i = 0; i < N; i++) {
		printf("%f です", b[i]);
	}

	return 0;
}

実行すると今のように出ます。

The coefficient matrix is :
5 -1 -1
2 1 -3
1 1 1

Use Gauss method to solve equations :
であるであるであるであるであるであるであるであるである1 0 0
-2 1 0
2 -1 1
0.355863 です-0.271378 です2.050691 です

Comments

[mkgrei]

プログラム以前に数式変形は把握していますか？3次元の連立方程式を書き下せていますか？

[carnage0216]

数式変形は把握しています。3次元の連立方程式は書き下せます。と言うことは、3は3次元のためのx,y,zを表しているのですか？

[mkgrei]

おっしゃる通り、3次元に対応して3回の引き算がそれぞれの要素に対して行われています。 3次元のベクトルを2つ渡して、その差を計算しています。in-placeの処理なので、返り値が必要ありません。ただ、値渡しが間違っているような…何か見落としたかな…https://qiita.com/agate-pris/items/05948b7d33f3e88b8967

[mkgrei]

これか。https://9cguide.appspot.com/15-06.html

[carnage0216]

解説ありがとうございます。あのなぜ引き算を行う必要があるのでしょうか？

[mkgrei]

下三角と上三角を0にするため？掃き出し法は理解してるんですよね？

[YomogiKOBO]

情報共有のため、他掲示板ＵＲＬを追記。https://dixq.net/forum/viewtopic.php?f=3&t=20251

Answer 1

「ふたつの等式: A = B と C = D が成り立つなら

• [1] tA = tB
• [2] A - C = B - D

も成り立つ」ことを利用して

• [3] b = mx (m:matrix, b,x:vector)

を解いている（ガウスの消去法,掃き出し法)。
操作[1],[2]により [3]を

• [4] b' = Ex (E:単位行列)

に変換できたとき、Ex = x だから b' = x すなわち b' が解となる。

[3]の右辺mxにおいてmの各行(とx)はベクトルであるから、
操作[2]はベクトルの引き算に相当する。

── 以上 ──

[ところで]
このコード、係数行列mの成分によってはコケる。
コード中、m[i][i]で割る処理が2か所に現れる。
したがって m[i][i]すなわち対角成分中に0があるとコケる。
これを回避するには必要に応じて行の入れ替え(pivot操作)を行わなければならない。

[さらに]
pivot操作によって対角成分を非0にできないとき、
それは逆行列が存在しないことを意味する。
すなわちこのとき連立方程式 b = mx は解を持たない。

Comments

[episteme]

>「何のため何を３回」行っているのかわかりません。

mの対角成分を1にするため/LU分解におけるUの左下部分を0にするため、操作[2]を繰り返している。
※ 行の各成分にcoeffを掛けてるのが 操作[1]に相当。

> 意図した３点の座標を入れるための準備のようなもの

チガウヨ。

[episteme]

> printf("である",b[i] -= a[i]);//付け加えた部分、b[i] -= a[i];の計算を表示しようとしたのですが、これをいれると計算は表示されないことはおろか間違った答えが出てきます。

プログラムは思ったとおりには動かない。作ったとおりに動く。
思ったとおりに動かないのは、思ったとおりに作ってないからだ。
printfに与えるフォーマットが間違っているし、同じ処理を二度やるから結果が狂う。

[carnage0216]

以下のように改良したところ結果は正しいものが出ました。
しかしまだ正しいb[i] -= a[i]が出ずにいます。printfのフォーマットはfloat型で実数を出力する%fにしたのですが、何がいけないのか…。
void vec_diff(float a[N], float b[N]) {

for (int i = 0; i < N; i++) {
printf("%f", b[i] -= a[i]);
}
}

実数が出力されるはずが、 b[i] -= a[i]においては
0.0000001.400000-2.6000000.0000001.2000001.2000000.0000000.0000003.4285711
と出ました。

[asm]

for(int i =0; i < N; i++){
printf("[%d]%f - %f =",i, b[i], a[i]);
printf("%f\n", b[i] -= a[i]);
}

なお、1回のprintfに纏めようとすると未定義動作となるので不可能です。
https://www.jpcert.or.jp/sc-rules/c-exp30-c.html

[carnage0216]

どうもありがとうございます。おかげさまでやっと理解できました。

Answer 2

#define N 3 //the dimension of equation
here 3 come in action

Answer 3

掃き出し法の各行の処理ですね。
掃き出しによって係数行列を単位行列にしようとしています。
まず、ガウス法、または、掃き出し法について調べる事をお勧めます。
この程度なら、手計算でもできると思いますので、その結果と比較してみたら、どうでしょうか。

---

なお、参考のサイトの「係数行列をガウス法で計算...」の最初の項の "3" ([2,2]の位置) は、"-3" の間違いでは無いかと思います。