杜利特尔分解法Doolittle(LU分解法)_解线性方程组的直接解法
标签:计算方法实验
#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], u[maxn][maxn];
int i, j, k, r, n, sum;
freopen("lu.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(i = 1; i <= n; i++) u[1][i] = a[1][i]; //U的第一行元素
for(i = 1; i <= n; i++) l[i][1] = a[i][1] / u[1][1]; //L的第一列元素
for(k = 2; k <= n; k++){
for(j = k; j <= n; j++){ //计算U的第k行元素
for(r = 1, sum = 0; r <= k - 1; r++) sum += (l[k][r] * u[r][j]);
u[k][j] = a[k][j] - sum;
}
for(i = k; i <= n; i++){ //计算L的第k列元素
for(r = 1, sum = 0; r <= k - 1; r++) sum += (l[i][r] * u[r][k]);
l[i][k] = (a[i][k] - sum) / u[k][k];
}
}
/*打印L U
for(i = 1; i <= n; i++){
for(j = 1; j <= n; j++) printf("%10f", l[i][j]);
printf("\t\t");
for(k = 1; k <= n; k++) printf("%10f", u[i][k]);
printf("\n");
}
*/
y[1] = b[1]; //求解Ly = b
for(i = 2; i <= n; i++){
for(k = 1, sum = 0; k <= i - 1; k++) sum += (l[i][k] * y[k]);
y[i] = b[i] - sum;
}
x[n] = y[n] / u[n][n]; //求解Ux = y
for(i = n - 1; i >= 1; i--){
for(k = i + 1, sum = 0; k <= n; k++) sum += (u[i][k] * x[k]);
x[i] = (y[i] - sum) / u[i][i];
}
for(int i = 1; i <= n; i++) printf("%10f\n", x[i]);
return 0;
}
数据文件
实验结果