#include <iostream> #include <cstdio> #include <algorithm> #include <cstdlib> using namespace std; #define N 20 double A[N][N],L[N][N],U[N][N],b[N],Y[N],X[N]; /// -------------------------------------------------------------------------文件处理 void saveLU(int n) { for(int i=0; i<n; i++) { for(int j=0; j<n; j++) { cout<<L[i][j]<<" "; } cout<<endl; } cout<<endl<<endl; for(int i=0; i<n; i++) { for(int j=0; j<n; j++) { cout<<U[i][j]<<" "; } cout<<endl; } } void saveT(double arr[], int n) { for(int i=0; i<n; i++) { cout<<arr[i]<<" "; } cout<<endl<<endl; } ///-------------------------------------------------------------------------初始化 void init(int n) { freopen("input.txt","r",stdin); freopen("lu.txt","w",stdout); for(int i=0; i<n; i++) { for(int j=0; j<n; j++) { cin>>A[i][j]; } } for(int i=0; i<n; i++) { cin>>b[i]; } } ///-------------------------------------------------------------------------直接法分解LU void breakdown(int n) { for (int i=0; i<n; i++) { U[0][i] = A[0][i]; ///U 第一行 L[i][0] = A[i][0]/U[0][0]; ///L 第一列 } ///U 第R行 L 第R列 double tmp = 0; for (int r=1; r<n; r++) { for (int i=r; i<n; i++) { tmp = A[r][i]; for(int k=0;k<=r-1; k++) { tmp -= L[r][k]*U[k][i]; } U[r][i] = tmp; tmp =A[i][r]; for(int k=0; k<=r-1; k++) { tmp -= L[i][k]*U[k][r]; } L[i][r] = tmp/U[r][r]; } } } ///-------------------------------------------------------------------------直接法计算Y void computeY(int n) { Y[0]=b[0]; ///自上往下 for (int i=1; i<n; i++) { Y[i] = b[i]; for (int j=0; j<=i-1; j++) { Y[i] -= L[i][j]*Y[j]; } } } ///-------------------------------------------------------------------------直接法计算X void computeX(int n) { int con = n; n--; X[n] = Y[n]/U[n][n]; ///自下往上 for (int i=n-1; i>=0; i--) { X[i] = Y[i]; for (int k=i+1; k<con; k++) { X[i] -= U[i][k]*X[k]; } X[i]/=U[i][i]; } } ///追赶法解三对角矩阵方程组{1.三对角矩阵LU分解 2.求y 3.求x } ///-------------------------------------------------------------------------1. LU分解 void TriangleBreakdown(int n) { for (int i=0; i<n; i++) { /// L的下对角线 U的主对角线可直接得出 U[i][i] = 1; if(i+1 < n) L[i+1][i] = A[i+1][i]; } L[0][0] = A[0][0]; U[0][1] = A[0][1]/L[0][0]; for (int i=1; i<n; i++) { ///L的下对角线 U的上对角线 L[i][i] = A[i][i] - L[i][i-1] * U[i-1][i]; if(i+1 < n) U[i][i+1] = A[i][i+1]/L[i][i]; } } ///------------------------------------------------------------------------- 计算X void TriangleY(int n) { Y[0] = b[0]/A[0][0]; for (int i=1; i<n; i++) { Y[i] = (b[i]-A[i][i-1]*Y[i-1])/L[i][i]; } } ///-------------------------------------------------------------------------计算Y void TriangleX(int n) { X[n-1] = Y[n-1]; for (int i=n-2; i>=0; i--) { X[i] = Y[i] - U[i][i+1] * X[i+1]; } } ///------------------------------------------------------三种方法整合 double AB[N][N+1]; void swap_cols(int x, int y, int n) ///交换两行 { double tmp = 0; for(int i=0; i<n+1; i++) { tmp = AB[x][i]; AB[x][i] = AB[y][i]; AB[y][i] = tmp; } } int find_max_col(int x, int n) /// 此列下方最大值 { double max1 = fabs(AB[x][x]); int res = x; for(int i=x+1; i<n; i++) { if(fabs(AB[i][x]) > max1) { max1 = AB[i][x]; res = i; } } return res; } void compute_gauss_X(int n) ///计算结果X { if(AB[n-1][n-1] == 0) cerr<<"wrong: divide 0 \n"; X[n-1] = AB[n-1][n]/AB[n-1][n-1]; double tmp =0; for(int i=n-2; i>=0; i--) { tmp = AB[i][n]; for(int j=i+1; j<n; j++){ tmp -= X[j]* AB[i][j]; //if(fabs(tmp)<10e-6) tmp = 0; } if(AB[i][i] != 0) X[i] = tmp/AB[i][i]; } } void solution_gauss(int n)/// 列主元高斯消元 { for(int i=0; i<n; i++) { for(int j=0; j<n; j++ ) { AB[i][j] = A[i][j]; } AB[i][n] = b[i]; } int pos = 0; double m = 0; for(int i=0; i<=n; i++) { ///标准行 pos =find_max_col(i, n); if( pos != i){ swap_cols(i, pos, n); } for(int i=0; i<n; i++) { for(int j=0; j<n+1; j++) { cout<<AB[i][j]<<" "; } cout<<endl; } cout<<"****************************\n"; for(int j=i+1; j<n; j++) { ///标准行以下 m = AB[j][i] / AB[i][i]; for(int k=i; k<n+1; k++) { ///此行所有数据 AB[j][k] -= m*AB[i][k]; } } cout<<"step# "<<i<<" :\n--------------------------\n"; for(int i=0; i<n; i++) { for(int j=0; j<n+1; j++) { cout<<AB[i][j]<<" "; } cout<<endl; } cout<<"----------------------------\n\n\n"; } compute_gauss_X(n); saveT(X,n); } void solution_LU(int n) ///直接法LU { breakdown(n); computeY(n); computeX(n); saveLU(n); saveT(Y,n); saveT(X,n); } void solution_triangle_chase(int n) ///追赶法 { TriangleBreakdown(n); TriangleY(n); TriangleX(n); saveT(Y,n); saveT(X,n); } int main() { int n = 3,choise=1; cout<<"选择方法: 1.Gauss 2.direct LU 3.triangle chase : \t\t"; cin>>choise; cout<<"输入矩阵大小\n"; cin>>n; init(n); switch(choise) { case 1: solution_gauss(n);break; case 2: solution_LU(n); break; case 3: solution_triangle_chase(n); } return 0; }