给定线性方程组的系数,求解方程组是否有解。
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;
}