【算法】高斯算法(1)

高斯算法(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");
}
 
阅读更多
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

关闭
关闭
关闭