ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ 立ったら歩きなはれ

C++ で連立一次方程式を解く(ガウス・ジョルダン法)

考え方は, ガウスの消去法とほぼ同じ.

連立一次方程式
9x+2y+z+u=20
2x+8y-2z+u=16
-x-2y+7z-2u=8
x-y-2z+6u=17
が与えられた場合, まず(1)式を9で割って(1)'とする.
x+\frac{2}{9}y+\frac{1}{9}z+\frac{1}{9}u=\frac{20}{9}

(2)式から(1)'式に2を掛けたものを引くとxの項が消える.
(3)式に(1)'式を足すとxの項が消える.
(4)式から(1)'式を引くとxの項が消える.

ここまでは, ガウスの消去法と同じ.

ガウスの消去法は,
2番目の式を利用して, 3番目以降の式のyを消す.
3番目の式を利用して, 4番目の式のzを消す...
という方法だったが,
今回は,
2番目の式を利用して, 1番目の式と, 3番目以降の式のyを消す.
3番目の式を利用して, 1番目, 2番目の式と, 4番目の式のzを消す...
というようにして, 対角行列にする.
あとは, 各係数で割れば解が得られる.

#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	
参考文献