平方根法 乔累斯基分解Cholesky_解线性方程组的直接解法
标签:计算方法实验
#include <stdio.h>
#include <math.h>
const int maxn = 15;
int main(){
double a[maxn][maxn], b[maxn], y[maxn], x[maxn], l[maxn][maxn];
int n, sum;
freopen("sqrt.txt", "r", stdin);
scanf("%d", &n);
for(int i = 1; i <= n; i++){
for(int j = 1; j <= n; j++) scanf("%lf", &a[i][j]);
scanf("%lf", &b[i]);
}
/*
for(int i = 1; i <= n; i++){
for(int j = 1; j <= n; j++) printf("%10f", a[i][j]);
printf("%10f\n", b[i]);
}
*/
for(int i = 1; i <= n; i++){ //算L: 对正定矩阵A分解,A = LL^T
sum = 0;
for(int k = 1; k <= i - 1; k++) sum += (l[i][k] * l[i][k]);
l[i][i] = sqrt(a[i][i] - sum);
sum = 0;
for(int j = i + 1; j <= n; j++){
for(int k = 1; k <= i - 1; k++) sum += (l[j][k] * l[i][k]);
l[j][i] = (a[j][i] - sum) / l[i][i];
}
}
for(int i = 1; i <= n; i++){ //求y: A = LL^T -> LL^Tx = b -> Ly = b
sum = 0;
for(int k = 1; k <= i - 1; k++) sum += (l[i][k] * y[k]);
y[i] = (b[i] - sum) / l[i][i];
}
for(int i = n; i >= 1; i--){ //求x: L^Tx = y
sum = 0;
for(int k = i + 1; k <= n; k++) sum += (l[k][i] * x[k]);
x[i] = (y[i] - sum) / l[i][i];
}
for(int i = 1; i <= n; i++) printf("%10f\n", x[i]);
return 0;
}
数据文件
实验结果