高斯消元法
高斯消元法主要有消元和回代两个过程,和其他求解线性方程组的方法如克莱姆法则,约当消除法相比运算量小,速度快。更适合在计算机内使用。
假设有n个方程,var个变量,时间:
for i=1 --> n-1
for j=i+1-->n
for k=j+1-->n+1
所以时间是O(n^3)
把方程组写成增广矩阵,进行行初等变换,系数矩阵形成上三角矩阵,然后回代逐一求解每一个未知数。
结果又分为:无解;无穷多的解(找出自由变量);浮点数解;唯一整数解等情况
高斯消元的模板:
#include#include#include#include#includeusing namespace std;const int MAXN=50;int a[MAXN][MAXN];//增广矩阵int x[MAXN];//解集bool free_x[MAXN];//不确定的变元int equ,var,free_num;void Debug(){ int i, j; for (i = 0; i < equ; i++) { for (j = 0; j < var + 1; j++) { cout << a[i][j] << " "; } cout << endl; } cout << endl;}int gcd(int a,int b){ return b?gcd(b,a%b):a;}int lcm(int a,int b){ return a/gcd(a,b)*b;}// 高斯消元法解方程组(free_num=-2表示有浮点数解,但无整数解,//-1表示无解,0表示唯一解,大于0表示无穷解,并返回自由变元的个数)//有equ个方程,var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var.int Gauss(int equ,int var){ int i,j,k; int max_r; int col;//当前处理的列 int ta,tb; int LCM; int temp; int free_x_num; int free_index; col=0; for(k = 0;k < equ && col < var;k++,col++) {// 找到该col列元素绝对值最大的那行与第k行交换.(为了处理无解的判断) max_r=k; for(i=k+1;iabs(a[max_r][col])) max_r=i; } if(max_r!=k) { for(j=k;j= 0; i--) { free_x_num = 0; // 用于判断该行中的不确定的变元的个数,如果超过1个,则无法求解,它为不确定的变元. for (j = 0; j < var; j++) { if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j; } if (free_x_num > 1) continue; temp = a[i][var]; for (j = 0; j < var; j++) { if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j]; } x[free_index] = temp / a[i][free_index]; free_x[free_index] = 0; } return var - k; // 自由变元有var - k个. } for (i = var - 1; i >= 0; i--) { temp = a[i][var]; for (j = i + 1; j < var; j++) { if (a[i][j] != 0) temp -= a[i][j] * x[j]; } if (temp % a[i][i] != 0) return -2; // 说明有浮点数解,但无整数解. x[i] = temp / a[i][i]; } return 0;}void in(){ memset(a, 0, sizeof(a)); memset(x,0,sizeof(x)); memset(free_x,1,sizeof(free_x)); int i,j; for (i = 0; i < equ; i++) { for (j = 0; j < var + 1; j++)scanf("%d", &a[i][j]); }}void out(){ int i,j; if (free_num == -1) printf("no solutions!\n"); else if (free_num == -2) printf("here are double but not Integer solutions!\n"); else if (free_num > 0) { printf("here are many solutions and var number is %d\n", free_num); for (i = 0; i < var; i++) { if (free_x[i]) printf("x%d is not certain\n", i + 1); else printf("x%d: %d\n", i + 1, x[i]); } } else for (i = 0; i < var; i++)printf("x%d: %d\n", i + 1, x[i]); printf("\n");}int main(void){ freopen("cin.txt", "r", stdin); int i, j; while (scanf("%d %d", &equ, &var) != EOF) { in(); free_num = Gauss(equ,var); out(); } return 0;}
测试用例:
3 3
2 1 -1 8
-3 -1 2 -11
-2 1 2 -3
3 3
3 2 0 0
1 1 0 -1
0 0 0 0
3 3
3 2 0 0
1 1 1 -1
0 1 1 -2
3 3
3 2 0 0
1 1 0 -1
1 1 0 0
对应的答案是:
2
3
-1
2
-3
任意数
1
-3/2
-1/2
程序的结果: x1: 2 x2: 3 x3: -1 here are many solutions and var number is 1 x1: 2 x2: -3 x3 is not certain here are double but not Integer solutions! no solutions! 中间的过程(以:
3 3 2 1 -1 8 -3 -1 2 -11 -2 1 2 -3 为例): -3 -1 2 -11 2 1 -1 8 -2 1 2 -3 -3 -1 2 -11 0 1 1 2 0 5 2 13 -3 -1 2 -11 0 5 2 13 0 1 1 2 -3 -1 2 -11 0 5 2 13 0 0 3 -3 -3 -1 2 -11 0 5 2 13 0 0 3 -3 -3 -1 2 -11 0 5 2 13 0 0 3 -3 x1: 2 x2: 3 x3: -1
高斯消元过程仅仅涉及行变换,所以得到的X[]解集顺序是不变的。
压缩的写法:
#include #include #include #include using namespace std;const int N=105;int gcd(int a,int b){ return b==0?a:gcd(b,a%b);}// a[][]:增广矩阵 n: 方程数 m:变量数void Gauss(int a[N][N],int n,int m,int &r,int &c){ for(r=0,c=0;rabs(a[maxi][c])) maxi=i; } if(maxi!=r){ for(int j=r;j=0;i--){ int dex=0; int cnt=0; for(int j=0;j1) continue; int t=a[i][m]; for(int j=0;j=0;i--){ int t=a[i][c]; for(int j=i+1;j>n>>m){ for(int i=0;i浮点数型的高斯消元法:
#include #include #include #include #include using namespace std;const int N=105;const double eps=1e-7;int myCmp(double a,double b){ if(fabs(a-b)b)-(a=0) a-=m; while(myCmp(a,0)<0) a+=m; return a;}double gcd(double a,double b){ return myCmp(b,0)==0?a:gcd(b,myMod(a,b));}// a[][]:增广矩阵 n: 方程数 m:变量数void Gauss(double a[N][N],int n,int m,int &r,int &c){ for(r=0,c=0;r0) maxi=i; } if(maxi!=r){ for(int j=r;j=0;i--){ int dex=0; int cnt=0; for(int j=0;j1) continue; double t=a[i][m]; for(int j=0;j=0;i--){ double t=a[i][c]; for(int j=i+1;j>n>>m){ for(int i=0;i用几个式子验证,结果让人满意。
用mathematica计算:
In[1]:= Solve[{0.001*x1 + 2 x2 == 1, 3*x1 + 5*x2 == 2}, {x1, x2}]Out[1]= {{x1 -> -0.166806, x2 -> 0.500083}}In[2]:= Solve[{0.0000000000000001*x1 + 20000*x2 == 10000, 300*x1 + 500*x2 == 2000}, {x1, x2}]Out[2]= {{x1 -> 5.83333, x2 -> 0.5}}In[3]:= Solve[{0.0000000000000000001*x1 + 20000*x2 == 10000, 30000000000000*x1 + 500*x2 == 2000}, {x1, x2}]Out[3]= {{x1 -> 5.83333*10^-11, x2 -> 0.5}}高斯运算结果:-0.166805671392827 0.5000834028356965.833333333333333 0.5000000000000000.000000000058333 0.500000000000000Process returned 0 (0x0) execution time : 0.375 sPress any key to continue.
版权声明:本文内容由网络用户投稿,版权归原作者所有,本站不拥有其著作权,亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容,请联系我们jiasou666@gmail.com 处理,核实后本网站将在24小时内删除侵权内容。
暂时没有评论,来抢沙发吧~