C++ で連立一次方程式を解く(ガウスの消去法)
連立一次方程式
が与えられた場合, まず(1)式を9で割って(1)'とする.
(2)式から(1)'式に2を掛けたものを引くとの項が消える.
(3)式に(1)'式を足すとの項が消える.
(4)式から(1)'式を引くとの項が消える.
同様にして,
2番目の式を利用して, 3番目以降の式のyを消す.
3番目の式を利用して, 4番目の式のzを消す.
このようにして、上三角型方程式を得る.
(前進消去)
(4)式から の値が求まるので, (3)式に代入すると の値が求まる.
, の値を (2)式に代入すると の値が求まる.
, , の値を (1)式に代入すると の値が求まる.
(後退代入)
#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(12) << 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; // 前進消去 for (int pivot = 0; pivot < 3; pivot++) { int p_row = p[pivot]; int p_col = pivot; for (int i = pivot + 1; i < 4; i++) { 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 i = 3; i >= 0; i--) { int row = p[i]; for (int col = 3; col >= i + 1; col--) b[row] -= a[row][col] * c[col]; c[i] = b[row] / a[row][i]; } }
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 9.0000000000 2.0000000000 1.0000000000 1.0000000000 20.0000000000 0.0000000000 7.5555555556 -2.2222222222 0.7777777778 11.5555555556 0.0000000000 0.0000000000 6.5882352941 -1.7058823529 12.9411764706 0.0000000000 0.0000000000 0.0000000000 5.3750000000 21.5000000000 1.0000000000 2.0000000000 3.0000000000 4.0000000000