即对增广矩阵进行初等行变换,将其转换为行最简阶梯形矩阵
AX = B与C(行最简阶梯形
)X = D是同解方程组。
来看算法实现过程:
有矩阵
\begin{pmatrix} 1 & 2 & -1 & -6\\ 2 & 1 & -3 & -9\\ -1& -1 & 2 & 7 \end{pmatrix}
我们来进行初等行变换
一、先来枚举每一列:
①枚举第一列
1.找到第一列中绝对值最大的一行,这里即2 1 -3 -9这行
2.将绝对值最大的一行与第一行交换得到\begin{pmatrix} 2 & 1 & -3 & -9\\ 1 & 2 & -1 & -6\\ -1& -1 & 2 & 7 \end{pmatrix}
3.将第一行第一个数变为1,得到\begin{pmatrix} 1 & 1/2 & -3/2 & -9/2\\ 1 & 2 & -1 & -6\\ -1& -1 & 2 & 7 \end{pmatrix}
4.将下面所有行的第一列消成0,得到\begin{pmatrix} 1 & 1/2 & -3/2 & -9/2\\ 0 & 3/2 & 1/2 & -3/2\\ 0& -1/2 & 1/2 & 5/2 \end{pmatrix}
这时候第一行已经固定
,所以接下来我们要从第二行开始的其他不固定的方程中寻找到第c列中绝对值最大的行。
②枚举第二列
1.找到第二列中绝对值最大的一行(注意此时第一行已经固定,因此我们只能从第二行开始寻找)即0 3/2 1/2 -3/2这行
2.将绝对值最大的一行与第二行交换得到\begin{pmatrix} 1 & 1/2 & -3/2 & -9/2\\ 0 & 3/2 & 1/2 & -3/2\\ 0& -1/2 & 1/2 & 5/2 \end{pmatrix}
3.将第二行的第二个数变为1,得到\begin{pmatrix} 1 & 1/2 & -3/2 & -9/2\\ 0 & 1 & 1/3 & -1\\ 0& -1/2 & 1/2 & 5/2 \end{pmatrix}
4.将第二行之下所有行的第二列消成0,得到\begin{pmatrix} 1 & 1/2 & -3/2 & -9/2\\ 0 & 1 & 1/3 & -1\\ 0& 0 & 2/3 & 2 \end{pmatrix}
这时候第二行已经固定
。
③枚举第三列
这里简化过程得到的结果是\begin{pmatrix} 1 & 1/2 & -3/2 & -9/2\\ 0 & 1 & 1/3 & -1\\ 0& 0 & 1 & 3 \end{pmatrix}
二、将每行首非零元所在列的其他元素化为0
得到的是\begin{pmatrix} 1 & 0 & 0 & 1\\ 0 & 1 & 0 & -2\\ 0& 0 & 1 & 3 \end{pmatrix}
代码实现即模拟手动进行初等行变换
将矩阵化为行阶梯形矩阵
注意:代码实现时,将行阶梯形矩阵
化为行最简阶梯形矩阵
的过程可以跳过,每个解直接用已经得到的解求解出来即可
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const int N = 110;
const double eps = 1e-6;
double a[N][N];
int n;
//debug
void out()
{
for(int i = 1; i <= n; i++){
for(int j = 1; j <= n + 1; j++)
printf("%5.2lf ", a[i][j]);
printf("\n");
}
printf("\n\n");
}
int gauss()
{
/****************************************************************************************************
具体实现:
1.枚举列,找到每列最大的绝对值,将其置换到当前未固定的最小行,然后将最小行之下的所有该列元素全部换为0
2.如果有唯一解,那么再按照从最高行的位置将每一个解求解出来。
****************************************************************************************************/
int r, c; //r为行,c为列
for(c = 1, r = 1; c <= n; c++){
int t = r; //初始化每次未固定的最小行为r
for(int i = r + 1; i <= n; i++)
if(fabs(a[i][c]) > fabs(a[t][c]))
t = i;
//如果这一列都为0,那么根据矩阵变换就无需进行这次的之后的操作了
if(fabs(a[t][c]) < eps) continue;
//如果最小行不是初始化的值
if(t != r){
for(int j = 1; j <= n + 1; j++) swap(a[t][j], a[r][j]);
}
//交换后将其变为1
for(int j = n + 1; j >= c; j--) a[r][j] /= a[r][c];
//把最小行之下的所有该列元素都置为0,同时置为0的列元素所在行的元素也要改变
for(int i = r + 1; i <= n; i++){
if(fabs(a[i][c]) > eps)
for(int j = n + 1; j >= c; j--)
a[i][j] = a[i][j] - a[i][c] * a[r][j];
}
r++;
}
if(r < n + 1){
//由于r < n,所以可知至少有一个解是0x=0,故右边必须是0,否则就是无解
for(int i = r; i <= n; i++)
if(fabs(a[i][n + 1]) > eps)
return 2;
return 1;
}
//当是唯一解时,我们需要将每个解都求解出来,因此我们需要从最后一个解开始,逐个借助已求出的解来解出其他的解
//由于最后一个解本就确定了,故我们从倒数第二个开始解
for(int i = n - 1; i >= 1; i--){
//每次是从该元素的后一个元素开始解
for(int j = i + 1; j <= n; j++)
a[i][n + 1] = a[i][n + 1] - a[i][j] * a[j][n + 1];
}
return 0;
}
int main()
{
scanf("%d", &n);
for(int i = 1; i <= n; i++)
for(int j = 1; j <= n + 1; j++)
scanf("%lf", &a[i][j]);
int t = gauss();
//0代表唯一解,1代表无穷多组解,2代表无解
if(t == 0){
for(int i = 1; i <= n; i++) printf("%.2lf\n", a[i][n + 1]);
}
else if(t == 1) puts("Infinite group solutions");
else puts("No solution");
return 0;
}
总体步骤:
当前已经到了第r行,第c列
1. 那么先找到[r,n]行的最大c列值,即max(a[r~n][c])
2. 将最大c列值所在行和r行交换
3. 将a[r][c]改成1,相应的r行其他元素都得除a[r][c]
4. 将a[r+1~n][c]改成0,相应的对于x$\in$[r+1,n],a[x][c+1~n] -= a[x][c]*a[r][c+1~n]
5. 当r行全部结束后,依次将第n行,n-1行消成除了a[n][n]外其余都为0
参考:y总标程-> https://www.acwing.com/activity/content/code/content/53389/