C語言 例題6-1 (b) 高斯消去法 Gauss Elimination

C語言 例題6-1 (b) 高斯消去法 Gauss Elimination


-1.0 x1 + 1.0  x2  + 2.0  x3 = 2.0
 3.0 x1  - 1.0  x2  + 1.0  x3 = 6.0
-1.0 x1 +  3.0 x2  + 4.0  x3 = 4.0

/*
 Gauss Elimination Method in C.
 For this, let us first consider the following three equations:

 a1x + b1y + c1z = d1
 a2x + b2y + c2z = d2
 a3x + b3y + c3z = d3

 Assuming a1 ≠ 0, x is eliminated from the second equation
 by subtracting (a2/ a1) times the first equation
 from the second equation.

 In the same way, the C code presented here eliminates x
 from third equation by subtracting (a3/a1) times the first
 equation from the third equation.

 Then we get the new equations as:

 a1x + b1y + c1z = d1
 b’2y + c’2z = d’2
 c’’3z = d’’3

 The elimination procedure is continued until only one
 unknown remains in the last equation.

 After its value is determined, the procedure is stopped.
 Now, Gauss Elimination in C uses back substitution to get
 the values of x, y and z as:

 z=  d’’3 / c’’3
 y=(d’2 – c’2z) / b’2
 x=( d1- c1z- b1y)/ a1
*/

#include <stdio.h>
#include <math.h>
int main(void)
{
    int i,j,k,n;
    n=3;    //the order of matrix:
    float c,x[n+1],sum=0.0;

    //the elements of augmented matrix row-wise
    float A[3][4]= { {  2.0 ,  -1.5 , 3.0 , 1.0 } ,
                     { -1.0 ,     0 , 2.0 , 3.0 } ,
                     {  4.0 ,  -4.5 , 5.0 , 1.0 } };
   
    for(j=0; j<n; j++) /* loop for the generation of upper triangular matrix*/
    {
        for(i=0; i<n; i++)
        {
            if(i>j)
            {
                c=(A[i][j])/ (A[j][j]);
                printf("%0.2lf\t--> \t ", c);
                for(k=0; k<=n; k++)
                {
                    A[i][k]=(A[i][k]) -c* (A[j][k]);
                    printf("%0.2lf\t ", A[i][k]);
                }
                printf("\n");
            }
        }
    }
    x[n-1]=A[n-1][n]/A[n-1][n-1];
   
    printf("x[%d]=%0.2lf",n-1,x[n-1]);
   
    /* this loop is for backward substitution*/
    for(i=n-2; i>=0; i--)
    {
        sum=0;
        for(j=i; j<n; j++)
        {
            sum=sum+A[i][j]*x[j];
        }
        x[i]=(A[i][n]-sum)/A[i][i];
    }
    printf("\nThe solution is: \n");
    for(i=0; i<=n-1; i++)
    {
        printf("\nx%d=\t  %0.2lf\t",i,x[i]); /* x1, x2, x3 are the required solutions*/
    }
    return(0);
}

輸出畫面
-0.50  --> 0.00 -0.75 3.50 3.50
2.00   --> 0.00 -1.50 -1.00 -1.00
2.00   --> 0.00 0.00 -8.00 -8.00

x[2]=1.00
The solution is:

x0=   -1.00
x1=   -0.00
x2=   1.00