标签:gauss
给定线性方程组的系数,求解方程组是否有解。
1,找到系数矩阵行列为(k, k)的块绝对值最大的数作为主元,记下行和列,分别与第k列交换,与第k行交换,(行跟行交换,列跟列交换)。
2,交换后,主元所在的行每一个元除以主元的值,使得主元所在的位置为1,常数列也除以主元的值。
3,进行初等行变换,使得主元所在的列的其他元为0。
4,判断系数矩阵和增广矩阵的秩是否相同,相同,有解,不相同,无解。
下面的是代码:
#include <cmath> #include <iostream> #include <iomanip> using namespace std; const int Max = 20; int s[Max]; //变量位置 double x[Max]; //解向量 double a[Max][Max], B[Max]; //系数矩阵和常数列 int rankB(double *b, int m) //常数列中非0的数的个数 { int i = m - 1; while(b[i] == 0) i--; return i; } int guass(int m, int n, double a[][Max], double b[]) { int js[Max], l, k, i, j, is, p; double d, t; for(j = 0; j < n; j++) //变量原始位置 s[j] = j; for(l = 1, k = 0; k <= m - 1; k++) //从第k+1行开始 { d = 0.0; for(i = k; i <= m - 1; i++) for(j = k; j <= n - 1; j++) { t = fabs(a[i][j]); if(t > d) //去最大元素和行列位置 { d = t; js[k] = j; is = i; } } if(d + 1.0 == 1.0) l = 0; else { if(js[k] != k) //不在处理块的第k列,将最大数所在的列调换到第k列 { for(i = 0; i <= m - 1; i++) { t = a[i][k]; a[i][k] = a[i][js[k]]; a[i][js[k]] = t; } p = s[k]; s[k] = s[js[k]]; s[js[k]] = p; //变量也作相应变动 } if(is != k) //不在处理块的第k行,将最大数所在的列调换到第k行 { for(j = 0; j <= n - 1; j++) { t = a[k][j]; a[k][j] = a[is][j]; a[is][j] = t; } t = b[k]; b[k] = b[is]; b[is] = t; } } d = a[k][k]; //(k,k)位置的元素最大 for(j = k; j <= n - 1; j++) a[k][j] = a[k][j] / d; b[k] = b[k] / d; for(i = k + 1; i <= m - 1; i++) //消去过程 { if(a[i][k] != 0) { for(j = k + 1; j <= n - 1; j++) a[i][j] = a[i][j] - a[i][k] * a[k][j]; b[i] = b[i] - a[i][k] * b[k]; a[i][k] = 0; } } } k = m > n ? n : m; //确定矩阵的秩 while(a[k - 1][k - 1] == 0) k--; if((rankB(b, m) > k - 1)) //增广矩阵与系数矩阵秩不同,无解 return 0; d = a[k - 1][k - 1]; for(j = n - 1; j > k - 1; j--) x[s[j]] = 0; for(i = k - 1; i >= 0; i--) { t = 0.0; for(j = i + 1; j <= n - 1; j++) t = t + a[i][j] * x[s[j]]; x[s[i]] = b[i] - t; } return 1; } int main() { int i, j, m, n, cnt = 1; while(cin >> m >> n) { cout << "Case" << cnt++ << endl; for(i = 0; i < m; i++) { for(j = 0; j < n; j++) cin >> a[i][j]; cin >> B[i]; //常数列 } if(guass(m, n, a, B)) { for(j = 0; j < n; j++) cout << "x[" << j << "]=" << setiosflags(ios::fixed) << setprecision(6) << x[j] << ' '; } else cout << "No solutions"; cout << endl; } return 0; }
版权声明:本文为博主原创文章,未经博主允许不得转载。
标签:gauss
原文地址:http://blog.csdn.net/qq_25425023/article/details/47831719