Loading [MathJax]/jax/output/CommonHTML/autoload/mtable.js

高斯消元 学习笔记

Gaussian EliminationGaussian Elimination



数学上,高斯消元法(或译:高斯消去法),是线性代数规划中的一个算法,可用来为线性方程组求解。

——百度百科

说实话,我不相信这是高斯发明的。感觉像是个小学生都学过的加减消元法。

它的时间复杂度与方程个数、未知数个数有关,一般来讲,是 O(n3)O(n3)

概念1:线性方程组

有多个未知数,且每个未知数的次数皆为一次,这样的多个未知数所组成的方程组被称为 线性方程组

其形式为:

{a1,1x1+a1,2x2++a1,nxn=b1a2,1x1+a2,2x2++a2,nxn=b2an,1x1+an,2x2++an,nxn=bn

概念2:增广矩阵

我们将上面线性方程组的所有系数记录在一个矩阵里,就叫此矩阵为增广矩阵

[a1,1a1,2a1,nb1a2,1a2,2a2,nb2an,1an,2an,nbn]

概念3:主元

算法实现中的当前未知数称为主元

高斯消元法:

例:

{3x+2y+z=6(1)2x+2y+2z=4(2)4x2y2x=2(3)

首先,进行消元操作:

将 (1)/3 ,将 x 的系数化为 1x+23y+13z=2   (1)

再令 (2)2×(1)(3)4×(1),消去 (2),(3)x

得到:

{x+23y+13z=2(1)0x+23y+43z=0(2)0x143y103z=6(3)

然后,令 (2) 除以 23 ,使得 y 的系数成为 1,得到:

y+2z=0  (2)

再令 (3)143×(2),使得 3y 被消去,得到:

{x+23y+13z=2(1)0x+y+2z=0(2)0x+0y+183z=6(3)

(3) 可以清晰地得到 z=1,将其带入 (2)便可快速地解出 y=2

最后将二者带回 (1) 即可清晰快速地得到 x=1 ,所以此线性方程组的解为:

{x=1y=2z=1

增广矩阵模拟:

初态:

[312622244222]

消去 x,得到:

[12313202343001431036]

消去 y,可以得到:

[1231320120001836]

回带解得:

[121]

判断无单解情况

  • 无解

消元完毕后,若有一行系数全部为 0,但最右端常数项不为 0,即无解,例:

[123401200006]

  • 多解 ( 前提:并非无解 )

消元完毕后,若有一行系数全部为 0,但最右端常数也为 0,此时存在多解,例:

[123400000000]

例题:

[SDOI2006] 线性方程组

题目描述

已知 n 元线性一次方程组。

{a1,1x1+a1,2x2++a1,nxn=b1a2,1x1+a2,2x2++a2,nxn=b2an,1x1+an,2x2++an,nxn=bn

请根据输入的数据,编程输出方程组的解的情况。

输入格式

第一行输入未知数的个数 n
接下来 n 行,每行 n+1 个整数,表示每一个方程的系数及方程右边的值。

输出格式

如果有唯一解,则输出解(小数点后保留两位小数)。

如果方程组无解输出 1
如果有无穷多实数解,输出 0

算法实现:

首先,将输入数据存成 double 类型的增广矩阵,然后进行模拟:

  1. 枚举两个变量 r,c 分别表示当前行和当前主元,找到主元系数不为 0 的第 t

  2. 交换当前行和第 t

  3. 把当前行主元系数变为 1

  4. 把从第 r+1 行道第 n 行的当前主元的系数变为 0

  • 最后,从第 n行开始往上解出答案。

  • 若循环完后的 r<=n ,则不存在单解,具体请在代码中看。

#include<bits/stdc++.h>
using namespace std;
#define int long long
int n;
int f;
double a[101][101];
signed main()
{
cin>>n;
for(int i=0;i<n;i++)
{
for(int j=0;j<=n;j++)
{
cin>>a[i][j];
}
}
int r=0,c=0;
for(;c<n;c++)
{
int t=r;
for(int i=r;i<n;i++)
{
if(fabs(a[i][c])>fabs(a[t][c]))
{
t=i;
}
}
if(fabs(a[t][c])<1e-6)
{
continue;
}
for(int i=c;i<=n;i++)
swap(a[r][i],a[t][i]);
/* 主元变为 1 */
for(int i=n;i>=c;i--)
{
a[r][i]/=a[r][c];
}
/* 清理后面的子矩阵 */
for(int i=r+1;i<n;i++)
{
if(fabs(a[i][c])>1e-6)
for(int j=n;j>=c;j--)
{
a[i][j]-=a[i][c]*a[r][j];
}
}
r++;
}
/* 回带三角矩阵求解 */
for(int i=n-1;i>=0;i--)
{
for(int j=i+1;j<n;j++)
{
a[i][n]-=a[i][j]*a[j][n];
// a[i][j]=0;
}
// a[i][i]=1;
}
/* 判断不存在单解的情况 */
if(r<n)
{
for(int i=r;i<n;i++)
{
if(fabs(a[i][n])>1e-6)
{
cout<<-1;return 0;
}
}
cout<<0;return 0;
}
for(int i=0;i<n;i++)
{
printf("x%d=%.2lf\n",i+1,a[i][n]);
}
}
posted @   ccjjxx  阅读(39)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
点击右上角即可分享
微信分享提示