计算方法 -- 解线性方程组直接法(LU分解、列主元高斯消元、追赶法)

#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 
";
    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<<"****************************
";
        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<<" :
--------------------------
";
        for(int i=0; i<n; i++) {
            for(int j=0; j<n+1; j++) {
                cout<<AB[i][j]<<" ";
            }
            cout<<endl;
        }
        cout<<"----------------------------


";
    }
    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 : 		";
    cin>>choise;
    cout<<"输入矩阵大小
";
    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;
}