C語言 例題6-1 (b) 利用Guassian Elimination 解行列式

C語言 例題6-1 (b)  利用Guassian Elimination 解行列式

 3 --> MAX 
 2.0 x0  -1.5 x1  + 3.0  x2 =  1.0
-1.0 x0   0.0 x1  + 2.0  x2 =  3.0
 4.0 x0 -4.5  x1  + 5.0  x2 =  1.0


/* ex6-1.c based on Gaussian Elimination method
 * for solving the n x n linear algebra system
 * a11 x1+a12 x2+...+a1n xn=b1
 * a21 x1+a22 x2+...+a2n xn=b2
 * .       .         .       .
 * .       .         .       .
 * an1 x1+an2 x2+...+ann xn=bn
 * Input number of unknowns and equations n
 * with coefficent a11,a12,...,ann and b1,b2,
 * ...bn. Output solution x1,x2,x3,...,xn.
 *

 3 --> MAX
 2.0  -1.5  3.0  1.0
-1.0   0.0  2.0  3.0
 4.0  -4.5  5.0  1.0
*/
#include <stdio.h>
#include <math.h>
#define MAX 3
void gaussh(int n,double a[MAX][MAX+1],double x[]);
void printmatrix(int n , double a[MAX][MAX+1]) ;
void backSub(double a[MAX][MAX+1]) ;

int main()
{   
    int i,j,k,m,n;
    n=MAX;
    double x[MAX];
    double a[MAX][MAX+1]= { {2.0 , -1.5 , 3.0 , 1.0},
                          {-1.0 ,  0.0 , 2.0 , 3.0},
                          {4.0 , -4.5 , 5.0 , 1.0}
                        };

    printf("原始行列式\n");
    printf("================================\n");
    printmatrix(n,a);
   

    gaussh(n,a,x); /* call the function gaussh() */
   

    return 0;
}


void gaussh(int n,double a[MAX][MAX+1],double x[])
{
    int i,j,k,m;
    double temp,bb,cc;
    for(k=0;k<=n-1;k++)
    {
       /* check if a[k][k]=0 is true then interchange */
       /* E(k) and E(k+1).............................*/
       if(a[k][k]==0)
       {
        for(m=0;m<=n;m++)
        {
            temp=a[k][m];
            a[k][m]=a[k+1][m];
            a[k+1][m]=temp;
        }
       }
     
       printf("\n");
       printf("================================\n");
       printmatrix(n,a);
     
       /* To reduce the matrix to triangular form */
       for(i=k;i<=n-1;i++)
       {
        bb=a[i+1][k]/a[k][k];
        for(j=k;j<=n+1;j++)
            a[i+1][j]=a[i+1][j]-bb*a[k][j];
       }
    }
   
   
    if(fabs(a[n-1][n-1])==0.0)
    {
        printf("NO UNIQUE SOLUTION!!!\n");
        return;
    }
   

   
    printf("\n\n上三角矩陣\n");
    printf("================================\n");
    printmatrix(n,a);
   
    backSub(a);
    return;
}

void printmatrix(int n , double a[MAX][MAX+1])
{
    for (int i=0; i<n; i++)
    {
        for (int j=0; j<=n; j++)
            printf("%4.2lf ", a[i][j]);
        printf("\n");
    }       
}


// function to calculate the values of the unknowns
void backSub(double a[MAX][MAX+1])
{
    double x[MAX];  // An array to store solution
 
    /* Start calculating from last equation up to the
       first */
    for (int i = MAX-1; i >= 0; i--)
    {
        /* start with the RHS of the equation */
        x[i] = a[i][MAX];
 
        /* Initialize j to i+1 since matrix is upper
           triangular*/
        for (int j=i+1; j<MAX; j++)
        {
            /* subtract all the lhs values
             * except the coefficient of the variable
             * whose value is being calculated */
            x[i] -= a[i][j]*x[j];
        }
 
        /* divide the RHS by the coefficient of the
           unknown being calculated */
        x[i] = x[i]/a[i][i];
    }
   
    printf("\n行列式 解 \n");
    printf("================================\n");
    for (int i=0; i<MAX; i++)
        printf("x%1d = %3.3lf\n", i, x[i]);
}

輸出畫面
原始行列式
================================
2.00 -1.50 3.00 1.00
-1.00 0.00 2.00 3.00
4.00 -4.50 5.00 1.00

================================
2.00 -1.50 3.00 1.00
-1.00 0.00 2.00 3.00
4.00 -4.50 5.00 1.00

================================
2.00 -1.50 3.00 1.00
0.00 -0.75 3.50 3.50
0.00 -1.50 -1.00 -1.00

================================
2.00 -1.50 3.00 1.00
0.00 -0.75 3.50 3.50
0.00 0.00 -8.00 -8.00


上三角矩陣
================================
2.00 -1.50 3.00 1.00
0.00 -0.75 3.50 3.50
0.00 0.00 -8.00 -8.00

行列式 解
================================
x0 = -1.000
x1 = -0.000
x2 = 1.000