高斯消元法,是线性代数中的一个算法,可用来求解线性方程组,并可以求出矩阵的秩,以及求出可逆方阵的逆矩阵。

在讲算法前先介绍些概念

矩阵的初等变换

矩阵的初等变换又分为矩阵的初等行变换和矩阵的初等列变换。矩阵的初等行变换和初等列变换统称为初等变换。另外:分块矩阵也可以定义初等变换。

等价

定义:如果B可以由A经过一系列初等变换得到,则称矩阵A与B称为等价

初等行变换

定义:所谓数域P上矩阵的初等行变换是指下列3种变换:
1)以P中一个非零的数乘矩阵的某一行
2)把矩阵的某一行的c倍加到另一行,这里c是P中的任意一个数
3)互换矩阵中两行的位置
一般来说,一个矩阵经过初等行变换后就变成了另一个矩阵,当矩阵A经过初等行变换变成矩阵B时,一般写作

 
可以证明:任意一个矩阵经过一系列初等行变换总能变成阶梯型矩阵。(方阵的话就是上三角矩阵)

初等列变换

同样地,定义初等列变换,即:
1)以P中一个非零的数乘矩阵的某一列
2)把矩阵的某一列的c倍加到另一列,这里c是P中的任意一个数
3)互换矩阵中两列的位置

初等矩阵

定义:由单位矩阵E经过一系列初等变换得到的矩阵称为初等矩阵。
引理:对一个

  

矩阵A作一初等行变换,就相当于在A的左边乘上相应的

  

的初等矩阵;

对A作一初等列变换,就相当于在A的右边乘上相应的 的初等矩阵。
 
 

高斯消元的思想——消元法(加减消元 & 代入消元)

消元法是将方程组中的一方程的未知数用含有另一未知数的代数式表示,并将其代人到另一方程中,这就消去了一未知数,得到一解;或将方程组中的一方程倍乘某个常数加到另外一方程中去,也可达到消去一未知数的目的。消元法主要用于二元一次方程组的求解。

步骤

首先要把方程转化为增广矩阵矩阵,比如
将上面的方程转为下面的增广矩阵

其中第4列是常数项,其它三列是方程的系数。

然后开始进行消元

高斯消元

首先找到一个第一列的数不为0的行(一般找第一列的数最大的行)(如果都为0就跳过当前步)

然后用它的第一列的数将下面行当前列的值化为0,变换过的初等矩阵与原矩阵等价,化为方程后依然成立

本矩阵第一列的数最大的为第三行就把第三行与第一行交换

然后下面行的当前列消去,如图

除了最后一列外,每一列都如此,最后得到上三角矩阵如图

这样我们很容易算出x3的值,再用x3的值算出x2,x1的值

int gauss(int n){//a为增广矩阵
int r,w;
for(int i=;w<n&&i<n;i++,w++){//进行到第i列,第w行
int r=w;
for(int j=w+;j<n;j++)if(fabs(a[j][i])>fabs(a[r][i]))r=j;//找到当前列绝对值最大的行
if(fabs(a[r][i])<eps){w--;continue;}//当前列值都为0,跨过当前步
if(w!=r)for(int j=i;j<=n;j++)swap(a[r][j],a[w][j]);//交换当前列绝对值最大的行和没计算过的第一行
for(int k=w+;k<n;k++){//消去当前列(只消下面行的)
double pro=a[k][i]/a[w][i];
for(int j=i;j<=n;j++)
a[k][j]-=pro*a[w][j];
}
}
return w;
}

高斯约旦消元

同高斯消元大体差不多,但消元时上面计算过的行也要消去当前列,最后得到的是对角矩阵而不是上三角矩阵。

第一次和高斯消元相同

第二次要顺带把第一行第二列消去

最后将最后一行上面所有行的倒数第二列消去

int gauss_jordan(int n){//a为增广矩阵
int r,w=;
for(int i=;i<n&&w<n;w++,i++){//进行到第i列,第w行
int r=w;
for(int j=w+;j<n;j++)if(fabs(a[j][i])>fabs(a[r][i]))r=j;//找到当前列绝对值最大的行
if(fabs(a[r][i])<eps){w--;continue;}//当前列值都为0,跨过当前步
if(r!=w)for(int j=;j<=n;j++)swap(a[r][j],a[w][j]);//交换当前列绝对值最大的行和没计算过的第一行
for(int k=;k<n;k++){//消去当前列(除本行外)
if(k!=w)
for(int j=n;j>=w;j--)a[k][j]-=a[k][i]/a[w][i]*a[w][j];
}
}
return w;
}

例题

[SDOI2007] 线性方程组

代码

