C++ で連立一次方程式を解く(修正コレスキー法)

行と列とを入れ替えても(転置行列)一致する行列を対称行列と言う.

修正コレスキー法では, 対称行列 を下三角行列 , 対角行列 , の転置行列 の積に分解し( ),

から を求め, から を求める.

#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

void decomp(double a[4][4], double b[4], double x[4]);

int main()
{
    double a[4][4] = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}}; 
    double b[4]    = {34,68,96,125};
    double x[4]    = {0,0,0,0};

    decomp(a,b,x);
    for (int i = 0;i<=3;i++)
        cout << setw(12) << fixed << setprecision(10) << x[i] << "\t";
    cout << endl;
   
    return 0;
}

void decomp(double a[4][4], double b[4], double x[4])
{
    for(int i = 0; i < 4; i++)
    {
        double s;

        // i < k の場合
        for (int j = 0; j < i; j++)
        {
            s = a[i][j];
            for (int k = 0; k < j; k++)
                s -= a[i][k] * a[j][k] * a[k][k];
            a[i][j] = s / a[j][j];
        }
 
        // i == k の場合
        s = a[i][i];
        for (int k = 0; k < i; k++)
            s -= a[i][k] * a[i][k] * a[k][k];
        a[i][i] = s;
    }
    // L と D
    for (int i = 0; i < 4; i++) 
    {
        for (int j = 0; j <= i; j++) 
            cout << setw(14) << fixed << setprecision(10) << a[i][j] << "\t";
        cout << endl;
    }
    cout << endl;

    // Ly=b から y を求める
    double y[4] = {0,0,0,0};
    for (int row = 0; row < 4; row++)
    {
        for (int col = 0; col < row; col++)
            b[row] -= a[row][col] * y[col];
        y[row] = b[row];
    }
    // y の 値
    for (int i = 0; i < 4; i++)
        cout << setw(12) << fixed << setprecision(10) << y[i] << "\t";
    cout << endl << endl;

    // DL^Tx=y から x を求める
    for (int row = 3; row >= 0; row--)
    {
        for (int col = 3; col >= row + 1; col--)
            y[row] -= a[col][row] * a[row][row] * x[col];
        x[row] = y[row] / a[row][row];
    }
}

  5.0000000000	
  0.4000000000	  9.2000000000	
  0.6000000000	  0.5217391304	 10.6956521739	
  0.8000000000	  0.5869565217	  0.3536585366	 12.2926829268	

34.0000000000	54.4000000000	47.2173913043	49.1707317073	

1.0000000000	2.0000000000	3.0000000000	4.0000000000	

参考文献