高斯·约当消元法

高斯·约当消元法

By CaptainChen 高斯消元法，用于解多元一次方程（几乎类似模拟手动解方程）。

思路：

通过等式的乘除，把方程1的x1x1系数a11a11分别化为方程2~方程n的x1x1系数，然后将方程2~方程n减去得到的新方程，从而消掉方程2~方程n中的x1x1。接着用方程2的x2x2继续把方程3~方程n中的x2x2消掉……

大概系数就成了这个样子↓

举个例子：

一个问题：遇到这种情况

0.0000000000001×1x1++x2x2==120.0000000000001×1+x2=1x1+x2=2

就会出现精度问题，处理成

x10++1000000000000x2999999999999x2==10000000000002x1+10000 00000000x2=10000000000000+999999999999x2=2

而当位数一多，我们的计算机就有可能存不了那么多9，从而忽略掉一堆，原本的误差就扩散得越来越大，这时候我们需要选主消元（比如把刚才两个式子交换一下），选择当前这一项系数绝对值最大的方程作为主元，来消掉其它方程。

高斯·约当消元：

高斯消元是当执行到第i个方程的x k xk时，消掉i以后的x k xk。而约当就同时消掉i以前的x k xk，使系数变成一条对角线

解得情况：高斯消元完成后，下面若有方程系数全部为0

0×1+0×2+0×3+...+0×n=00×1+0×2+0×3+...+0×n=0

则说明有多解

如果为

0×1+0×2+0×3+...+0×n=非00×1+0×2+0×3+...+0×n=非0 则无解

代码：

#include

#include

#include

using std::swap;

#define MAXN 405

#define EPS 1e-8

typedef double Matrix[MAXN][MAXN];//系数矩阵

int n,m;

int Rank;//有效方程的行数，Rank之后的方程x系数为0

double X[MAXN];//解

bool Free[MAXN];//是否为自由变量

Matrix A;//系数矩阵

//浮点数比较

int fcmp(double a,double b)

{

if((a-b)>EPS)

return1;

else if((a-b)>-EPS)

return0;

return -1;

}

//高斯·约当消元

void Gauss()

{

int r,c,mxr;

for(r=1,c=1;r<=n&&c

{

//寻找绝对值最大（选主）

mxr=r;

for(int i=r+1;i<=n;i++)

if(fcmp(fabs(A[i][c]),fabs(A[mxr][c]))>0) mxr=i;

//第c项在第r个方程之后系数都为0

if(fcmp(A[mxr][c],0.0)==0)

{r--;continue;}//执行下一项

//选好主，交换方程

if(mxr!=r)swap(A[mxr],A[r]);

//消元

for(int i=1;i<=n;i++)

if(i!=r&&fcmp(A[i][c],0.0)!=0)

for(int j=m;j>=c;j--)

A[i][j]-=A[r][j]/A[r][c]*A[i][c];

}

Rank=r-1;

}

//判断是否有解

bool check()

{

//判断是否无解

for(int i=Rank+1;i<=n;i++)

if(fcmp(A[i][m],0.0)!=0)//在x系数为0时，结果不为0 return0;

//初始所有都是未知变量

for(int i=1;i

Free[i]=1;

//计算解

for(int i=Rank,cnt=0,pos;i>0;i--)

{

cnt=0;//记数当前方程未知变量个数

for(int j=1;j

if(Free[j]&&fcmp(A[i][j],0.0)!=0) cnt++,pos=j;

if(cnt==1)//一个未知变量，可计算

{

Free[pos]=0;

X[pos]=A[i][m]/A[i][pos];

}

}

return1;

}

int main()

{

freopen("Gauss_data.in","r",stdin);

scanf("%d%d",&n,&m);

m++;

for(int i=1;i<=n;i++)

for(int j=1;j<=m;j++)

scanf("%lf",&A[i][j]);

Gauss();

if(!check())

printf("No solution\n");

else

{

if(Rank

{

printf("Multiple solution, free_num: %d\n",m-1-Rank);

for(int i=1;i

if(Free[i])

printf("X[%d] not determined\n",i);

else

printf("X[%d] = %.4lf\n",i,X[i]);

}

else

for(int i=1;i

printf("X[%d] = %.4lf\n",i,X[i]);

}

return0; }