C++ - 連立方程式解法（ガウスの消去法（ピボット選択））！
Updated:
かつて、連立方程式を「ガウスの消去法」で解くアルゴリズムを C++ で実装したことを紹介しました。
しかし、計算途中で対角成分がゼロになるケースでは計算ができませんでした。
今回はその問題を解決すべく、「ガウスの消去法（ピボット選択）」で解くアルゴリズムを実装してみました。（少し前に作成していたが、ブログ記事として記録していなかった）
ちなみに、前々回、前回は Ruby, Fortran95 で実装しました。
0. 前提条件
- Debian GNU/Linux 11.3 での作業を想定。
- GCC 11.2.0 (G++ 11.2.0) (C++17) でのコンパイルを想定。
- 連立方程式の解法（ガウスの消去法（ピボット選択））についての説明は割愛。（Web 上等で容易に確認可能）
1. ソースコードの作成
File: gauss_elimination_pivot.cpp
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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
/***********************************************************
連立方程式の解法 ( ガウスの消去法（ピボット選択） )
DATE AUTHOR VERSION
2021.09.16 mk-mode.com 1.00 新規作成
Copyright(C) 2021 mk-mode.com All Rights Reserved.
***********************************************************/
#include <cmath> // for fabs
#include <cstdlib> // for EXIT_XXXX
#include <iostream>
#include <vector>
namespace gauss_elimination_pivot {
/*
* @brief ピボット選択
*
* @param[ref] 行列（配列） mtx (vector<vector<double>>)
* @param[in] 対象行
* @return true|false
*/
bool pivot(std::vector<std::vector<double>>& mtx, unsigned int k) {
unsigned int n; // 行数
unsigned int pv; // ピボット行
unsigned int i; // loop index
double v_max; // k 列絶対値の最大値
try {
n = mtx.size();
pv = k;
v_max = std::fabs(mtx[k][k]);
for (i = k + 1; i < n; ++i) {
if (std::fabs(mtx[i][k]) > v_max) {
pv = i;
v_max = std::fabs(mtx[i][k]);
}
}
if (k != pv) { swap(mtx[k], mtx[pv]); }
} catch (...) {
return false;
}
return true;
}
/*
* @brief ガウスの消去法
* * 値が NaN の場合は 0 とする
*
* @param[ref] 行列（配列） mtx (vector<vector<double>>)
* @return true|false
*/
bool solve_ge(std::vector<std::vector<double>>& mtx) {
int i; // loop インデックス
int j; // loop インデックス
int k; // loop インデックス
int n; // 元（行）の数
double d; // 計算用
try {
n = (int)mtx.size();
// 前進消去
for (k = 0; k < n - 1; ++k) {
if (!pivot(mtx, k)) {
std::cout << "[ERROR] Failed to pivot!" << std::endl;
return false;
}
for (i = k + 1; i < n; ++i) {
d = mtx[i][k] / mtx[k][k];
for (j = k + 1; j <= n; ++j)
mtx[i][j] -= mtx[k][j] * d;
}
}
// 後退代入
for (i = n - 1; i >= 0; --i) {
d = mtx[i][n];
for (j = i + 1; j < n; ++j)
d -= mtx[i][j] * mtx[j][n];
mtx[i][n] = d / mtx[i][i];
//if (std::isnan(mtx[i][n])) { mtx[i][n] = 0.0;}
}
} catch (...) {
return false;
}
return true;
}
} // namespace gauss_elimination_pivot
int main(int argc, char* argv[]) {
namespace ns = gauss_elimination_pivot;
std::vector<std::vector<double>> mtx; // 計算用行列
// テストデータ（右端は b, その他左側正方行列は a）
// * 3x3; 対角成分が 0 にならない例
//mtx.push_back({ 2.0, -3.0, 1.0, 5.0});
//mtx.push_back({ 1.0, 1.0, -1.0, 2.0});
//mtx.push_back({ 3.0, 5.0, -7.0, 0.0});
// * 4x4; 対角成分が 0 にならない例
//mtx.push_back({ 1.0, -2.0, 3.0, -4.0, 5.0});
//mtx.push_back({-2.0, 5.0, 8.0, -3.0, 9.0});
//mtx.push_back({ 5.0, 4.0, 7.0, 1.0, -1.0});
//mtx.push_back({ 9.0, 7.0, 3.0, 5.0, 4.0});
// * 4x4; 対角成分が 0 になる例
mtx.push_back({1.0, 2.0, 7.0, 6.0, 6.0});
mtx.push_back({2.0, 4.0, 4.0, 2.0, 2.0});
mtx.push_back({1.0, 8.0, 5.0, 2.0, 12.0});
mtx.push_back({2.0, 4.0, 3.0, 3.0, 5.0});
try {
// 元行列確認
for (auto &r: mtx) {
for (auto &c: r) { std::cout << c << " "; }
std::cout << std::endl;
}
// 計算（ガウスの消去法（ピボット選択））
if (!ns::solve_ge(mtx)) {
std::cout << "[ERROR] Failed to solve by the Gauss-Elimination(Pivot) method!"
<< std::endl;
return EXIT_SUCCESS;
}
// 計算結果確認
std::cout << "---" << std::endl;
for (auto &r: mtx) {
for (auto &c: r) { std::cout << c << " "; }
std::cout << std::endl;
}
} catch (...) {
std::cerr << "EXCEPTION!" << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
2. ソースコードのコンパイル
$ g++ -std=c++17 -Wall -O2 -o gauss_elimination_pivot gauss_elimination_pivot.cpp
3. 動作確認
$ ./gauss_elimination_pivot
1 2 7 6 6
2 4 4 2 2
1 8 5 2 12
2 4 3 3 5
---
2 4 4 2 -3
1 6 3 1 2
1 0 5 5 -1
2 0 -1 2 2
計算結果の最右列が解。
ちなみに、上の例は、ピボット選択しない場合には計算ができないケース。
試しにソースコードの 70行目〜73行目 をコメントアウトすると、以下のように計算できなくなる。
$ ./gauss_elimination_pivot
1 2 7 6 6
2 4 4 2 2
1 8 5 2 12
2 4 3 3 5
---
1 2 7 6 -nan
2 0 -10 -10 -nan
1 6 inf inf -nan
2 0 -nan -nan -nan
以上。
Comments