运用Doolitle分解法解线性方程组

Doolittle分解法是将系数矩阵A分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，即
A=L*U，其中L和U的形式为
L=，U=

然后通过公式L*Y=b（顺代）解得Y；
最后通过公式Y=UX（回代）解得X。
运用Dolittle分解法求解线性方程组的基本步骤为：
（1）输入方程组的阶数n，系数矩阵A和右端的常系数矩阵b；
（2）L矩阵对角元素赋值为1；
（3）按顺序（k=1,2,…,n）先行后列交替计算U和L的元素：
①对j = 1,2,…,n，计算，
对i = 2,3,…,n，计算
②计算U的第二行， L的第二列，U的第三行，L的第三列…
计算U的第k行和L的第k列元素公式为：
，j = k,k+1,…,n
，i = k+1,k+2…,n
（4）计算Y和X：
顺代计算
，k =2 ,3,…,n
回代计算
，k=n-1,…,2,1
（5）输出方程组的解，结束。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define INF 0x3f3f3f3f
#pragma GCC optiize(2)
#define maxn 0x3f3f3f3f
double A[11][11],b[11];
double l[11][11],u[11][11];
double c[11][11];
double x[11],y[11];
int main()
{
    int n;
    cin>>n;
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            cin>>A[i][j];
        }
    }
    for(int i=1;i<=n;i++)
    {
        cin>>b[i];
    }
    for(int i=1;i<=n;i++)
    {
        l[i][i]=1;
    }
    for(int i=1;i<=n;i++)
    {
        for(int j=i;j<=n;j++)
        {
            if(i==1)
            {
                u[i][j]=A[i][j];
            }
            else
            {
                int sum=0;
                for(int t=1;t<=i-1;t++)
                {
                    sum+=u[t][j]*l[i][t];
                }
                u[i][j]=A[i][j]-sum;
            }
        }
        for(int j=i+1;j<=n;j++)
        {
            if(i==1)
            {
                l[j][i]=A[j][i]/u[1][1];
            }
            else
            {
                int sum=0;
                for(int t=1;t<=i-1;t++)
                {
                    sum+=l[j][t]*u[t][i];
                }
                l[j][i]=(A[j][i]-sum)/u[i][i];
            }
        }
    }
    cout<<"U :"<<endl;
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            cout<<u[i][j]<<' ';
        }
        cout<<endl;
    }
    cout<<"L :"<<endl;
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            cout<<l[i][j]<<' ';
        }
        cout<<endl;
    }

    for(int i=1;i<=n;i++)
    {
        if(i==1)
        {
            y[i]=b[i];
        }
        else
        {
            double sum=0;
            for(int t=1;t<=i-1;t++)
            {
                sum+=l[i][t]*y[t];
            }
            y[i]=b[i]-sum;
        }
    }
    for(int i=n;i>=1;i--)
    {
        if(i==n)
        {
            x[i]=y[n]/u[n][n];
        }
        else
        {
            double sum=0;
            for(int t=i+1;t<=n;t++)
            {
                sum+=u[i][t]*x[t];
            }
            x[i]=(y[i]-sum)/u[i][i];
        }
    }
    cout<<"Xi :"<<endl;
    for(int i=1;i<=n;i++)
        cout<<x[i]<<' ';
    cout<<endl;
    cout<<"Yi :"<<endl;
    for(int i=1;i<=n;i++)
        cout<<y[i]<<' ';
    cout<<endl;
    return 0;
}