数值分析-克劳特分解法C语言代码

Crout(克劳特)分解

如果在三角分解中，L为下三角矩阵，U为单位上三角矩阵，则成为克劳特分解。

定理：矩阵A有唯一克劳特分解的充分必要条件为A的前n-1个顺序主子式不为0.

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define Max 10

double A[Max][Max];
double L[Max][Max];
double U[Max][Max];
double x[Max];
double b[Max];
double y[Max];
int Num;

//输入系数矩阵A和右端项b
void Input_Matrix()
{
    int i,j;
    printf("系数矩阵A的阶数:\n");
    scanf("%d",&Num);
    for(i=1; i<=Num; i++)
    {
        printf("系数矩阵A的第%d行元素:\n",i);
        for(j=1; j<=Num; j++)
            scanf("%lf",&A[i-1][j-1]);
    }
    printf("右端项b:\n");
    for(i=1; i<=Num; i++)
    {
        scanf("%lf",&b[i-1]);
    }
    printf("输入的系数矩阵A:\n");
    for(i=0; i<Num; i++)
    {
        for(j=0; j<Num; j++)
            printf("%.4lf\t",A[i][j]);//\t表示加上之前的内容总共空8个字符
        printf("\n");
    }
    printf("输入的右端项b:\n");
    for(i=0; i<Num; i++)
        printf("%.4lf\n",b[i]);
}

//对系数矩阵进行Crout分解
void Crout()
{
    int i,j,k;
    int t;
    double sum_L,sum_U;
    for(i=0;i<Num;i++)
    {
        L[i][0]=A[i][0];
        U[i][i]=1;
    }
    for(j=1;j<Num;j++)
    {
        U[0][j]=A[0][j]/L[0][0];
    }
    for(k=0;k<Num;k++)
    {
        for(i=k;i<Num;i++)
        {
            sum_L=0;
            for(t=0;t<k;t++)
            {
                sum_L+=L[i][t]*U[t][k];
            }
            L[i][k]=A[i][k]-sum_L;
        }
        for(j=k+1;j<Num;j++)
        {
            sum_U=0;
            for(t=0;t<k;t++)
            {
                sum_U+=L[k][t]*U[t][j];
            }
            U[k][j]=(A[k][j]-sum_U)/L[k][k];
        }
    }
    printf("输出矩阵L\n");
    for(i=0;i<Num;i++)
    {
        for(j=0;j<Num;j++)
        {
            printf("%.4lf\t",L[i][j]);
        }
        printf("\n");
    }
    printf("输出矩阵U\n");
    for(i=0;i<Num;i++)
    {
        for(j=0;j<Num;j++)
        {
            printf("%.4lf\t",U[i][j]);
        }
        printf("\n");
    }
    double sum_y,sum_x;
    //计算y和x
    y[0]=b[0]/L[0][0];
    for(i=1;i<Num;i++)
    {
        sum_y=0;
        for(t=0;t<i;t++)
        {
            sum_y+=L[i][t]*y[t];
        }
        y[i]=(b[i]-sum_y)/L[i][i];
    }
    x[Num-1]=y[Num-1];
    for(i=Num-2;i>=0;i--)
    {
        sum_x=0;
        for(t=i+1;t<Num;t++)
        {
            sum_x+=U[i][t]*x[t];
        }
        x[i]=y[i]-sum_x;
    }
    printf("输出向量y\n");
    for(i=0;i<Num;i++)
    {
        printf("%.4lf\n",y[i]);
    }
    printf("输出向量x\n");
    for(i=0;i<Num;i++)
    {
        printf("%.4lf\n",x[i]);
    } 
}

int main()
{
    Input_Matrix();
    Crout();
    return 0;
}