C語言 例題6-2 利用 高斯-喬登 方法 Gauss-Jordan Method 求聯立方程式解

C語言 例題6-2 利用 高斯-喬登 方法 Gauss-Jordan Method  求聯立方程式解

/* ex6-3.c based on Gauss-Jordan Method to
 * solve n x n system of linear algebraic  equations.
 *
 * 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.
 *
      input matrix
    -x1 + x2 +  2x3 = 2
    3x1 - x2 +   x3 = 6
    -x1 + 3x2 + 4x3 = 4
 */
#include <stdio.h>
#include <math.h>

#define MAX 3
void gaussjd(int n,double a[MAX][MAX+1],double x[MAX]);
void printmatrix(int n , double a[MAX][MAX+1]) ;

int  main()
{
    int i,j,k,m,n;
    n=MAX;
    double x[MAX];
    double a[MAX][MAX+1]= { {-1.0 ,   1.0  , 2.0 ,  2.0},
                            { 3.0 ,  -1.0  , 1.0 ,  6.0},
                            {-1.0 ,   3.0  , 4.0 ,  4.0}
                        };
   
    printf("原始行列式\n");
    printf("================================\n");
    printmatrix(n,a);
   
    gaussjd(n,a,x);  /* call function gaussjd() */
   

   
    return 0;
 }

void gaussjd(int n,double a[MAX][MAX+1],double x[MAX])
{
    int i,j,k,m , cn=1;
    double temp,bb,cc,jj;
    for(k=0;k<n;k++)
    {
        /* check if a[k][k]=0 is true, then *
    * interchange E(k) and E(k+1)......*/
        if(a[k][k]==0)
        {
        for(m=1;m<=n+1;m++)
        {
            temp=a[k][m];
            a[k][m]=a[k+1][m];
            a[k+1][m]=temp;
        }
        }
    /* 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];
        }
        printf("\n");
        printf("===========%1d===================\n",cn);
        printmatrix(n,a);
        cn++;
    }
   
    printf("\n");
    printf("===========上三角矩陣================\n");
    printmatrix(n,a);
   
   
    if(fabs(a[n][n])==0.0)
    {
       printf("NO UNIQUE SOLUTION EXISTS!!!\n");
       return ;
    }
   
    for(k=n-1;k>=0;k--)
    {
        for(i=k-1;i>=0;i--)
        {
        jj=a[i][k]/a[k][k];
        for(j=n;j>=k;j--)
            a[i][j]=a[i][j]-jj*a[k][j];
        }
        printf("\n");
        printf("===========%1d===================\n",cn);
        printmatrix(n,a);
        cn++;
    }
   
    printf("\n");
    printf("===========對角矩陣================\n");
    printmatrix(n,a);
   
    printf("\n");
    /* starting to collect the solutions. */
   
    printf("\n行列式 解 \n");
    printf("================================\n");
    for(i=0;i<MAX;i++)
    {
        x[i]=a[i][n]/a[i][i];
        printf("x%d=%4.2lf\n",i,x[i]);
    }
    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");
    }     
}


輸出畫面
原始行列式
================================
-1.00 1.00 2.00 2.00
3.00 -1.00 1.00 6.00
-1.00 3.00 4.00 4.00

===========1===================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 2.00 2.00 2.00

===========2===================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 0.00 -5.00 -10.00

===========3===================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 0.00 -5.00 -10.00

===========上三角矩陣================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 0.00 -5.00 -10.00

===========4===================
-1.00 1.00 0.00 -2.00
0.00 2.00 0.00 -2.00
0.00 0.00 -5.00 -10.00

===========5===================
-1.00 0.00 0.00 -1.00
0.00 2.00 0.00 -2.00
0.00 0.00 -5.00 -10.00

===========6===================
-1.00 0.00 0.00 -1.00
0.00 2.00 0.00 -2.00
0.00 0.00 -5.00 -10.00

===========對角矩陣================
-1.00 0.00 0.00 -1.00
0.00 2.00 0.00 -2.00
0.00 0.00 -5.00 -10.00


行列式 解
================================
x0=1.00
x1=-1.00
x2=2.00