您的位置:首页 > 编程语言 > PHP开发

全选主元Gauss 消元法

2011-03-20 14:26 344 查看
这个算法是健壮的(数学术语)。我们通过选主元的方法来弥补顺序Gauss消元法的弊端,增加了pivrow和pivcol来记录包含主元素的行和列。由线性代数行列变换的知识知道,增广益矩阵行变换时,方程的解不变,但是列变换时,未知数的顺序发生了交换,因此需要记录列交换的信息在pivcol中,以便求解后还原。

全选主元高斯消元法的代码如下:

 

#include<iostream.h>
#include<fstream.h>
#include<math.h>
#include<process.h>

class gauss
{
private:
int i,j,k,n,*pivrow,**pivcol;
double eps,pivot,aik,al,sum,*x,**a;
public:
void gauss_input();
void gauss_elimination();
void gauss_output();

~gauss()
{
delete [] x;
delete [] pivrow;
for(i = 0;i < n; i++)
{
delete [] a[i];
}
delete [] a;

for(i=0;i<n;i++)
{
delete [] pivcol[i];
}
}
};

void main()
{
gauss solution;
solution.gauss_input();
solution.gauss_elimination();
solution.gauss_output();
}

void gauss::gauss_input()
{
ifstream fin("gauss.dat");
fin>>n;

x=new double
;
pivrow = new int
;
pivcol = new int *
;

for(i=0;i<n;i++)
{
pivcol[i] = new int[2];
}

a=new double*
;

for(i=0;i<n;i++)
{
a[i] = new double[n+1];
}

for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
fin>>a[i][j];
}

for(i=0;i<n;i++)
{
fin>>a[i]
;
}
fin>>eps;		//输入段结束
fin.close();
}

void gauss::gauss_elimination()	//gauss 消去函数
{
for(k=0;k<(n-1);k++)
{
pivot = fabs(a[k][k]);
pivrow[k] = k;
pivcol[k][0] = k;
pivcol[k][1] = k;

//下面是扫描所有行的最大值,并保存位置信息
for(i=k;i<n;i++)
for(j=k;j<n;j++)
{
if(pivot<fabs(a[i][j]))
{
pivot = fabs(a[i][j]);
pivrow[k] = i;
pivcol[k][1] = j;
}
}

if(pivot<eps)
{
cout<<"/n矩阵可能奇异.失败..."<<endl;
exit(0);

}

//交换行
if(pivrow[k]!=k)
{
for(j=k;j<(n+1);j++)
{
al = a[pivrow[k]][j];
a[pivrow[k]][j] = a[k][j];
a[k][j] = al;
}
}

if(pivcol[k][1]!=k)//交换列
{
for(i=0;i<n;i++)
{
al = a[i][pivcol[k][1]];
a[i][pivcol[k][1]] = a[i][k];
a[i][k] = al;
}
}
if(k!=(n-1))
{
for(i=(k+1);i<n;i++)
{
aik = a[i][k];
for(j=k;j<(n+1);j++)
{
a[i][j] -= aik*a[k][j]/a[k][k];
}
}
}
}
x[n-1] = a[n-1]
/a[n-1][n-1];//回代

for(i=(n-2);i>=0;i--)
{
sum = 0.0;
for(j=(i+1);j<n;j++)
{
sum += a[i][j] * x[j];
}
x[i] = (a[i]
- sum)/a[i][i];
}

for(k=(n-2);k>=0;k--)
{
al = x[pivcol[k][1]];
x[pivcol[k][1]]=x[pivcol[k][0]];
x[pivcol[k][0]] = al;
}
}

void gauss::gauss_output()	//结果输出函数
{
ofstream fout("gauss.out");
fout<<"/n结果是:"<<endl;
for(i=0;i<n;i++)
{
fout<<"/nx["<<i<<"]="<<x[i]<<endl;
}
fout.close();
for(i=0;i<n;i++)
{
cout<<"/nx["<<i<<"]="<<x[i]<<endl;
}
}


 
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息