高斯算法(1)

#include "stdio.h"

#include "math.h"

#include "stdlib.h"

void max_ele(double af[][21],int n,int k)  /* 選主元 */

  {double max1,t;

   int i,j;

   max1=af[k][k]; i=k;

   for (j=k+1;j<=n;j++)

     if(fabs(af[j][k])>fabs(max1)) i=j;

   if(i>k)

     for (j=k;j<=n+1;j++)

       {t=af[i][j],af[i][j]=af[k][j],af[k][j]=t;}

  }

int gauss(double a[][20],double f[],double x[],int n)

/* Guass 消去法,引入af[20[21],x[]是為了不改變矩陣a[][],f[]

   及簡化程序  */

 {int i,j,k;

  double af[20][21],del,t;

  for(i=1;i<=n;i++)

    { for(j=1;j<=n;j++)  af[i][j]=a[i][j];

      af[i][n+1]=f[i];

    }

  for (k=1;k<n;k++)

    { max_ele(af,n,k);

      del=af[k][k];

      if (fabs(del)<1e-7)

       {printf("Single! press any key return...\n");

        getchar();

        return 0;}

      for (j=k;j<=n+1;j++)  af[k][j]/=del;

      for (i=k+1;i<=n;i++)

       { t=af[i][k];

         for (j=k;j<=n+1;j++)

           af[i][j]-=af[k][j]*t;

       }

     }

  del=af[n][n];

  for (j=n;j<=n+1;j++) af[n][j]/=del;

  for (j=1;j<=n;j++) x[j]=af[j][n+1];

  for (i=n-1;i>0;i--)

    for (j=n;j>i;j--)  x[i]-=x[j]*af[i][j];

  return 1;

 }

void test1()

 /* 檢驗 Guass消去法 */

 { double a[20][20], f[20],x[20];

   int i,j,n;

   printf("Input n(<20)");

   scanf("%d",&n);

   /* 隨機生成系數矩陣及右端項,進行檢驗 */

   for (i=1;i<=n;i++)

      for (j=1;j<=n;j++) a[i][j]=(float)rand()/100.0;

   for (i=1;i<=n;i++)

     { for (j=1;j<=n;j++) printf("%8.4f ",a[i][j]);

       printf("\n");}

   for (i=1;i<=n;i++)

     { f[i]=0;

       for (j=1;j<=n;j++) f[i]+=a[i][j];

      }

   gauss(a,f,x,n);

   printf("the system solution is:\n");

   for (i=1;i<=n;i++)  printf("%lf ",x[i]);

   printf("\n");

 }

int inver(double a[20][20],double e[20][20],int n)

 /* 用追趕法求三對角陣 a[][] 的逆 */

 {double ae[20][40],gamma,del;

  int i,j;

  gamma=-a[2][1];

  for (i=1;i<=n;i++)

    {for (j=1;j<=n;j++)

      { ae[i][j]=a[i][j]; ae[i][n+j]=0; }

      ae[i][n+i]=1;

     }

  for(i=1;i<=n-1;i++)

    {if((del=ae[i][i])==0) return 0;

     for (j=i;j<=n*2;j++)

       { ae[i][j]/=del; ae[i+1][j]+=ae[i][j]*gamma; }

     }

  del=ae[n][n];

  for(j=n;j<=n*2;j++)

    ae[n][j]/=del;

  for(i=1;i<=n;i++)

    {for (j=1;j<=2*n;j++)

       printf("%8.4f ",ae[i][j]);

    printf("\n");

    }

  for(i=n;i>1;i--)

    { del=ae[i-1][i];

      for (j=i;j<=n*2;j++)

       ae[i-1][j]-=ae[i][j]*del;

     }

  for (i=1;i<=n;i++)

    for (j=1;j<=n;j++)

      e[i][j]=ae[i][n+j];

  return 1;

  }

void mult(double a[][20],double e[][20],double f[][20],int n)

 {double s;

  int i,j,k;

  for (i=1;i<=n;i++)

    for (j=1;j<=n;j++)

      { s=0;

       for (k=1;k<=n;k++)

        s+=a[i][k]*e[k][j];

       if (fabs(s)<1e-7) s=fabs(s);

       f[i][j]=s;

      }

 }

void minus(double a[][20],double b[][20],double c[][20],int n)

 {int i,j;

  for (i=1;i<=n;i++)

    for (j=1;j<=n;j++)

      c[i][j]=a[i][j]-b[i][j];

  }

main()

{

 int i,j,n,m;

 double gamma,beta,beta1,a[20][20],b[20][20],c[20][20],d[20][20],

       e[20][20],f[20][20],deltag,deltax;

 double gx[20]={0,1,0,1,0,-1}, r[20][20],x[20];

 printf("請輸入主梁數目,gamma,beta,beta1:\n");

 scanf("%d,%lf,%lf,%lf",&n,&gamma,&beta,&beta1);

 m=n-1;

 for (i=1;i<=m;i++)

   for (j=1;j<=m;j++)

     a[i][j]=b[i][j]=c[i][j]=d[i][j]=0;

 deltag=2*(1+gamma+beta);

 deltax=2*(gamma+3*beta1);

 printf("%8.4f  %8.4f\n",deltag,deltax);

 for (i=1;i<=m;i++)

  { a[i][i]=deltax;

    b[i][i]=deltag;}

 for (i=1;i<=m-1;i++)

  { a[i][i+1]=a[i+1][i]=-gamma;

    b[i][i+1]=b[i+1][i]=gamma-1;

    c[i][i+1]=d[i+1][i]=gamma;

    c[i+1][i]=d[i][i+1]=-gamma;}

 inver(a,e,m);

 for (i=1;i<=m;i++)

   { for (j=1;j<=m;j++)

     printf("%6.3f ", e[i][j]);

     printf("\n"); }

 mult(c,e,f,m);

 mult(f,d,e,m);

 minus(b,e,r,m);

 gauss(r,gx,x,m);

 printf("The result is:\n");

 for (j=1;j<=m;j++)

   printf("%8.4f ", x[j]);

 printf("\n");

}