C言語で連立一次方程式の解き方の解き方を教えて！
こういった悩みにお答えします．
本記事の信頼性
- リアルタイムシステムの研究歴12年．
- 東大教員の時に，英語でOS（Linuxカーネル）の授業．
- 2012年9月～2013年8月にアメリカのノースカロライナ大学チャペルヒル校（UNC）コンピュータサイエンス学部で客員研究員として勤務．C言語でリアルタイムLinuxの研究開発．
- プログラミング歴15年以上，習得している言語: C/C++，Python，Solidity/Vyper，Java，Ruby，Go，Rust，D，HTML/CSS/JS/PHP，MATLAB，Assembler (x64，ARM)．
- 東大教員の時に，C++言語で開発した「LLVMコンパイラの拡張」，C言語で開発した独自のリアルタイムOS「Mcube Kernel」をGitHubにオープンソースとして公開．
- 2020年1月～現在はアメリカのノースカロライナ州チャペルヒルにあるGuarantee Happiness LLCのCTOとしてECサイト開発やWeb/SNSマーケティングの業務．2022年6月～現在はアメリカのノースカロライナ州チャペルヒルにあるJapanese Tar Heel, Inc.のCEO兼CTO．
- 最近は自然言語処理AIとイーサリアムに関する有益な情報発信に従事．
- （AI全般を含む）自然言語処理AIの論文の日本語訳や，AIチャットボット（ChatGPT，Auto-GPT，Bardなど）の記事を30本以上執筆．アメリカのサンフランシスコ（広義のシリコンバレー）の会社でプロンプトエンジニア・マネージャーの業務委託の経験あり．
- （スマートコントラクトのプログラミングを含む）イーサリアムや仮想通貨全般の記事を250本以上執筆．イギリスのロンドンの会社で仮想通貨の英語の記事を日本語に翻訳する業務委託の経験あり．
こういった私から学べます．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
本記事は行列を習得していることを前提とします．
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でなければならない場合もある．）
ガウス・ジョルダンの消去法のコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 |
/* * 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; } |
実行結果は以下になります．
こちらと同じ結果になることを確認しましょう．
|
1 2 3 4 5 |
$ gcc gauss_jordan.c $ a.out x1 = 1.000000 x2 = 0.000000 x3 = -1.000000 |
クラメルの公式（クラメルの法則）
クラメルの公式（クラメルの法則）は，未知数の数と方程式の数が同じで，かつ一意的に解ける線型方程式系を明示的に解く行列式の公式です．
クラメルの公式は逆行列が存在しない場合（det(A)が0の場合）は利用できないことに注意して下さい．
クラメルの公式のコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 |
/* * 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; } |
実行結果は以下になります．
|
1 2 3 4 5 |
$ gcc cramer.c $ a.out x1 = 1.000000 x2 = 0.000000 x3 = -1.000000 |
ガウス・ザイデル法
ガウス・ザイデル法は，連立一次方程式を反復法で解くアルゴリズムです．
ガウス・ザイデル法は，行列の対角要素の絶対値が非対角要素の絶対値より相対的に小さい場合，反復が収束せず解を得られないことに注意して下さい．
ガウス・ザイデル法のコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 |
/* * 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; } |
実行結果は以下になります．
|
1 2 3 4 5 |
$ gcc gauss_seidel.c $ a.out x1 = 1.000000 x2 = -0.000000 x3 = -1.000000 |
ヤコビ法
ヤコビ法は，連立一次方程式を反復法で解くアルゴリズムです．
ガウス・ザイデル法と同様に，ヤコビ法は，行列の対角要素の絶対値が非対角要素の絶対値より相対的に小さい場合，反復が収束せず解を得られないことに注意して下さい．
ヤコビ法は，直列計算ではガウス・ザイデル法よりも遅いデメリットがありますが，ガウス・ザイデル法と異なり並列計算による高速化ができるメリットがあります．
ヤコビ法のコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 |
/* * 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; } |
実行結果は以下になります．
|
1 2 3 4 5 |
$ gcc jacobi.c $ a.out x1 = 1.000000 x2 = 0.000000 x3 = -1.000000 |
まとめ
C言語で連立一次方程式の解法として，ガウスの消去法（掃き出し法），クラメルの公式（クラメルの法則），ガウス・ザイデル法，ヤコビ法を紹介しました．
どれも良く使うので，きちんと理解しておきましょう！
方程式の解の公式を学びたいあなたはこちらからどうぞ．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
