高斯消元解线性方程组
2015-11-17 12:30
423 查看
在线性代数中我们会经常遇到诸如求以下解线性方程组的问题:
x1*a11+x2*a12+...xn*a1n=b1
x1*a21+x2*a22+...xn*a2n=b2
...
x1*am1+x2*am2+...xn*amn=bm
我们可以用高斯消元法来求解上述方程组的解,具体实现代码如下:
整型系数代码如下:
高精度实现代码如下:
x1*a11+x2*a12+...xn*a1n=b1
x1*a21+x2*a22+...xn*a2n=b2
...
x1*am1+x2*am2+...xn*amn=bm
我们可以用高斯消元法来求解上述方程组的解,具体实现代码如下:
整型系数代码如下:
/* Author:Ibsen Date:2015.11.16 */ #include <cstdio> #include <iostream> using namespace std; const int M=100; //处理的矩阵阶数的最大值 int a[M][M]; //系数矩阵和增广矩阵 int x[M]; //存放最后的解 int equ,var; //equ是系数阵的行数,var是系数阵的列数 void Debug() { cout<<"消元结果为:"<<endl; for(int i=0;i<equ;i++) { for(int j=0;j<var+1;j++) cout<<a[i][j]<<" "; cout<<endl; } } int abs(int x) { if(x>=0) return x; else return -x; } inline int Gcd(int a,int b) { return !b?a:Gcd(b,a%b); } inline int Lcm(int a,int b) { return a*b/Gcd(a,b); } void swap(int &a,int &b) { int tmp=a; a=b; b=tmp; } int Gauss() { int k,col=0; //k表示当前处理的行,col表示当前处理的列 for(k=0;k<equ&&col<var;k++,col++) { int max_r=k; for(int i=k+1;i<equ;i++) if(abs(a[i][col])>abs(a[max_r][col])) //这里要找绝对值最大的列,避免代码尾部样例第二组(只有负数和0,最大值为0)的情况的情况 max_r=i; if(max_r!=k) { for(int i=k;i<var+1;i++) swap(a[k][i],a[max_r][i]); } if(!a[k][col]) //(第k行往后)第col列元素全为0,直接进入下一列 { k--; continue; } for(int i=k+1;i<equ;i++) { if(a[i][col]) { int lcm=Lcm(a[i][col],a[k][col]); int ta=lcm/a[i][col], tb=lcm/a[k][col]; for(int j=col;j<var+1;j++) a[i][j]=a[i][j]*ta-a[k][j]*tb; } } } for(int i=k;i<equ;i++) if(a[i][col]) return -1; //无解 Debug(); for(int i=0;i<equ;i++) { if(!a[i][i]) { int j=i+1; for(;j<var;j++) if(a[i][j]) break; if(j==var) break; for(int k=0;k<equ;k++) swap(a[k][i],a[k][j]); } } Debug(); cout<<"col="<<col<<endl; if(var-k==0) //有唯一解:回带求解 { for(int i=k-1;i>=0;i--) //找出最后一个非零行 { int tmp=a[i][var]; for(int j=i+1;j<var;j++) if(a[i][j]) { tmp=tmp-a[i][j]*x[j]; } x[i]=tmp/a[i][i]; } } } void display() { cout<<"解向量X为:"; for(int i=0;i<var;i++) cout<<x[i]<<" "; cout<<endl; } int main() { while(cin>>equ>>var) { for(int i=0;i<equ;i++) for(int j=0;j<var+1;j++) cin>>a[i][j]; if(Gauss()==-1) puts("No solution!"); else display(); } return 0; } /* 2 2 4 5 13 -6 3 -9 解向量为(2,1) 4 4 2 -1 0 0 1 -1 2 -1 0 0 0 -1 2 -1 0 0 0 -1 2 6 解向量为(2,3,4,5) */
高精度实现代码如下:
/* Author:Ibsen Date:2015.11.16 */ #include <cstdio> #include <iostream> using namespace std; const int M=100; //处理的矩阵阶数的最大值 const double esp=1e-8; //精度 double a[M][M]; //系数矩阵和增广矩阵 double x[M]; //存放最后的解 int equ,var; //equ是系数阵的行数,var是系数阵的列数 double abs(double x) { if(x>0) return x; else return -x; } void Debug() { cout<<"消元结果为:"<<endl; for(int i=0;i<equ;i++) { for(int j=0;j<var+1;j++) cout<<a[i][j]<<" "; cout<<endl; } } void swap(double &a,double &b) { double tmp=a; a=b; b=tmp; } int Gauss() { int k,col=0; //k表示当前处理的行,col表示当前处理的列 for(k=0;k<equ&&col<var;k++,col++) { int max_r=k; for(int i=k+1;i<equ;i++) if(abs(a[i][col])>abs(a[max_r][col])) max_r=i; if(max_r!=k) { for(int i=k;i<var+1;i++) swap(a[k][i],a[max_r][i]); } if(abs(a[k][col])<esp) //(第k行往后)第col列元素全为0,直接进入下一列 { k--; continue; } for(int i=k+1;i<equ;i++) { if(abs(a[i][col])>esp) { double tmp=a[k][col]/a[i][col]; for(int j=col;j<var+1;j++) a[i][j]=a[i][j]*tmp-a[k][j]; } } } for(int i=k;i<equ;i++) if(abs(a[i][col])<esp) return -1; //无解 //Debug(); for(int i=0;i<equ;i++) { if(abs(a[i][i])<esp) { int j=i+1; for(;j<var;j++) if(abs(a[i][j])>esp) break; if(j==var) break; for(int k=0;k<equ;k++) swap(a[k][i],a[k][j]); } } Debug(); cout<<"col="<<col<<endl; if(var-k==0) //有唯一解:回带求解 { for(int i=k-1;i>=0;i--) //找出最后一个非零行 { double tmp=a[i][var]; for(int j=i+1;j<var;j++) if(abs(a[i][j])>esp) { tmp=tmp-a[i][j]*x[j]; } x[i]=tmp/a[i][i]; } } } void display() { cout<<"解向量X为:"; for(int i=0;i<var;i++) cout<<x[i]<<" "; cout<<endl; } int main() { while(cin>>equ>>var) { for(int i=0;i<equ;i++) for(int j=0;j<var+1;j++) cin>>a[i][j]; if(Gauss()==-1) puts("No solution!"); else display(); } return 0; } /* 3 3 1.0 2.0 3.0 6.6 2.0 3.0 4.0 9.9 1.0 3.0 2.0 6.6 解向量为(1.1,1.1,1.1) */
相关文章推荐
- Mysql免安装版安装配置及常用操作
- linux中获取系统时间的几种方法
- 二叉树的递归遍历与非递归改写
- FastJson和AsyncHttpCLient
- Android4.2以后,多屏幕的支持 学习(一)
- python 数据结构之-列表
- Python中文编码
- spring aop
- assign,copy,retain的区别以及weak和strong的区别
- 运笔中要强调注意的问题
- 安卓SDK下载国内镜像网站 http://gmirror.org/#android-sdk-tools-only
- NSDate和NSString之间的相互转换
- 六、Maven依赖机制
- jdk环境变量配置
- cocos2d-lua class 方法解释
- 全志科技的一道笔试题
- UIControl中的相关属性和方法
- StringBuffer一些知识
- Win10 TH2正式版10586(1511)免费自动升级的图文教程
- Volley