gaussian.c

数值分析 gauss jacobi消去法

/*  Gaussian elimitation */

#include "stdio.h"
#include "math.h"
#define N 11
#define M 10

void main()
{
    static float a[M][N];

 float x[M]={0};
    float r,s;
    int i,j,k;
    for (i=0;i<M;i++)
       {for(j=0;j<N;j++)
          {if (i==j) a[i][j]=3;
           else if (i==j+1) a[i][j]=9;
           else if (i==j-1&j!=N-1) a[i][j]=1;
           else if (j==N-1)
           {for (k=0;k<M;k++)
             {if (k!=0&k!=M-1)
              a[k][N-1]=13;
              else if (k==0)a[k][N-1]=4;
              else a[k][N-1]=12; }
              }
           else a[i][j]=0;
           printf("%f",a[i][j]);

           }
           printf("\n");
           }

   for (k=0;k<M-1;k++)
       {if (a[k][k]!=0)
        {   for (i=k+1;i<M;i++)
              r=a[i][k]/a[k][k];
              if (a[i][k]!=0)
                {for (j=k;j<N;j++)
                   a[i][j]=a[i][j]-r*a[k][j]; }
        else
         {for(i=0;i<N-1;i++)
           a[k][i]=a[k+1][i];

         }
         }

      else
           {for (i=k+1;i<M;i++)
              r=a[i][k]/a[k][k];
              if (a[i][k]!=0)
                {for (j=k;j<N;j++)
                   a[i][j]=a[i][j]-r*a[k][j];

                }

           }
        }
    for (i=0;i<N;i++)
        for(j=0;j<N;j++)
          printf("%f",a[i][j]);

       x[M-1]=a[M-1][N-1]/a[M-1][M-1];
          for(i=M-2;i>=0;i--)
            {
             s=a[i][N-1];
             for(j=i+1;j<M;j++)
             {
              s=s-a[i][j]*x[j];
              x[i]=s/a[i][i];
             }

            }
          for (i=0;i<M;i++)
                {printf("%f\n",x[i]);}

    getch();
}