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	

参考文献