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 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