浅谈高斯消元法
2023.6.16:发布
2023.8.29:修缮,加上自己觉得通俗易懂的理解,更新矩阵求逆。
高斯消元
高斯消元可以用于线性方程组求解或者行列式计算,求矩阵的逆等等,也算是比较基础的一个思想。
消元法#
定义#
消元法是将方程组中的一方程的未知数用含有另一未知数的代数式表示,并将其带入到另一方程中,这就消去了一未知数,得到一解;或将方程组中的一方程倍乘某个常数加到另外一方程中去,也可达到消去一未知数的目的。消元法主要用于二元一次方程组的求解。
举例:利用消元法求解二元一次线性方程组:
我们很容易由二式得到:
我们把这个式子代入一式,然后得到:
展开得到:
然后一式的未知数由两个变为了一个,我们可以解出:
然后代入一式或者二式就可以解出来
消元法理论核心#
-
两方程互换,解不变
-
一方程乘以非零数
,解不变 -
一方程乘以数
加上另一方程,解不变(此时的 可以是 )
高斯消元法思想概念#
德国数学家高斯对消元法进行了思考分析,得出了如下结论:
-
在消元法中,参与计算和发生改变的是方程中各变量的系数;
-
各变量并未参与计算,且没有发生改变;
-
可以利用系数的位置表示变量,从而省略变量;
-
在计算中将变量简化省略,方程的解不变。
高斯在这些结论的基础上,提出了高斯消元法,首先将方程的增广矩阵利用行初等变换化为行最简形,然后以线性无关为准则对自由未知量赋值,最后列出表达方程组通解。
高斯消元五步法#
-
增广矩阵行初等变换为行最简形
-
还原线性方程组
-
求解第一个变量
-
补充自由未知量
-
列表示方程组的通解
过程#
我看到 OI Wiki 上写的太过复杂,所以我尽量用通俗易懂的语言来叙述。
例如我们需要解下面的线性方程组:
增广矩阵初等变换为行最简形#
最简形矩阵
定义形如下面矩阵
其判断的依据可以简记为是否能画出一条阶梯状的线将矩阵分为上下两部分,满足在其下面的元素都是
增广矩阵是由方程组系数矩阵
上面的线性方程组可以变换为下面的形式。
我们用第三行减去两倍的第二行,得到:
然后第一行除以
然后我们用第一行减去
至此化为了最简形矩阵
还原线性方程组#
也就是把最简形矩阵重新书写成线性方程组的形式。
求解第一个变量#
也就是将每一行第一个系数不为
补充自由未知量#
第三步中得到的方程式只有
列表示方程组的通解#
其中的
由于
其实可以看做化为一个上三角矩阵然后从下往上依次回代的过程。
我们就是使每一行的对角线上的元素都是
在处理矩阵的时候,我们一般用到以下几种操作:
-
交换两行
-
系数化为
-
某行加上或减去另一行的
倍
首先我们确定第
然后往下遍历,下面行只要第
整个高斯消元的过程:首先我们钦定每一个变量系数为
的行在第 行,我们找到一行 系数不为 的一行,然后再枚举把下面每一行的所有 的系数都化为 ,然后到了最后得到了一个上三角矩阵,这样我们从最后一行,依次把值往上面的方程组中回代,逐步求得所有方程的解。
参考代码:
/*
* @Author: Aisaka_Taiga
* @Date: 2023-08-29 11:02:40
* @LastEditTime: 2023-08-29 16:48:06
* @LastEditors: Aisaka_Taiga
* @FilePath: \Desktop\P3389.cpp
* 心比天高,命比纸薄。
*/
#include<bits/stdc++.h>
#define DB double
#define N 1100
using namespace std;
int n;
DB a[N][N], x[N];
inline int work()
{
for(int i = 1; i <= n; i ++)
{
int r = i;//当前变量
for(int k = i; k <= n; k ++)//找到此列下面第一个不为0的值所在的行
if(fabs(a[k][i]) > 0){r = k; break;}
if(r != i) swap(a[r], a[i]);//交换行
if(fabs(a[i][i]) == 0)return 1;//如果当前变量的系数是0则此变量可以取任意值,无数个解
for(int j = n + 1; j >= i; j --)//将第i个变量所在的方程变成i变量系数为1的方程
a[i][j] /= a[i][i];
for(int k = i + 1; k <= n; k ++)//遍历后面的每一行
for(int j = n + 1; j >= i; j --)//遍历后面的每一列
a[k][j] -= a[k][i] * a[i][j];//把后面所有的第i个变量系数都变成0
}
for(int i = n - 1; i; i --)//遍历每一行
for(int j = i + 1; j <= n; j ++)//遍历每一列
a[i][n + 1] -= a[i][j] * a[j][n + 1];//a[i][j]是系数,a[j][n+1]是当前第j列变量的值,因为j-n的所有变量都算出来了所以直接用
return 0;//有解
}
signed main()
{
cin >> n;
for(int i = 1; i <= n; i ++)
for(int j = 1; j <= n + 1; j ++)
cin >> a[i][j];
int flag = work();
if(flag != 0) return cout << "No Solution" << endl, 0;
for(int i = 1; i <= n; i ++)
printf("%.2lf\n", a[i][n + 1]);
return 0;
}
高斯约旦消元法#
他和普通的高斯消元的区别就是这个最后得到的矩阵是只有对角线是有值的,其余元素为
最后得到的矩阵就像下面这样:
其中
这样做的好处就是让每一个变量都直接乘上自己的系数得到一个值,这样就便于计算了。
最后直接循环一下计算即可。
当然弊端就是无法判断无解是没有解还是无数个解。
高斯约旦消元法和朴素的高斯消元法最大的区别就是最后得到的矩阵的不同,朴素的高斯消元法得到的就是一个上三角矩阵,我们想要得到具体的值还需要倒着一个一个往回带,但是高斯约旦最后得到的是一个对角矩阵,且每一个方程对应一个未知量,系数都为
具体的过程就是,还是钦定第
参考代码
/*
* @Author: Aisaka_Taiga
* @Date: 2023-08-29 16:22:52
* @LastEditTime: 2023-08-29 17:10:20
* @LastEditors: Aisaka_Taiga
* @FilePath: \Desktop\P3389高斯约旦消元法.cpp
* 心比天高,命比纸薄。
*/
#include<bits/stdc++.h>
#define DB double
#define N 1100
using namespace std;
int n;
DB a[N][N];
signed main()
{
cin >> n;
for(int i = 1; i <= n; i ++)
for(int j = 1; j <= n + 1; j ++)
cin >> a[i][j];
for(int i = 1; i <= n; i ++)//? 枚举每一个变量
{
int pl = i;//? 当前行的变量编号
for(int j = i; j <= n; j ++)//? 枚举每一行
if(a[j][i] > a[pl][i]) pl = j;//?找到系数最大的一行
if(a[pl][i] == 0) return cout << "No Solution" << endl, 0;//?如果是零就是无解
swap(a[i], a[pl]);//? 交换两列
DB k = a[i][i];//? 取出当前变量所在行的系数
for(int j = 1; j <= n + 1; j ++)//? 枚举当前行的系数
a[i][j] = a[i][j] / k;//? 系数化为1
for(int j = 1; j <= n; j ++)//? 枚举每一行
{
if(i == j) continue;//? 如果是当前变量所在行就跳过,因为我们处理过了
DB kk = a[j][i];//? 取出当前行的第i个变量的系数
for(int o = i; o <= n + 1; o ++)//? 枚举当前行的每一个元素的系数
a[j][o] = a[j][o] - kk * a[i][o];//? 将除第i行以外的所有的i的行都化为1
}
}
for(int i = 1; i <= n; i ++)
printf("%.2lf\n", a[i][n + 1]);
return 0;
}
高斯消元法对矩阵求逆#
对于一个方阵
求解方法如下:
-
构造一个
的矩阵,左半边是 ,右半边是 。 -
然后对
的左半部分进行高斯消元,最后的右半边矩阵就成了 。
要注意的是,如果左半部分最后得到的不是单位矩阵就说明不存在
下面来看一道模板题:
P4783 【模板】矩阵求逆#
/*
* @Author: Aisaka_Taiga
* @Date: 2023-08-29 17:34:42
* @LastEditTime: 2023-08-29 19:17:43
* @LastEditors: Aisaka_Taiga
* @FilePath: \Desktop\P4783.cpp
* 心比天高,命比纸薄。
*/
#include <bits/stdc++.h>
#define int long long
#define P 1000000007
#define DB double
#define N 1010
using namespace std;
int n, a[N][N << 1];
inline int ksm(int a, int b)
{
int res = 1;
while(b)
{
if(b & 1) res = (res * a) % P;
a = (a * a) % P;
b >>= 1;
}
return res % P;
}
inline int Aisaka_Taiga()
{
for(int i = 1; i <= n; i ++)
{
int p = i;
for(int j = i + 1; j <= n; j ++)
if(a[j][i] > a[p][i]) p = j;
if(a[p][i] == 0) return 1;
if(p != i) swap(a[p], a[i]);
int k = ksm(a[i][i], P - 2);
for(int j = 1; j <= 2 * n; j ++)
a[i][j] = (a[i][j] * k) % P;
for(int j = 1; j <= n; j ++)
{
if(i == j) continue;
int k = a[j][i];
for(int o = i; o <= 2 * n; o ++)
a[j][o] = (a[j][o] - k * a[i][o] % P + P) % P;
}
}
return 0;
}
signed main()
{
cin >> n;
for(int i = 1; i <= n; i ++)
for(int j = 1; j <= n; j ++)
cin >> a[i][j];
for(int i = 1; i <= n; i ++) a[i][n + i] = 1;
int ff = Aisaka_Taiga();
if(ff == 1) return cout << "No Solution" << endl, 0;
for(int i = 1; i <= n; i ++)
{
for(int j = n + 1; j <= 2 * n; j ++)
cout << a[i][j] << " ";
cout << endl;
}
return 0;
}
作者: 北烛青澜
出处:https://www.cnblogs.com/Multitree/p/17426691.html
本站使用「CC BY 4.0」创作共享协议,转载请在文章明显位置注明作者及出处。
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 单线程的Redis速度为什么快?
· 展开说说关于C#中ORM框架的用法!
· Pantheons:用 TypeScript 打造主流大模型对话的一站式集成库
· SQL Server 2025 AI相关能力初探
· 为什么 退出登录 或 修改密码 无法使 token 失效