LU分解法 EX6-5

* LU Decomposition  

源自於 

http://www.stat.nctu.edu.tw/misg/SUmmer_Course/C_language/Ch06/LUdecomposition.htm

將矩陣 A 分解成一上三腳矩陣 U 與一下三角矩陣 L所以方程式為：

令

則 x 由以下兩個方程式求解:

 

先算出b',再代回算出x.

　

(1)矩陣的分解

    4X4的矩陣為例 子:

　

將 L 矩陣每一列與與 U 矩陣 的1,2,3,4行相乘則得到:

L₁₁ = a₁₁ , L₂₁ = a₂₁, L₃₁ = a₃₁, L₄₁ = a₄₁

　

將 L 矩陣的第一列和 U 矩陣的2,3,4行相乘則得到:

同理 U 的第二行乘以 L 的2,3,4列可得:

同理可得:

其一般化的式子為：

　

當 j=1 時,   L_i1 = a_i1

當 i=1 時,

LU 分解法受歡迎的原因是不需要儲存 LU 中的 0 與 1,原矩陣A能表示LU兩個矩陣：

(2)演算法

(a) for i = 1 to n :do

L[i,1]=a[i,1]

(b) for j = 1 to n: do

U[1,j]=a[1,j]/L[1,1]

　

End i , j.

(c) for  j = 2 to n : do

           for  i = j to n : do

              sum=0.0

               for k = 1 to j-1 : do

                 sum = sum + L[i,k]*U[k,j]

               End k.

             L[i,j] = a[i,j] - sum

           End i.

         U[j,j] = 1.

          for i = j+1 to n : do

            sum=0.0

            for k = 1 to j-1 : do

                sum = sum + L[j,k]*U[k,i]

            End k.

          U[j,i] = (A[j,i] - sum)/L[j,j]

         End i.

    End j.

(d)反代

(3)程式

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#define    dim 10

double     a[dim][dim], L[dim][dim], U[dim][dim], x[dim];

double     b[dim], bp[dim];

mprint( int , int , double [][dim] , double [] );

main()

{

 int            n,i,j,k;

 double      temp;

　

  printf("input the dimension of matrix\n");

    scanf("%d",&n);

   if(n>=dim){

     printf("dimension over preset value %d\n",dim);

     exit(0);

     }

   printf("Input the a matrix\n\n");

  for(i=1; i<=n; i++){

     for(j=1;j<=n;j++){

       printf("input the element of a[%d,%d]\n",i,j);

       scanf("%lf",&a[i][j]);

      }

   }

  printf("\n Input the b vector\n\n");

   for(i=1; i<=n ; i++){

      printf("input the value of b[%d]\n",i);

      scanf("%lf",&b[i]);

    }

  mprint(n,n,a,b);

/* start to the procedure of LU decomposition */

   for(i=0; i<dim; i++){

     for(j=0; j<dim;j++){

       L[i][j]=0.0; U[i][j]=0.0;

       }

     }

  for(i=1; i<=n; i++){L[i][1]=a[i][1]; }

  for(j=1; j<=n;j++){U[1][j]=a[1][j]/L[1][1];}

　

 for(j=2; j<=n; j++){

      for(i=j; i<=n;i++){

         temp=0.0;

           for(k=1;k<=j-1;k++){

               temp=temp+L[i][k]*U[k][j];

            }

           L[i][j]=a[i][j]-temp;

        }

      U[j][j]=1.0;

     for(i=j+1;i<=n;i++){

         temp=0.0;

        for(k=1; k<=j-1;k++){

          temp=temp+L[j][k]*U[k][i];

         }

        U[j][i]=(a[j][i]-temp)/L[j][j];

       }

   }

/* Backward substitution */

    bp[1]=b[1]/L[1][1];

    for(i=2; i<=n; i++){

          temp=0.0;

        for(k=1; k<=i-1;k++){

          temp=temp+L[i][k]*bp[k];

         }

       bp[i]=(b[i]-temp)/L[i][i];

     }

  x[n]=bp[n];

  for(j=n-1; j>=1; j--)  {

        temp=0.0;

        for(k=j+1; k<=n; k++){

           temp=temp+U[j][k]*x[k];

         }

       x[j]=bp[j]-temp;

     }

mprint(n,n,l,b);

mprint(n,n,u,bp);

   for(j=1; j<=n; j++){

        printf(" x[%d]=%10.5lf\n",j,x[j]);

      }

}/*End of main() */

mprint(int m,int n,double aa[][dim],double bb[])

{

int i,j;

  for(i=1; i<=m; i++){

      for(j=1; j<=n;j++){

        printf("%10.5e ",aa[i][j]);

       }

    printf(" %10.5e \n",bb[i]);

    }

}