#include<iostream>
#include<math.h>
using namespace std;
#define MAX 100
double A[MAX+1][MAX+1];
double B[MAX+1];
double X[MAX+1];
double Z[MAX+1];
int D[MAX+1]; //未知變量位置的變化
int n;
int e;
void input()
{
int i,j;
printf("n:");
scanf("%d",&n);
printf("A[][]:\n");
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
scanf("%lf",&A[i][j]);
printf("B[]:\n");
for(i=1;i<=n;i++)
scanf("%lf",&B[i]);
printf("e:");
scanf("%lf",&e);
}
void SwapE(double a,double b) //swap elements
{
double T;
T=a;
a=b;
b=T;
}
void SwapR(int k,int kmi)
{
int j;
for(j=k;j<=n;j++)
SwapE(A[k][j],A[kmi][j]);
}
void SwapC(int k,int kmj)
{
int i;
for(i=k;i<=n;i++)
SwapE(A[i][k],A[i][kmj]);
}
void AllGaussianElimination()
{
int kmi,kmj;//the i and j of the max(abs) element when k
int i,j,k;
double T;
for(i=1;i<=n;i++)//初始化位置變量位置
D[i]=i;
for(k=1;k<=n-1;k++)
{
//選主元
T=0;
for(i=k;i<=n;i++)
for(j=k;j<=n;j++)
if(fabs(A[i][j])>T) { T=fabs(A[i][j]); kmi=i;kmj=j;}
if(T<=e) {printf("Error!\n"); return ;}
if(kmi!=k) { SwapR(k,kmi); SwapE(B[k],B[kmi]); }
if(kmj!=k) { SwapC(k,kmj); SwapE(D[k],D[kmj]); }
//消元
for(i=k+1;i<=n;i++)
{
T=A[i][k]/A[k][k];
B[i]-=T*B[k];
for(j=k;j<=n;j++)
A[i][j]-=T*A[k][j];
}
//回代
if(A[n][n]<=e) {printf("Error!\n");return ;}
Z[n]=B[n]/A[n][n];
double S_Aij_Zj;
for(i=n-1;i>=1;i--)
{
S_Aij_Zj=0;
for(j=i+1;j<=n;j++)
S_Aij_Zj+=A[i][j]*Z[j];
Z[i]=(B[i]-S_Aij_Zj)/A[i][i];
}
for(j=1;j<=n;j++)
X[D[j]]=Z[j];
}
}
void print(double X[])
{
int i;
printf("X[]:\n");
for(i=1;i<=n;i++)
printf("%f\n",X[i]);
}
int main()
{
input();
AllGaussianElimination();
print(X);
system("pause");
return 0;
}posted on 2009-05-16 16:37
wyiu 閱讀(802)
評論(0) 編輯 收藏 引用 所屬分類:
數值分析