C++实现Cholesky分解-创新互联
题目:
编制程序求解矩阵 A 的 Cholesky 分解,并用程序求解方程组 Ax=b,其中
代码实现:
#include#include#includeusing namespace std;
//Cholesky 分解
void Cholesky(double a[4][4], double L[4][4], int m, int n) {for (int j = 0; j< n; j++)
for (int i = j; i< m; i++) {if (i == j) { double sum = 0;
for (int k = 0; k< j; k++) { sum = sum + L[j][k] * L[j][k];
}
L[i][j] = sqrt(a[i][j] - sum);
}
else { double sum = 0;
for (int k = 0; k< j; k++) {sum = sum + L[i][k] * L[j][k];
}
L[i][j] = (a[i][j] - sum) / L[j][j];
}
}
}
//初始化矩阵 L
void InitL(double L[4][4], int m, int n) {for (int i = 0; i< m; i++)
for (int j = 0; j< n; j++) { if (i< j) { L[i][j] = 0;
}
}
}
//显示 L 和 LT 矩阵
void Display(double L[4][4], double LT[4][4], int m, int n) {cout<< "Cholesky 分解的 L 矩阵为:"<< endl;
for (int i = 0; i< m; i++) {for (int j = 0; j< n; j++) { cout<< setprecision(8)<< setw(12)<< L[i][j]<< "";
}
cout<< endl;
}
cout<< "Cholesky 分解的 LT 矩阵为:"<< endl;
for (int i = 0; i< m; i++) {for (int j = 0; j< n; j++) { cout<< setprecision(8)<< setw(12)<< L[j][i]<< "";
LT[i][j] = L[j][i];
}
cout<< endl;
}
}
//Ax=b 的解答 中间步骤求 y
double* SolveOne(double L[4][4], double b[4], int m, int n) {static double y[4];
y[0] = b[0] / L[0][0];
for (int i = 1; i< m; i++) {double sum = 0;
for (int k = 0; k< i; k++) { sum = sum + L[i][k] * y[k];
}
y[i] = (b[i] - sum) / L[i][i];
}
cout<< "解答的中间结果 y 为:"<< endl;
for (int i = 0; i< m; i++) {cout<< setprecision(8)<< setw(12)<< y[i]<< endl;
}
cout<< endl;
return y;
}
//Ax=b 的解答 最终结果
void SolveTwo(double L[4][4], double y[4], int m, int n) {double x[4];
x[m - 1] = y[m - 1] / L[n - 1][n - 1];
for (int i = m - 2; i >= 0; i--) {double sum = 0;
for (int k = i + 1; k< m; k++) { sum = sum + L[k][i] * x[k];
}
x[i] = (y[i] - sum) / L[i][i];
}
cout<< "解答的最终结果 x 为:"<< endl;
for (int i = 0; i< m; i++) {cout<< setprecision(8)<< setw(12)<< x[i]<< endl;
}
cout<< endl;
}
//主函数
int main() {double a[4][4] = {{2,1,-1,1},{1,5,2,7},{-1,2,10,1},{1,7,1,11} };
double b[4] = {13,-9,6,0 };
double L[4][4];//L 矩阵
double LT[4][4];//LT 矩阵
double* y;//中间矩阵 y
InitL(L, 4, 4);//初始化矩阵 L
Cholesky(a, L, 4, 4);//Cholesky 分解矩阵 A
Display(L, LT, 4, 4);//显示矩阵 L 和 LT
//求解原方程 Ax=b 中间步骤:L*y = b; LT*x = y;
y = SolveOne(L, b, 4, 4);
SolveTwo(L, y, 4, 4);
return 0;
}
数值结果:
分析:
首先对 A 矩阵进行 Cholesky 分解,得到 L 矩阵,然后通过中间矩阵 y 求解矩阵 x。
你是否还在寻找稳定的海外服务器提供商?创新互联www.cdcxhl.cn海外机房具备T级流量清洗系统配攻击溯源,准确流量调度确保服务器高可用性,企业级服务器适合批量采购,新人活动首月15元起,快前往官网查看详情吧
网站题目:C++实现Cholesky分解-创新互联
分享路径:http://pwwzsj.com/article/deeojc.html