【C言語】連立一次方程式の解法【ガウスの消去法（掃き出し法），クラメルの公式（クラメルの法則），ガウス・ザイデル法，ヤコビ法】

悩んでいる人

C言語で連立一次方程式の解き方の解き方を教えて！

こういった悩みにお答えします．

本記事の信頼性

リアルタイムシステムの研究歴12年．
東大教員の時に，英語でOSの授業．
2012年9月～2013年8月にアメリカのノースカロライナ大学チャペルヒル校コンピュータサイエンス学部（2021年の世界大学学術ランキングで20位）で客員研究員として勤務．C言語でリアルタイムLinuxの研究開発．
プログラミング歴15年以上，習得している言語: C/C++，Solidity/Vyper，Java，Python，Ruby，HTML/CSS/JS/PHP，MATLAB，Assembler (x64，ARM)．
東大教員の時に，C++言語で開発した「LLVMコンパイラの拡張」，C言語で開発した独自のリアルタイムOS「Mcube Kernel」をGitHubにオープンソースとして公開．

こういった私から学べます．

本記事は行列を習得していることを前提とします．

: 【C言語】2次元配列で行列の四則演算（足し算，引き算，掛け算，割り算）【逆行列】

こういった悩みにお答えします．こういった私から学べます．目次1 行列の四則演算（足し算，引き算，掛け算，割り算）2 行列の足し算と引き算3 行列の掛け算4 行列の割り算（逆行列の掛け算）4.1 2 ...

続きを見る

1 C言語で連立一次方程式の解法
2 ガウスの消去法（掃き出し法）
3 クラメルの公式（クラメルの法則）
4 ガウス・ザイデル法
5 ヤコビ法
6 まとめ

C言語で連立一次方程式の解法

C言語で連立一次方程式の解法を紹介します．

具体的には，以下の解法を紹介します．

ガウスの消去法（掃き出し法）
クラメルの公式（クラメルの法則）
ガウス・ザイデル法
ヤコビ法

本記事では，以下の連立一次方程式を解きます．

