高斯·约当消元法
说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。
高斯·约当消元法 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+1000000000000x2=10000000000000+999999999999x2=2 而当位数一多,我们的计算机就有可能存不了那么多9,从而忽略掉一堆,原本的误差就扩散得越来越大,这时候我们需要选主消元(比如把刚才两个式子交换一下),选择当前这一项系数绝对值最大的方程作为主元,来消掉其它方程。 高斯·约当消元: 高斯消元是当执行到第i个方程的xkxk时,消掉i以后的xkxk。而约当就同时消掉i以前的xkxk,使系数变成一条对角线 解得情况:高斯消元完成后,下面若有方程系数全部为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) return 1; else if((a-b)>-EPS) return 0; 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 return 0; //初始所有都是未知变量 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]; } }
return 1; }
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(Rank1) {
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]);
}
return 0; }
本文来源:https://www.wddqw.com/doc/8aa3cd6db007e87101f69e3143323968011cf4fa.html