高斯消去法

高斯消去法

Ax=B

A=
  x1      x2      x3      x4
 4.01  1.23   1.43  -0.73   
 1.23  7.41   2.41   3.02 
 1.43  2.41   5.79  -1.11 
-0.73  3.02  -1.11   6.41   



B=
   5.94
   14.07
   8.52
   7.59

輸入資料

 n=4

 4.01 x1   1.23 x2   1.43 x3  -0.73 x4  =   5.94
 1.23 x1   7.41 x2   2.41 x3   3.02 x4  =  14.07
 1.43 x1   2.41 x2   5.79 x3  -1.11 x4  =   8.52
-0.73 x1   3.02 x2  -1.11 x3   6.41 x4   =  7.59


/* 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.

 4
 4.01  1.23   1.43  -0.73   5.94
 1.23  7.41   2.41   3.02  14.07
 1.43  2.41   5.79  -1.11   8.52
-0.73  3.02  -1.11   6.41   7.59
 */

 #include <stdio.h>
 #include <math.h>
 #define MAX 20
 void gaussh(int n,double a[MAX][MAX],double x[]);
 void main()
 {
    int i,j,k,m,n;
    double a[MAX][MAX],x[MAX];
    n=4;
    a[1][1]=4.01;   a[1][2]= 1.23;  a[1][3]=1.43;  a[1][4]=-0.73;  a[1][5]=5.94;
    a[2][1]=1.23;   a[2][2]= 7.41;  a[2][3]=2.41;  a[2][4]=3.02;   a[2][5]=14.07;
    a[3][1]=1.43;   a[3][2]= 2.41;  a[3][3]=5.79;  a[3][4]=-1.11;  a[3][5]=8.52;
    a[4][1]=-0.73;  a[4][2]= 3.02;  a[4][3]=-1.11;  a[4][4]=6.41;  a[4][5]=7.59;

    gaussh(n,a,x); /* call the function gaussh() */
    for(i=1;i<=n;i++)
       printf("x%d=%6.3lf\n",i,x[i]);
    return;
 }

 void gaussh(int n,double a[MAX][MAX],double x[])
 {
    int i,j,k,m;
    double temp,bb,cc;
    for(k=1;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=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];
       }
    }
   if(fabs(a[n][n])==0.0)
   {
      printf("NO UNIQUE SOLUTION!!!\n");
      exit(1);
   }
/* To start backward substitution */
   x[n]=a[n][n+1]/a[n][n];
   for(i=n-1;i>=1;i--)
   {
      cc=0.0;
      for(j=i+1;j<=n;j++)
         cc=cc+a[i][j]*x[j];
      x[i]=(a[i][n+1]-cc)/a[i][i];
   }
   return;
}

輸出畫面

x1= 1.000
x2= 1.000
x3= 1.000
x4= 1.000