数值作业：Crout克洛特分解法之C语言代码

柔情只为你懂 2022-06-17 06:25 328阅读 0赞

矩阵crout分解的次序与Doolittle分解的次序不同，其他差不多．在矩阵A=LU的crout分解形式中，L是下三角矩阵，U是单位上三角矩阵  
![这里写图片描述][SouthEast]  
　　　当A的各阶主子式不为0,依次计算L的第一列，U的第一行，L的第K列，Ｕ的第K行，Doolitlle三角分解法跟这个恰好相反．依据Doolittle分解的手法，crout的算法可以比较容易实现，Crout分解题目： 不知道是分解形式Ax=b的形式还是，分解A=L\*U的形式，我这里是根据实验报告书，分解A=L\*U的形式，弄成Ax=b的形式修改下代码就行了．  
\[ 2 10 0 -3 \]  
\[-3 -4 -12 13\]  
\[ 1 2 3 -4 \]＝　Ｌ(单位下三角矩阵）　＊　　　Ｕ（上三角矩阵）  
\[ 4 14 9 -13\]  
下面上代码．

/********************************************
        > File Name: Crout.c
        > Author:chendiyang
        > School:WUST_CST_1501班
        > Myblog:www.chendsir.com
        > Mail:1441353519@qq.com 
        > Created Time: 2017年04月18日 星期二 22时25分45秒
     *****************************************************************
    #include <stdio.h>
    #include <math.h>
    #define MAX 20
    
    int main()
    {
        int n;
        int i,j,k;
        int mi;
        double mx,tmp;
        static double a[MAX][MAX],b[MAX],x[MAX],y[MAX];//方程组Ly=b和Ux=y
        static double l[MAX][MAX],u[MAX][MAX]; //l是矩阵L,u是矩阵U
        printf("\n 输入方程组的维数:");//输入AX=b的维数
        scanf("%d",&n);
        if(n>MAX)
        {
            return 1;
         } 
         if(n<=0)
         {
            return 1;
         }
         //输入矩阵的值
         printf("\n请输入A矩阵的值:"); 
         for(i=0;i<n;i++)
           for(j=0;j<n;j++) 
             scanf("%lf",&a[i][j]);
          //输入L矩阵的对角线的值
          for(i=0;i<n;i++)
          u[i][i]=1; //注意这里U是单位上三角矩阵，Doolittle中L是单位下三角矩阵
          for(k=0;k<n;k++)
          {
            //同时计算L的第1到K列元素
              for(i=k;i<n;i++)
              {
                 l[i][k]=a[i][k];
                 for(j=0;j<=k-1;j++)
                 l[i][k]-=(l[i][j]*u[j][k]);
              }
              //同时计算U的第1到K行元素
                for(j=k+1;j<n;j++)
              {
                  u[k][j]=a[k][j];
                  for(i=0;i<=k-1;i++)
                  u[k][j]-=(l[k][i]*u[i][j]);
                  u[k][j]/=l[k][k];
              }
          }
          for(i=0;i<n;i++)  //解方程组Ly=b
          {
              y[i]=b[i];
              for(j=0;j<=i-1;j++)
              y[i]-=(l[i][j]*y[j]);
              y[i]/=l[i][i];
          }
          for(i=n-1;i>=0;i--) //解方程组Ux=y
          {
              x[i]=y[i];
              for(j=i+1;j<n;j++)
              x[i]-=(u[i][j]*x[j]);
          }
          //输出运行结果 
              printf("solution is :\n");
              printf("矩阵L的值如下:\n");
              for(i=0;i<n;i++)
            {  
                printf("[");
                for(j=0;j<n;j++)
              {
                  printf("%lf ",l[i][j]);
              }
              printf("]\n");
            }
             printf("矩阵u的值如下:\n");
              for(i=0;i<n;i++)
            {  
                printf("[");
                for(j=0;j<n;j++)
              {
                  printf("%lf ",u[i][j]);
              }
              printf("]\n");
            }
              return 0;    
    
    }

运行结果：![这里写图片描述][SouthEast 1]

如果算法没有问题的话，L和U的结果应该是正确的，会算的可以手算验证一下，我是懒得算，如果你算的结果和我不同，可以在下方给我留言，方便及时改正．

[SouthEast]: /images/20220617/d4472c6ca0fe44fbade6e8e8333d5b70.png
[SouthEast 1]: /images/20220617/9c609f768d2f43a0a42b21f6f1a2c504.png