mifa.zip_反幂法_幂法

#include <stdio.h> #include <math.h> #define MAX_N 20 #define MAXREPT 100 #define epsilon 0.0001 void main() { int n,i,j,k; double err,xmax,oxmax; static double a[MAX_N][MAX_N]; static double l[MAX_N][MAX_N],u[MAX_N][MAX_N]; static double x[MAX_N],nx[MAX_N]; printf("INput n="); scanf("%d",&n); printf("Please input the matrix a[i][j].\n"); for(i=0;i<n;i++) for(j=0;j<n;j++) scanf("%lf",&a[i][j]); for(i=0;i<n;i++) { x[i]=1; printf("Y[0]=%f\t",x[i]); } oxmax=0; for (i=0;i<MAXREPT;i++) { for(j=0;j<n;j++) { nx[j]=0; for(k=0;k<n;k++) nx[j]+=a[j][k]*x[k]; printf("X[%-2d]=%f\t",i,nx[j]); } printf("\n"); xmax=0.0; for(j=0;j<n;j++) if(fabs(nx[j])>xmax) xmax=fabs(nx[j]); for(j=0;j<n;j++) { nx[j]/=xmax; printf("Y[%-2d]=%lf\t",i+1,nx[j]); } for(j=0;j<n;j++) x[j]=nx[j]; if (fabs(xmax-oxmax)<epsilon) { printf("\n按模最大：\n"); printf("solve...max lamda=%lf\n",xmax); printf("The vector is:\n"); for(i=0;i<n;i++) printf("%lf\n",nx[i]); break; } oxmax=xmax; } if(i==MAXREPT-1) printf("After %d repeat,max no result.\n",MAXREPT); for(i=0;i<n;i++) u[i][i]=1; for (k=0;k<n;k++) { for (i=k;i<n;i++) { l[i][k]=a[i][k]; for(j=0;j<=k-1;j++) l[i][k]-=l[i][j]*u[j][k]; } for (j=k+1;j<n;j++) { u[k][j]=a[k][j]; for(i=0;i<=k-1;i++) u[k][j]-=l[k][i]*u[i][j]; u[k][j]/=l[k][k]; } } for(i=0;i<n;i++) { x[i]=1; printf("Y[0]=%lf\t",x[i]); } oxmax=0; for (i=0;i<MAXREPT;i++) { for(j=0;j<n;j++) { nx[j]=x[j]; for(k=0;k<=j-1;k++) nx[j]-=l[j][k]*nx[k]; nx[j]/=l[j][j]; } for (j=n-1;j>=0;j--) { x[j]=nx[j]; for(k=j+1;k<n;k++) x[j]-=u[j][k]*x[k]; } for(j=0;j<n;j++) printf("X[%-2d]=%lf\t",i,x[j]); printf("\n"); xmax=0.0; for(j=0;j<n;j++) if(fabs(x[j])>xmax) xmax=fabs(x[j]); for(j=0;j<n;j++) { x[j]/=xmax; printf("Y[%-2d]=%lf\t",i+1,x[j]); } if (fabs(xmax-oxmax)<epsilon) { printf("\n按模最小：\n"); printf("solve...max lamda=%lf\n",1/xmax); printf("The vector is:\n"); for(i=0;i<n;i++) printf("%lf\n",x[i]); break; } oxmax=xmax; } if(i==MAXREPT-1) printf("After %d repeat,max no result.\n",MAXREPT); }