ONLY DO WHAT ONLY YOU CAN DO

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

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

行と列とを入れ替えても(転置行列)一致する行列を対称行列と言う.
f:id:fornext1119:20140728234117p:plain

コレスキー法では, 対称行列 A を下三角行列 L と, L の転置行列 L^T との積に分解し( A=LL^T ),
f:id:fornext1119:20140729223952p:plain

f:id:fornext1119:20140729223802p:plain

Ly=b から y を求め, L^Tx=y から x を求める.

#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 j = 0; j < 4; j++)
    {
        double s = 0;
        for (int k = 0; k < j; k++)
            s += a[j][k] * a[j][k];
        a[j][j] = sqrt(a[j][j] - s);

        for (int i = j + 1; i < 4; i++)
        {
            s = 0;
            for (int k = 0; k < j; k++)
                s += a[i][k] * a[j][k];
            a[i][j] =  (a[i][j] - s) / a[j][j];
            a[j][i] =  a[i][j];
        } 
    }
    // L と L^T
    for (int i = 0; i < 4; i++) 
    {
        for (int j = 0; j < 4; 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] / a[row][row];
    }

    // y の 値
    for (int i = 0;i<=3;i++)
        cout << setw(12) << fixed << setprecision(10) << y[i] << "\t";
    cout << endl << endl;

    // Ux=y から x を求める
    for (int row = 3; row >= 0; row--)
    {
        for (int col = 3; col >= row + 1; col--)
            y[row] -= a[row][col] * x[col];
        x[row] = y[row] / a[row][row];
    }
}
  2.2360679775	  0.8944271910	  1.3416407865	  1.7888543820	
  0.8944271910	  3.0331501776	  1.5825131361	  1.7803272782	
  1.3416407865	  1.5825131361	  3.2704207946	  1.1566122322	
  1.7888543820	  1.7803272782	  1.1566122322	  3.5060922587	

15.2052622470	17.9351488764	14.4377113129	14.0243690350	

1.0000000000	2.0000000000	3.0000000000	4.0000000000	
参考文献