#include<cstdio>
#include<cmath>
#include<algorithm>
#define eps 1e-8
using namespace std;
double a[][],ans[];
int d;
int gauss_jordan(int n){//a为增广矩阵
int r,w=;
for(int i=;i<n&&w<n;w++,i++){//进行到第i列,第w行
int r=w;
for(int j=w+;j<n;j++)if(fabs(a[j][i])>fabs(a[r][i]))r=j;//找到当前列绝对值最大的行
if(fabs(a[r][i])<eps){w--;continue;}//当前列值都为0,跨过当前步
if(r!=w)for(int j=;j<=n;j++)swap(a[r][j],a[w][j]);//交换当前列绝对值最大的行和没计算过的第一行
for(int k=;k<n;k++){//消去当前列(除本行外)
if(k!=w)
for(int j=n;j>=w;j--)a[k][j]-=a[k][i]/a[w][i]*a[w][j];
}
}
return w;
}
int gauss(int n){//a为增广矩阵
int r,w;
for(int i=;w<n&&i<n;i++,w++){//进行到第i列,第w行
int r=w;
for(int j=w+;j<n;j++)if(fabs(a[j][i])>fabs(a[r][i]))r=j;//找到当前列绝对值最大的行
if(fabs(a[r][i])<eps){w--;continue;}//当前列值都为0,跨过当前步
if(w!=r)for(int j=i;j<=n;j++)swap(a[r][j],a[w][j]);//交换当前列绝对值最大的行和没计算过的第一行
for(int k=w+;k<n;k++){//消去当前列(只消下面行的)
double pro=a[k][i]/a[w][i];
for(int j=i;j<=n;j++)
a[k][j]-=pro*a[w][j];
}
}
return w;
}
int main(){
freopen("gaess.in","r",stdin);
// freopen("gaess.out","w",stdout);
int n;scanf("%d",&n);
for(int i=;i<n;i++){
for(int j=;j<=n;j++){
scanf("%lf",a[i]+j);
}
}
#if 0 //guass
d=gauss(n);d--;
for(int j=d;j<n;j++){
if(fabs(a[j][n-])<eps&&fabs(a[j][n])>eps){printf("-1");return ;}//最后一个方程未知数最少,方程=右边不为0=左边为0,则无解
}
if(d<n-){printf("");return ;}//一个方程解一个未知数,有效方程少于n个,则多个解
for(int i=d;i>=;i--){
for(int k=i+;k<n;k++)a[i][n]-=a[i][k]*ans[k];
ans[i]=a[i][n]/a[i][i];
}
#endif
#if 1 //guass_jordan
d=gauss_jordan(n);d--;
for(int i=;i<n;i++){//有一个方程=右边不为0=左边为0,则无解
bool d=;
for(int j=i;j<n;j++)d&=(fabs(a[i][j])<eps);
if(d&&fabs(a[i][n])>eps){
printf("-1");return ;
}
}
for(int i=;i<n;i++){//消元后有变量在多个方程中出现,则有多个解
int max1=;
for(int j=i;j<n;j++)if(fabs(a[i][j])>eps)max1++;
if(max1>){printf("");return ;}
}
for(int i=;i<n;i++){
ans[i]=a[i][n]/a[i][i];
}
#endif
for(int i=;i<n;i++){
if(fabs(ans[i])<eps)printf("x%d=0\n",i+);
else printf("x%d=%.2lf\n",i+,ans[i]);
}
return ;
}

最新文章

  1. HTTP状态码302、303和307的故事
  2. 8.mvc core上传文件
  3. JAVA小知识
  4. 更改apache网站根目录导致localhost不能访问
  5. Java中的Set, List, Map漫谈
  6. java 21 - 11 IO流的标准输入流和标准输出流
  7. Java迭代 : Iterator和Iterable接口
  8. android Activity的启动模式
  9. Emmet(前身是zen coding)介绍和使用
  10. 分布式文件系统 FastDFS Ceph
  11. [RxJS] Adding Conditional Logic with Filter
  12. Android - 位置定位(Location)服务(Service)类的基本操作
  13. iOS之NSDictionary和NSArray以及NSMutableDictionary和NSMutableArray:将不再是问题
  14. 跨进程通信之Messenger
  15. [HDU4669]Editor (栈)
  16. hibernate 和mybatis
  17. 【传输协议】什么是CA证书
  18. centos上安装elasticsearch 5.5.1 遇到的各种坑
  19. 【深入理解javascript】this的用法
  20. 配置mysql允许远程链接

热门文章

  1. vue导航条选中项样式
  2. Java编译与反编译命令记录
  3. Pyinstaller打包Web项目
  4. 认识DOM 文档对象模型DOM(Document Object Model)定义访问和处理HTML文档的标准方法。元素、属性和文本的树结构(节点树)。
  5. 锋利的Jquery(Table,Checkbox)
  6. docker 个人遇到问题日志记录
  7. Type.GetType(string.contains(&#39;,&#39;))
  8. 洛谷P3376【模板】网络最大流 ISAP
  9. VS 断点不会命中的情况
  10. PAT甲级——A1080 Graduate Admission