C語言 例題6-8 LU分解法求線性代數解

C語言 例題6-5 LU分解法求線性代數解

-1.00x1    1.00x2  2.00x3  = 2.00
 3.00x1   -1.00x2  1.00x3  = 6.00
-1.00x1    3.00x2  4.00x3  = 4.00


/************** LU Decomposition for solving linear equations ***********/
/*
[A] = [L][U]
where, L is the lower triangular matrix and
U is the upper triangular matrix.
The elements of these three matrices row-wise are:

[A] = { a11, a12, a13,
        a21, a22, a23,
        a31, a32, a33 }

[L] = {  1,   0,  0,
       l21,   1,  0,
       l31, l32,  1}
     
[U] = {u11, u12, u13,
       0  , u22, u23,
       0  ,   0, u33}


-1.00x1    1.00x2  2.00x3  = 2.00
 3.00x1   -1.00x2  1.00x3  = 6.00
-1.00x1    3.00x2  4.00x3  = 4.00

x3 = 2.000   x2 = -1.000   x1 = 1.000
*/

#include<stdio.h>
#include<math.h>
int main()
{
    int n,i,k,j,p;
    float l[10][10]={0},u[10][10]={0},sum,z[10]={0},x[10]={0};
    n=3;
 
    printf("The order of square matrix: %2d\n",n);
    float a[4][4]= { {  0,  0 , 0, 0 },
                     {  0, -1,  1, 2 },
                     {  0,  3, -1, 1 },
                     {  0, -1,  3, 4 } };
                   
 
    float b[4]= {0,2,6,4};
 
    for(i=1;i<=n;i++)
    {
        printf("\nRow%1d\t",i);
        for(j=1;j<=n;j++)
            printf("%0.2lf\t\t",a[i][j]);

    }
 
    printf("\n\nelements of b matrix\n");
    for(i=1;i<=n;i++)
        printf("%0.2lf\t\t",b[i]);
     
    //********** LU decomposition *****//
    for(k=1;k<=n;k++)
    {
        u[k][k]=1;
        for(i=k;i<=n;i++)
        {
            sum=0;
            for(p=1;p<=k-1;p++)
                sum+=l[i][p]*u[p][k];
            l[i][k]=a[i][k]-sum;
        }

        for(j=k+1;j<=n;j++)
        {
            sum=0;
            for(p=1;p<=k-1;p++)
                sum+=l[k][p]*u[p][j];
            u[k][j]=(a[k][j]-sum)/l[k][k];
        }
    }
    //******** Displaying LU matrix**********//
    printf("\n\nLU matrix is : \n");
    for(i=1;i<=n;i++)
    {
        for(j=1;j<=n;j++)
            printf("%0.2lf\t",l[i][j]);
        printf("\n");
    }
    printf("\n\n");
 
    for(i=1;i<=n;i++)
    {
        for(j=1;j<=n;j++)
            printf("%0.2lf\t",u[i][j]);
        printf("\n");
    }

    //***** FINDING Z; LZ=b*********//

    for(i=1;i<=n;i++)
    {                                        //forward subtitution method
        sum=0;
        for(p=1;p<i;p++)
        sum+=l[i][p]*z[p];
        z[i]=(b[i]-sum)/l[i][i];
    }
    //********** FINDING X; UX=Z***********//

    for(i=n;i>0;i--)
    {
        sum=0;
        for(p=n;p>i;p--)
            sum+=u[i][p]*x[p];
        x[i]=(z[i]-sum)/u[i][i];
    }
    //*********** DISPLAYING SOLUTION**************//
    printf("\n\nSet of solution is\n");
    for(i=1;i<=n;i++)
        printf("%0.2lf\t",x[i]);

    return 0;
}

輸出畫面
The order of square matrix:  3

Row1 -1.00 1.00 2.00
Row2 3.00 -1.00 1.00
Row3 -1.00 3.00 4.00

elements of b matrix
2.00 6.00 4.00

LU matrix is :
-1.00 0.00 0.00
3.00 2.00 0.00
-1.00 2.00 -5.00


1.00 -1.00 -2.00
0.00 1.00 3.50
0.00 0.00 1.00


Set of solution is
1.00      -1.00   2.00