\begin{equation} \left\{ \, \begin{aligned} 4x_1 + x_2 - 2x_3 &= 6 \\ x_1 + 6x_2 + 3x_3 &= -2 \\ 2x_1 + x_2 + 9x_3 &= -7 \end{aligned} \right. \end{equation}

ガウスの消去法（掃き出し法）

ガウスの消去法（掃き出し法）は，連立一次方程式を解くための多項式時間アルゴリズムです．

特に，行基本変形を用いて行列を行階段形に変形することをガウス・ジョルダンの消去法と呼びます．

行列が行階段形であるとは，以下の2つが成立することを意味します．

0でない成分を持つ行（少なくとも1つの成分が0でない行）が，0しか成分に持たない行よりも上に位置すること（もし0成分だけからなる行が存在する場合，それらは行列の最下部に配置）
主成分（行の最も左にある0でない成分で，ピボットとも呼ばれる）が，その行の上にある行の主成分よりも真に右側に位置すること（主成分は必ず1でなければならない場合もある．）

ガウス・ジョルダンの消去法のコードは以下になります．

/*
 * Author: Hiroyuki Chishiro
 * License: 2-Clause BSD
 */
#include <stdio.h>
#define N 3

int main(void)
{
  double a[N][N + 1] = {
    {4.0, 1.0, -2.0, 6.0},
    {1.0, 6.0, 3.0, -2.0},
    {2.0, 1.0, 9.0, -7.0}
  };
  double p, d;
  int i, j, k;

  for (i = 0; i < N; i++) {
    p = a[i][i];

    for (j = 0; j < (N + 1); j++) {
      a[i][j] = a[i][j] / p;
    }

    for (j = 0; j < N; j++) {
      if (i != j) {
        d = a[j][i];

        for (k = i; k < (N + 1); k++) {
          a[j][k] = a[j][k] - d * a[i][k];
        }
      }
    }
  }

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

  return 0;
}

* Author: Hiroyuki Chishiro

* License: 2-Clause BSD

#include <stdio.h>

#define N 3

int main(void)

{

double a[N][N + 1] = {

{4.0, 1.0, -2.0, 6.0},

{1.0, 6.0, 3.0, -2.0},

{2.0, 1.0, 9.0, -7.0}

};

double p, d;

int i, j, k;

for (i = 0; i < N; i++) {

p = a[i][i];

for (j = 0; j < (N + 1); j++) {

a[i][j] = a[i][j] / p;

}

for (j = 0; j < N; j++) {

if (i != j) {

d = a[j][i];

for (k = i; k < (N + 1); k++) {

a[j][k] = a[j][k] - d * a[i][k];

}

for (i = 0; i < N; i++) {

printf("x%d = %f\n", i + 1, a[i][N]);

}

return 0;

}

実行結果は以下になります．

こちらと同じ結果になることを確認しましょう．

$ gcc gauss_jordan.c 
$ a.out 
x1 = 1.000000
x2 = 0.000000
x3 = -1.000000

$ gcc gauss_jordan.c

$ a.out

x1 = 1.000000

x2 = 0.000000

x3 = -1.000000

クラメルの公式（クラメルの法則）

クラメルの公式（クラメルの法則）は，未知数の数と方程式の数が同じで，かつ一意的に解ける線型方程式系を明示的に解く行列式の公式です．

クラメルの公式は逆行列が存在しない場合（det(A)が0の場合）は利用できないことに注意して下さい．

クラメルの公式のコードは以下になります．

/*
 * Author: Hiroyuki Chishiro
 * License: 2-Clause BSD
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 3
#define DBL_EPSILON 2.2204460492503131e-16

double get_det(double a[N][N])
{
  double det = a[0][0] * a[1][1] * a[2][2]
               + a[0][1] * a[1][2] * a[2][0]
               + a[0][2] * a[1][0] * a[2][1]
               - a[0][2] * a[1][1] * a[2][0]
               - a[0][1] * a[1][0] * a[2][2]
               - a[0][0] * a[1][2] * a[2][1];

  return det;
}

int main(void)
{
  double a[N][N + 1] = {
    {4.0, 1.0, -2.0, 6.0},
    {1.0, 6.0, 3.0, -2.0},
    {2.0, 1.0, 9.0, -7.0}
  };
  double det = a[0][0] * a[1][1] * a[2][2]
               + a[0][1] * a[1][2] * a[2][0]
               + a[0][2] * a[1][0] * a[2][1]
               - a[0][2] * a[1][1] * a[2][0]
               - a[0][1] * a[1][0] * a[2][2]
               - a[0][0] * a[1][2] * a[2][1];
  double b[N][N];
  double x[N];
  int i, j, k;

  if (fabs(det) <= DBL_EPSILON) {
    fprintf(stderr, "Error: cannot find inverse matrix.\n");
    exit(1);
  }

  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      for (k = 0; k < N; k++) {
        if (i == k) {
          b[j][k] = a[j][N];
        } else {
          b[j][k] = a[j][k];
        }
      }
    }

    x[i] = get_det(b) / det;
  }

  for (i = 0; i < N; i++) {
    printf("x%d = %f\n", i + 1, x[i]);
  }

  return 0;
}

* Author: Hiroyuki Chishiro

* License: 2-Clause BSD

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#define N 3

#define DBL_EPSILON 2.2204460492503131e-16

double get_det(double a[N][N])

{

double det = a[0][0] * a[1][1] * a[2][2]

+ a[0][1] * a[1][2] * a[2][0]

+ a[0][2] * a[1][0] * a[2][1]

- a[0][2] * a[1][1] * a[2][0]

- a[0][1] * a[1][0] * a[2][2]

- a[0][0] * a[1][2] * a[2][1];

return det;

}

int main(void)

{

double a[N][N + 1] = {

{4.0, 1.0, -2.0, 6.0},

{1.0, 6.0, 3.0, -2.0},

{2.0, 1.0, 9.0, -7.0}

};

double det = a[0][0] * a[1][1] * a[2][2]

+ a[0][1] * a[1][2] * a[2][0]

+ a[0][2] * a[1][0] * a[2][1]

- a[0][2] * a[1][1] * a[2][0]

- a[0][1] * a[1][0] * a[2][2]

- a[0][0] * a[1][2] * a[2][1];

double b[N][N];

double x[N];

int i, j, k;

if (fabs(det) <= DBL_EPSILON) {

fprintf(stderr, "Error: cannot find inverse matrix.\n");

exit(1);

}

for (i = 0; i < N; i++) {

for (j = 0; j < N; j++) {

for (k = 0; k < N; k++) {

if (i == k) {

b[j][k] = a[j][N];

} else {

b[j][k] = a[j][k];

}

x[i] = get_det(b) / det;

}

for (i = 0; i < N; i++) {

printf("x%d = %f\n", i + 1, x[i]);

}

return 0;

}

実行結果は以下になります．

$ gcc cramer.c 
$ a.out 
x1 = 1.000000
x2 = 0.000000
x3 = -1.000000

$ gcc cramer.c

$ a.out

x1 = 1.000000

x2 = 0.000000

x3 = -1.000000

ガウス・ザイデル法

ガウス・ザイデル法は，連立一次方程式を反復法で解くアルゴリズムです．

ガウス・ザイデル法は，行列の対角要素の絶対値が非対角要素の絶対値より相対的に小さい場合，反復が収束せず解を得られないことに注意して下さい．

ガウス・ザイデル法のコードは以下になります．

/*
 * Author: Hiroyuki Chishiro
 * License: 2-Clause BSD
 */
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>

#define TOLERABLE_ERROR 1.0e-10
#define N 3

bool is_more_than_tolerable_error(double e[N])
{
  int i;

  for (i = 0; i < N; i++) {
    if (e[i] <= TOLERABLE_ERROR) {
      return false;
    }
  }

  return true;
}

int main(void)
{
  double x[N] = {0.0, 0.0, 0.0};
  double e[N] = {0.0, 0.0, 0.0};
  int i, j;
  double a[N][N + 1] = {
    {4.0, 1.0, -2.0, 6.0},
    {1.0, 6.0, 3.0, -2.0},
    {2.0, 1.0, 9.0, -7.0}
  };
  double sum, val;

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

      for (j = 0; j < N; j++) {
        if (i != j) {
          sum += a[i][j] * x[j];
        }
      }

      val = (1.0 / a[i][i]) * (a[i][N] - sum);
      e[i] = fabs(val - x[i]);
      x[i] = val;
    }
  } while (is_more_than_tolerable_error(e));

  for (i = 0; i < N; i++) {
    printf("x%d = %f\n", i + 1, x[i]);
  }

  return 0;
}

* Author: Hiroyuki Chishiro

* License: 2-Clause BSD

#include <stdio.h>

#include <stdlib.h>

#include <stdbool.h>

#include <math.h>

#define TOLERABLE_ERROR 1.0e-10

#define N 3

bool is_more_than_tolerable_error(double e[N])

{

int i;

for (i = 0; i < N; i++) {

if (e[i] <= TOLERABLE_ERROR) {

return false;

}

return true;

}

int main(void)

{

double x[N] = {0.0, 0.0, 0.0};

double e[N] = {0.0, 0.0, 0.0};

int i, j;

double a[N][N + 1] = {

{4.0, 1.0, -2.0, 6.0},

{1.0, 6.0, 3.0, -2.0},

{2.0, 1.0, 9.0, -7.0}

};

double sum, val;

do {

for (i = 0; i < N; i++) {

sum = 0.0;

for (j = 0; j < N; j++) {

if (i != j) {

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

}

val = (1.0 / a[i][i]) * (a[i][N] - sum);

e[i] = fabs(val - x[i]);

x[i] = val;

}

} while (is_more_than_tolerable_error(e));

for (i = 0; i < N; i++) {

printf("x%d = %f\n", i + 1, x[i]);

}

return 0;

}

実行結果は以下になります．

$ gcc gauss_seidel.c
$ a.out 
x1 = 1.000000
x2 = -0.000000
x3 = -1.000000

$ gcc gauss_seidel.c

$ a.out

x1 = 1.000000

x2 = -0.000000

x3 = -1.000000

ヤコビ法

ヤコビ法は，連立一次方程式を反復法で解くアルゴリズムです．

ガウス・ザイデル法と同様に，ヤコビ法は，行列の対角要素の絶対値が非対角要素の絶対値より相対的に小さい場合，反復が収束せず解を得られないことに注意して下さい．

ヤコビ法は，直列計算ではガウス・ザイデル法よりも遅いデメリットがありますが，ガウス・ザイデル法と異なり並列計算による高速化ができるメリットがあります．

ヤコビ法のコードは以下になります．

/*
 * Author: Hiroyuki Chishiro
 * License: 2-Clause BSD
 */
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>

#define TOLERABLE_ERROR 1.0e-10
#define N 3

double jacobi(int i, double aa[N], double v[N])
{
  if (i == 0) {
    return (aa[3] - aa[1] * v[1] - aa[2] * v[2]) / aa[0];
  } else if (i == 1) {
    return (aa[3] - aa[0] * v[0] - aa[2] * v[2]) / aa[1];
  } else if (i == 2) {
    return (aa[3] - aa[0] * v[0] - aa[1] * v[1]) / aa[2];
  } else {
    fprintf(stderr, "Error: unknown index %d\n", i);
    exit(1);
  }
}

bool is_more_than_tolerable_error(double e[N])
{
  int i;

  for (i = 0; i < N; i++) {
    if (e[i] <= TOLERABLE_ERROR) {
      return false;
    }
  }

  return true;
}

int main(void)
{
  double v[N] = {0.0, 0.0, 0.0};
  double x[N], e[N];
  int i;
  double a[N][N + 1] = {
    {4.0, 1.0, -2.0, 6.0},
    {1.0, 6.0, 3.0, -2.0},
    {2.0, 1.0, 9.0, -7.0}
  };

  do {
    for (i = 0; i < N; i++) {
      x[i] = jacobi(i, a[i], v);
    }

    for (i = 0; i < N; i++) {
      e[i] = fabs(v[i] - x[i]);
      v[i] = x[i];
    }
  } while (is_more_than_tolerable_error(e));

  for (i = 0; i < N; i++) {
    printf("x%d = %f\n", i + 1, x[i]);
  }

  return 0;
}

* Author: Hiroyuki Chishiro

* License: 2-Clause BSD

#include <stdio.h>

#include <stdlib.h>

#include <stdbool.h>

#include <math.h>

#define TOLERABLE_ERROR 1.0e-10

#define N 3

double jacobi(int i, double aa[N], double v[N])

{

if (i == 0) {

return (aa[3] - aa[1] * v[1] - aa[2] * v[2]) / aa[0];

} else if (i == 1) {

return (aa[3] - aa[0] * v[0] - aa[2] * v[2]) / aa[1];

} else if (i == 2) {

return (aa[3] - aa[0] * v[0] - aa[1] * v[1]) / aa[2];

} else {

fprintf(stderr, "Error: unknown index %d\n", i);

exit(1);

}

bool is_more_than_tolerable_error(double e[N])

{

int i;

for (i = 0; i < N; i++) {

if (e[i] <= TOLERABLE_ERROR) {

return false;

}

return true;

}

int main(void)

{

double v[N] = {0.0, 0.0, 0.0};

double x[N], e[N];

int i;

double a[N][N + 1] = {

{4.0, 1.0, -2.0, 6.0},

{1.0, 6.0, 3.0, -2.0},

{2.0, 1.0, 9.0, -7.0}

};

do {

for (i = 0; i < N; i++) {

x[i] = jacobi(i, a[i], v);

}

for (i = 0; i < N; i++) {

e[i] = fabs(v[i] - x[i]);

v[i] = x[i];

}

} while (is_more_than_tolerable_error(e));

for (i = 0; i < N; i++) {

printf("x%d = %f\n", i + 1, x[i]);

}

return 0;

}

実行結果は以下になります．

$ gcc jacobi.c 
$ a.out 
x1 = 1.000000
x2 = 0.000000
x3 = -1.000000

$ gcc jacobi.c

$ a.out

x1 = 1.000000

x2 = 0.000000

x3 = -1.000000

まとめ

C言語で連立一次方程式の解法として，ガウスの消去法（掃き出し法），クラメルの公式（クラメルの法則），ガウス・ザイデル法，ヤコビ法を紹介しました．

どれも良く使うので，きちんと理解しておきましょう！

方程式の解の公式を学びたいあなたはこちらからどうぞ．

: 【C言語】一次方程式，二次方程式，三次方程式，四次方程式の解の公式

こういった悩みにお答えします．こういった私から学べます．目次1 C言語で一次方程式，二次方程式，三次方程式，四次方程式の解の公式2 一次方程式の解の公式3 ニ次方程式の解の公式3.1 ニ次方程式の ...

続きを見る

C言語を独学で習得することは難しいです．

私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．

私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！

独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！

: 元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社

こういった悩みにお答えします．本記事の内容：おすすめのC言語を学べるオンラインプログラミングスクール5社 SAMURAI ENGINEER（侍エンジニア） Tech Teacher KENスクール ...

続きを見る