C++ で連立一次方程式を解く(ガウス・ジョルダン法)
考え方は, ガウスの消去法とほぼ同じ.
連立一次方程式
が与えられた場合, まず(1)式を9で割って(1)'とする.
(2)式から(1)'式に2を掛けたものを引くとの項が消える.
(3)式に(1)'式を足すとの項が消える.
(4)式から(1)'式を引くとの項が消える.
ここまでは, ガウスの消去法と同じ.
ガウスの消去法は,
2番目の式を利用して, 3番目以降の式のを消す.
3番目の式を利用して, 4番目の式のを消す...
という方法だったが,
今回は,
2番目の式を利用して, 1番目の式と, 3番目以降の式のを消す.
3番目の式を利用して, 1番目, 2番目の式と, 4番目の式のを消す...
というようにして, 対角行列にする.
あとは, 各係数で割れば解が得られる.
#include <iostream> #include <iomanip> #include <math.h> using namespace std; void gauss(double a[4][4], double b[4], double c[4]); int main() { double a[4][4] = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}}; double b[4] = {8,17,20,16}; double c[4] = {0, 0, 0, 0}; gauss(a,b,c); for (int i = 0;i<=3;i++) cout << setw(14) << fixed << setprecision(10) << c[i] << "\t"; cout << endl; return 0; } void gauss(double a[4][4], double b[4], double c[4]) { // ピボット選択 int p[4] = {0,1,2,3}; for(int pivot = 0; pivot < 4; pivot++) { // 各列で 一番値が大きい行を 探す int col = pivot; int max_row = pivot; double max_val = 0; for (int i = p[pivot]; i < 4; i++) { int row = p[i]; if (fabs(a[row][col]) > max_val) { // 一番値が大きい行 max_val = fabs(a[row][col]); max_row = row; } } // 一番値が大きい行と入れ替え if (max_row != p[pivot]) { int tmp = p[pivot]; p[pivot] = max_row; p[max_row] = tmp; } } // 並べ替え後の状態を確認 for (int i = 0; i < 4; i++) { for (int j = 0; j < 4; j++) cout << setw(14) << fixed << setprecision(10) << a[p[i]][j] << "\t"; cout << setw(14) << fixed << setprecision(10) << b[p[i]] << "\t"; cout << endl; } cout << endl; // 対角線上の係数を 1 にする for (int pivot = 0; pivot < 4; pivot++) { int p_row = p[pivot]; int p_col = pivot; double s = a[p_row][p_col]; for (int col = 0; col < 4; col++) { a[p_row][col] /= s; } b[p_row] /= s; } for (int i = 0; i < 4; i++) { for (int j = 0; j < 4; j++) cout << setw(14) << fixed << setprecision(10) << a[p[i]][j] << "\t"; cout << setw(14) << fixed << setprecision(10) << b[p[i]] << "\t"; cout << endl; } cout << endl; // 対角行列にする for (int pivot = 0; pivot < 4; pivot++) { int p_row = p[pivot]; int p_col = pivot; for (int i = 0; i < 4; i++) { if (i == pivot) continue; int row = p[i]; double s = a[row][p_col] / a[p_row][p_col]; for (int col = pivot; col < 4; col++) a[row][col] -= a[p_row][col] * s; b[row] -= b[p_row] * s; } } // 状態を確認 for (int i = 0; i < 4; i++) { for (int j = 0; j < 4; j++) cout << setw(14) << fixed << setprecision(10) << a[p[i]][j] << "\t"; cout << setw(14) << fixed << setprecision(10) << b[p[i]] << "\t"; cout << endl; } cout << endl; // 単位行列にする for (int pivot = 0; pivot < 4; pivot++) { int p_row = p[pivot]; int p_col = pivot; b[p_row] /= a[p_row][p_col]; // 解 c[pivot] = b[p_row]; } }
9.0000000000 2.0000000000 1.0000000000 1.0000000000 20.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 16.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 8.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 17.0000000000 1.0000000000 0.2222222222 0.1111111111 0.1111111111 2.2222222222 0.2500000000 1.0000000000 -0.2500000000 0.1250000000 2.0000000000 -0.1428571429 -0.2857142857 1.0000000000 -0.2857142857 1.1428571429 0.1666666667 -0.1666666667 -0.3333333333 1.0000000000 2.8333333333 1.0000000000 0.0000000000 0.0000000000 0.0000000000 1.0000000000 0.0000000000 0.9444444444 0.0000000000 0.0000000000 1.8888888889 0.0000000000 0.0000000000 0.9411764706 -0.0000000000 2.8235294118 0.0000000000 0.0000000000 0.0000000000 0.8958333333 3.5833333333 1.0000000000 2.0000000000 3.0000000000 4.0000000000