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

Solve the following system of equations using LU Decomposition method:

  

Solution: Here, we have

A =  and  such that A X = C.

Mathematics | L U Decomposition of a System of Linear Equations

源自於 https://www.geeksforgeeks.org/l-u-decomposition-system-linear-equations/

Example:
Solve the following system of equations using LU Decomposition method:

  

Solution: Here, we have

A =  and  such that A X = C.

Now, we first consider  and convert it to row echelon form using Gauss Elimination Method.

So, by doing

(1) 

(2) 

we get

 

Now, by doing

(3) 

we get

(Remember to always keep ‘ – ‘ sign in between, replace ‘ + ‘ sign by two ‘ – ‘ signs)

Hence, we get L =  and U = 

(notice that in L matrix,  is from (1),  is from (2) and  is from (3))

Now, we assume Z  and solve L Z = C.

 

So, we have   

Solving, we get ,  and  .

Now, we solve U X = Z

 

Therefore, we get   , 

Thus, the solution to the given system of linear equations is , ,  and hence the matrix X = 

/************** 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[3][3]= { {  1,  1,  1 },

                     {  4,  3, -1 },

                     {  3,  5,  3 } };

                    

  

    float b[3]= {1,6,4};

  

    for(i=0;i<n;i++)

    {

        printf("\nRow%1d\t",i);

        for(j=0;j<n;j++)

            printf("%0.2lf\t\t",a[i][j]);

    }

  

    printf("\n\nelements of b matrix\n");

    for(i=0;i<n;i++)

        printf("%0.2lf\t\t",b[i]);

      

    //********** LU decomposition *****//

    for(k=0;k<n;k++)

    {

        u[k][k]=1;

        for(i=k;i<n;i++)

        {

            sum=0;

            for(p=0;p<k;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=0;p<k;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");

    //==================L matrix===============

    for(i=0;i<n;i++)

    {

        for(j=0;j<n;j++)

            printf("%0.2lf\t",l[i][j]);

        printf("\n");

    }

    printf("\n\n");

    //==================U matrix===============

    for(i=0;i<n;i++)

    {

        for(j=0;j<n;j++)

            printf("%0.2lf\t",u[i][j]);

        printf("\n");

    }

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

    for(i=0;i<n;i++)

    {                                        //forward subtitution method

        sum=0;

        for(p=0;p<i;p++)

        sum+=l[i][p]*z[p];

        z[i]=(b[i]-sum)/l[i][i];

    }

    //********** FINDING X; UX=Z***********//

    for(i=n-1;i>=0;i--)

    {

        sum=0;

        for(p=n-1;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=0;i<n;i++)

        printf("%0.2lf\t",x[i]);

    return 0;

}

輸出畫面

The order of square matrix:  3

Row0 1.00 1.00 1.00

Row1 4.00 3.00 -1.00

Row2 3.00 5.00 3.00

elements of b matrix

1.00 6.00 4.00

LU matrix is : 

1.00 0.00 0.00

4.00 -1.00 0.00

3.00 2.00 -10.00

1.00 1.00 1.00

0.00 1.00 5.00

0.00 0.00 1.00

Set of solution is

1.00 0.50 -0.50