版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/Richard__Ting/article/details/84949372
数值分析实验报告 Lab5 LU分解算法
#include<stdio.h>
#include<iostream>
#include<math.h>
#define MAX_SIZE 100 /* 矩阵最大维数 */
#define ZERO 0.000000001 /* 当一个正数小于ZERO就认为该数是0 */
using namespace std;
/*
LU方法计算 Ax=b
*/
bool Direct( int n, double a[][MAX_SIZE], double b[] ){
/*初始化矩阵,向量*/
double L[n][n];
double U[n][n];
for(int i = 0;i < n;i++){
for(int j = 0;j < n;j++){
L[i][j] = 0;
if(i == j)L[i][j] = 1;
U[i][j] = 0;
}
}
double y[n];
for(int i = 0;i < n;i++){
y[i] = 0;
}
for(int k = 0;k < n;k++){
for(int j1 = k;j1 < n;j1++){
double sum_U = 0;
for(int x1 = 0;x1 < k;x1++){
sum_U += L[k][x1]*U[x1][j1];
}
U[k][j1] = a[k][j1] - sum_U;
//cout<<U[k][j1]<<endl;
if((j1==k)&&(fabs(U[k][j1])<ZERO)){return 0;}
}
//cout<<endl;
for(int y2 = k+1;y2 < n;y2++){
double sum_L = 0;
for(int x2 = 0;x2 < k;x2++){
sum_L += L[y2][x2]*U[x2][k];
}
L[y2][k] = (a[y2][k] - sum_L)/U[k][k];
//cout<<L[y2][k]<<endl;
}
//cout<<endl;
}
//计算Ly=b中y的值
for(int m1 = 0;m1 < n;m1++){
double tmp1 = 0;
for(int t1 = 0;t1 < m1;t1++){
tmp1 += y[t1]*L[m1][t1];
}
y[m1] = b[m1] - tmp1;
//cout<<y[m1]<<endl;
}
//计算Ux=y中x的值
b[n-1] = y[n-1]/U[n-1][n-1];
//cout<<b[n-1]<<endl;
for(int m2 = n-2;m2 >= 0;m2--){
double tmp2 = 0;
for(int t2 = m2+1;t2 < n;t2++){
tmp2 += U[m2][t2]*b[t2];
}
b[m2] = (y[m2] - tmp2)/U[m2][m2];
}
return b;
}
int main(){
int n, i, j;
double a[MAX_SIZE][MAX_SIZE], b[MAX_SIZE];
freopen("input.txt", "r", stdin);
while ( scanf("%d", &n) != EOF ) { /* 读取裁判测试用例 */
for ( i=0; i<n; i++ ) {
for ( j=0; j<n; j++ )
scanf("%lf", &a[i][j]);
scanf("%lf", &b[i]);
}
/*--- 输出直接法的解 ---*/
if ( Direct(n, a, b) ) {
printf("Result of direct method:\n");
for ( j=0; j<n; j++ )
printf("%.8lf\n", b[j]);
}
else
printf("Doolittle factorization failed.\n");
printf("\n");
}
//system("pause");
fclose(stdin);
return 0;
}