2019年1月16日 星期三

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行相乘則得到:
L11 = a11 , L21 = a21, L31 = a31, L41 = a41
 
 L 矩陣的第一列和 U 矩陣的2,3,4行相乘則得到:
同理 U 的第二行乘以 L 的2,3,4列可得:

同理可得:
其一般化的式子為:
 
 j=1 時,   Li1 = ai1
 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]);
    }

}

沒有留言:

張貼留言

Messaging API作為替代方案

  LINE超好用功能要沒了!LINE Notify明年3月底終止服務,有什麼替代方案? LINE Notify將於2025年3月31日結束服務,官方建議改用Messaging API作為替代方案。 //CHANNEL_ACCESS_TOKEN = 'Messaging ...