数值分析-杜利特尔分解法C语言代码

直接三角分解法

如果方程组的系数矩阵可以分解为一个下三角矩阵和一个上三角矩阵的乘积，即：A=LU。

Doolittle(杜利特尔)分解

在对系数矩阵进行三角分解时，如果L为单位下三角矩阵(主对角线系数全为1，对角线左下方的系数全不为0），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 y[Max];
double b[Max];
int num;

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]);
}

void Doolittle()
{
    int i,j,k;
    int t;
    double s,q;
    for(i=0;i<num;i++)//确定主对角线元素和0元素
    {
        for(j=0;j<num;j++)
    {
        if(i==j)
        {
            L[i][j]=1;
            U[i][j]=A[i][i];
        }
        if(i<j)
            L[i][j]=0;
        if(i>j)
            U[i][j]=0;
    }
    }
    for(k=0;k<num;k++)
    {
        for(j=k;j<num;j++)//计算矩阵u
        {
            s=0;
            for(t=0;t<k;t++)
            {
                s+=L[k][t]*U[t][j];
            }
            U[k][j]=A[k][j]-s;
        }
        for(i=k+1;i<num;i++)//计算矩阵L
        {
            q=0;
            for(t=0;t<k;t++)
            {
                q+=L[i][t]*U[t][k];
            }
            L[i][k]=(A[i][k]-q)/U[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");
    }
    //计算向量y
    y[0]=b[0];
    for(i=1;i<num;i++)
    {
        s=0;
        for(t=0;t<i;t++)
        {
            s+=L[i][t]*y[t];
        }
        y[i]=b[i]-s;
    }
    //计算向量x
    x[num-1]=y[num-1]/U[num-1][num-1];
    for(i=num-2;i>=0;i--)
    {
        q=0;
        for(t=i+1;t<num;t++)
        {
            q+=U[i][t]*x[t];
        }
        x[i]=(y[i]-q)/U[i][i];
    }
    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();
   Doolittle();
   return 0;
}