gaosi.rar_高斯算法

#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"); }