高斯消元求多元一次方程(讲解+板子)
一、我们解多元一次方程需要什么?
因为未知数有多个,所以我们需要方程的数量也不同
如果你要求n元一次方程,那么你至少需要给出n个方程才可能会求出来所有未知量的大小
二、增广矩阵是什么?
我们求多元一次方程的方法就是按照大学线性代数课程中的方法。
我们首先需要构造一个增广矩阵,然后我们把这个增广矩阵化为上三角矩阵,如果顺利的话就可以求出来所有未知数的解,但也有可能只能求出来一些未知数的解。无解情况也会处理。
1、增广矩阵:就是在系数矩阵的右边添上一列,这一列是线性方程组的等号右边的值。(下面给出n个未知量,m个方程构成的增广矩阵)
a11x1+a12x2+a13x3+...+a1nxn=b1
a21x1+a22x2+a23x3+...+a2nxn=b2
......
am1x1+am2x2+am3x3+...+amnxn=bm
2、上三角矩阵:由增广矩阵经过一系列初等变化而成。(最终形式如下,以4个未知数,4个方程为例)
a11 a12 a13 a14
0 a22 a23 a24
0 0 a33 a34
0 0 0 a44
就是第一行第一个元素下面的值都要为0、第二行第二个元素下面的值都要为0、第三行第三个元素下面的值都要为0......
简单说一下初等变化中比较实用的两个规则(使用初等变化之后矩阵所有性质都是不改变的):
1、矩阵中任意两行或者两列可以交换
2、矩阵中一行中每一个元素,可以加上另一行对应的位置k倍的元素,也就是
a11+k*a21 a12+k*a22 ... a1n+k*a2n 可以替换掉下面这一行(k可以取任意值)
a11 a12 a1n
三、矩阵的秩
就是将矩阵转化为上三角矩阵之后不是全0行的行数。(全0行,也就是矩阵这一行所有的值都是0)
示例(R就代表矩阵的秩):
上图截取自百度文库:点我跳转
四、什么是列主元消去法
列主元消去法是在系数矩阵中按列选取元素绝对值得最大值作为主元素(可以通过交换行的方式来满足这个需求)
这样做的原因就是可以减小误差,因为求解方程的时候不能避免出现相除的情况,而计算机内部不能保存分数,也就导致了结果和预想的不一样。
例如,你让变量a=1/3,b=3.你输出a*b的结果就不是1
代码:
col=0; // 当前处理的列 for(k = 0;k < equ && col < var;k++,col++) {// 枚举当前处理的行. // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差) max_r=k; for(i=k+1;i<equ;i++) { if(abs(a[i][col])>abs(a[max_r][col])) max_r=i; } if(max_r!=k) {// 与第k行交换. for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]); } }
五、怎么把系数矩阵转换成上三角矩阵
首先,按照列主元消去法把这一列中最大元素那一行和第一行交换。之后让其他行的第一个元素变成0,那么就会满足公式
a11*k+ay1=0 或者 a11+k*ay1=0 (y>1)
k=-ay1/a11 或者 k=-a11/ay1
这样的话就让第y行每一个元素加上a11*k 或者 ay1*k (这就看你选择那个公式了,其实都差不多)
代码:
for(i=k+1;i<equ;i++) {// 枚举要删去的行. if(a[i][col]!=0) { LCM = lcm(abs(a[i][col]),abs(a[k][col])); ta = LCM/abs(a[i][col]); tb = LCM/abs(a[k][col]); if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加 for(j=col;j<var+1;j++) { a[i][j] = a[i][j]*ta-a[k][j]*tb; } } }
六、判断矩阵的解
如果你是n元一次方程,而且你给了n个方程,且这n个方程化成的上三角矩阵每一行都至少有一个元素不为0(要排除0 0 ... 0 0 a,a!=0,因为着一种情况就代表着每一个系数都是0,但是结果却不是0,这样的情况是无解的)
那么这样的方程n个未知数都可以求出来
如果你有n个未知数,但是你只给出来m各方程(m<n),那么这个注定是求不出来所有未知数的解的(可能会求出来一部分未知数)
如果有k个未知数没有求出来,那就以为着有k-1个变元,就比如x+y=1,如果y确定了,那么x也就确定了。所以变元就只有一个(从这里我们也看出来了,如果存在变元,那么这个方程的解也就有无数个)
代码:
// 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解,但无整数解, //-1表示无解,0表示唯一解,大于0表示无穷解,并返回自由变元的个数) //有equ个方程,var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var. // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0). for (i = k; i < equ; i++) { // 对于无穷解来说,如果要判断哪些是自由变元,那么初等行变换中的交换就会影响,则要记录交换. if (a[i][col] != 0) return -1; } // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行,即说明没有形成严格的上三角阵. // 且出现的行数即为自由变元的个数. if (k < var) { // 首先,自由变元有var - k个,即不确定的变元至少有var - k个. for (i = k - 1; i >= 0; i--) { // 第i行一定不会是(0, 0, ..., 0)的情况,因为这样的行是在第k行到第equ行. // 同样,第i行一定不会是(0, 0, ..., a), a != 0的情况,这样的无解的. 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; // 无法求解出确定的变元. // 说明就只有一个不确定的变元free_index,那么可以求解出该变元,且该变元是确定的. 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个. } // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵. // 计算出Xn-1, Xn-2 ... X0. 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;
七、全部代码
#include<stdio.h> #include<algorithm> #include<iostream> #include<string.h> #include<math.h> using namespace std; const int MAXN=50; int a[MAXN][MAXN];//增广矩阵 int x[MAXN];//解集 bool free_x[MAXN];//标记是否是不确定的变元 /* void Debug(void) { int i, j; for (i = 0; i < equ; i++) { for (j = 0; j < var + 1; j++) { cout << a[i][j] << " "; } cout << endl; } cout << endl; } */ inline int gcd(int a,int b) { int t; while(b!=0) { t=b; b=a%b; a=t; } return a; } inline int lcm(int a,int b) { return a/gcd(a,b)*b;//先除后乘防溢出 } // 高斯消元法解方程组(Gauss-Jordan elimination).(-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; for(int i=0; i<=var; i++) { x[i]=0; free_x[i]=true; } //转换为阶梯阵. col=0; // 当前处理的列 for(k = 0; k < equ && col < var; k++,col++) { // 枚举当前处理的行. // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差) max_r=k; for(i=k+1; i<equ; i++) { if(abs(a[i][col])>abs(a[max_r][col])) max_r=i; } if(max_r!=k) { // 与第k行交换. for(j=k; j<var+1; j++) swap(a[k][j],a[max_r][j]); } if(a[k][col]==0) { // 说明该col列第k行以下全是0了,则处理当前行的下一列. k--; continue; } for(i=k+1; i<equ; i++) { // 枚举要删去的行. if(a[i][col]!=0) { LCM = lcm(abs(a[i][col]),abs(a[k][col])); ta = LCM/abs(a[i][col]); tb = LCM/abs(a[k][col]); if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加 for(j=col; j<var+1; j++) { a[i][j] = a[i][j]*ta-a[k][j]*tb; } } } } // Debug(); // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0). for (i = k; i < equ; i++) { // 对于无穷解来说,如果要判断哪些是自由变元,那么初等行变换中的交换就会影响,则要记录交换. if (a[i][col] != 0) return -1; } // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行,即说明没有形成严格的上三角阵. // 且出现的行数即为自由变元的个数. if (k < var) { // 首先,自由变元有var - k个,即不确定的变元至少有var - k个. for (i = k - 1; i >= 0; i--) { // 第i行一定不会是(0, 0, ..., 0)的情况,因为这样的行是在第k行到第equ行. // 同样,第i行一定不会是(0, 0, ..., a), a != 0的情况,这样的无解的. 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; // 无法求解出确定的变元. // 说明就只有一个不确定的变元free_index,那么可以求解出该变元,且该变元是确定的. 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个. } // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵. // 计算出Xn-1, Xn-2 ... X0. 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; } int main(void) { // freopen("in.txt", "r", stdin); // freopen("out.txt","w",stdout); int i, j; int equ,var; while (scanf("%d %d", &equ, &var) != EOF) { memset(a, 0, sizeof(a)); for (i = 0; i < equ; i++) { for (j = 0; j < var + 1; j++) { scanf("%d", &a[i][j]); } } // Debug(); int free_num = Gauss(equ,var); if (free_num == -1) printf("无解!\n"); else if (free_num == -2) printf("有浮点数解,无整数解!\n"); else if (free_num > 0) { printf("无穷多解! 自由变元个数为%d\n", free_num); for (i = 0; i < var; i++) { if (free_x[i]) printf("x%d 是不确定的\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"); } return 0; }
八、样例输入
示例:求方程x+y=1和x-y=1的解
2方程 2未知数
系数矩阵:
1 1 1
1 -1 1