gaussion.c

数值分析 gauss jacobi消去法

       #include "stdio.h"
       #include "math.h"
       #define  n  10
       main()
      {
       float a[n+1][n+1];
       float b[n],x[n],h[n+1];
       int i,j,k,c;
       float s=0,MAX;
       printf("This is the Gaussion elimination method\n\n");
       for(i=0;i<n;i++)                 /*初始化矩阵a[n][n]*/
         for(j=0;j<n;j++)
         {
          if (i==j)
                a[i][j]=3;
           else if (i==j-1)
                a[i][j]=1;
           else if (i==j+1)
                a[i][j]=9;
           else  a[i][j]=0;
          }
       b[0]=4;                       /*初始化b[n]*/
       b[n-1]=12;
       for(i=1;i<n-1;i++)
           b[i]=13;
       for(i=0;i<n;i++)
       a[i][n]=b[i];
        /*for(i=0;i<n;i++)
         for(j=0;j<=n;j++)
       printf("%f ",a[i][j]);
       printf("\n\n");*/
       for(i=0;i<n-1;i++)           /*高斯行主元素消去法*/
       {
        MAX=fabs(a[i][0]);          /*比较*/
        c=i;
        for(j=0;j<=n;j++)
          if(fabs(a[i][j])>MAX)
            c=j;
        if(c!=i)                    /*交换*/
          for(j=0;j<n;j++)
            { h[j]=a[j][c];
              a[j][c]=a[j][i];
              a[j][i]=h[j];
             }
        for(j=i+1;j<=n;j++)         /*消元*/
          h[j]=a[j][i]/a[i][i];
       for(j=i+1;j<=n;j++)
          for(k=i+1;k<=n;k++)
            a[j][k]=a[k][k]-h[j]*a[j-1][k];
        if(a[i][i]==0)
            {for(j=0;j<n;j++)
              {h[j]=a[i][j];
               a[i][j]=a[i+1][j];
               a[i+1][j]=h[j];
              }
             h[1]=b[i];
             b[i]=b[i+1];
             b[i+1]=h[1];
            i=i+2;
            }
        }
        for(i=0;i<n;i++)          /*输出a的增广矩阵*/
         for(j=0;j<=n;j++)
       printf("%f ",a[i][j]);
       printf("\n\n");
       for(i=0;i<n;i++)
         b[i]=a[i][n];
       for(i=n-1;i>=0;i--)      /*求解方程组的解*/
       {
         for(j=i+1;j<=n-1;j++)
           s=s+a[i][j]*x[j];
         x[i]=(b[i]-s)/a[i][i];
         s=0;
       }
       for(i=0;i<n;i++)           /*输出结果*/
         printf("%f   ",x[i]);
         getch();
       }