初等变换法求解线性方程组

简介: 初等变换法求解线性方程组

初等变换:通过初等行(列同理)变换把增广矩阵变为简化阶梯型矩阵的线性方程组求解算法

具体步骤:
枚举每一列,找到枚举的当前列绝对值最大数的所在行
将该行换到最上面一行(第r行)
将该行第一个数(该行第c列元素)消成1
将该行第一个元素(该行第c列元素)下面的所有非零元素消成0
重复上面的操作,直到枚举到最后一列,把矩阵消成简化阶梯型矩阵
目录

1.初等变换基本操作

2.具体函数实现

枚举每一列

找到枚举当前列绝对值最大的数的所在行 (ps:fabs函数:求绝对值函数)

将绝对值最大的所在行换到最上面一行(第r行)

将该行第一个数消成1

将该行第一个元素(该行第c列元素)下面的所有非零元素消成0

整体代码

1.初等变换基本操作

   复习:

用一个非零的数乘某一行
把其中一行的若干倍加到另一行上
交换两行的位置
在函数实现的时候无形之中会用到如上的三条操作,以下就不再证明。

2.具体函数实现
枚举每一列
for (c = 0, r = 0; c < n; c++)
找到枚举当前列绝对值最大的数的所在行 (ps:fabs函数:求绝对值函数)
int t = r;

for (int i = r + 1; i < n; i++)
    if (fabs(a[t][c]) < fabs(a[i][c]))
        t = i;

如果该列全为0,continue直接让c++(跳过这一列,开始处理下一列)

ps:由于浮点数加减乘除时候会有精度问题,所以这里认定:如果浮点数 < eps(eps是一个很小的数,其值为1e-6),就认为该浮点数的值为0

< eps 代表 == 0 ,> eps 代表 非零

if (fabs(at) < eps) continue;
将绝对值最大的所在行换到最上面一行(第r行)
for (int i = c; i <= n; i++) swap(at, ar);
将该行第一个数消成1
注意:这里需要倒着消,最后在处理第一个数 (如果先处理第一个数,后面的数就都是除1)

for (int i = n; i >= c; i--) ar /= ar;
将该行第一个元素(该行第c列元素)下面的所有非零元素消成0
for (int i = r + 1; i < n; i++)

if (fabs(a[i][c]) > eps) // 保证消的都是非零元素
    for (int j = n; j >= c; j--)
        a[i][j] -= a[i][c] * a[r][j];

最后让r++,代表这一行已经处理完了,不再参与后面的运算(非常重要)

r++;
进行完如上操作之后,该矩阵就变成了阶梯型矩阵(共有如下几种情况)

1.完美的阶梯型矩阵

例如:

最右侧一列为常数列

由最后一行可以直接得出x3 == 1

将x3往上代入就可以依次求解出x2和x1

可以将向上代入得过程等价于把矩阵化成简化阶梯型矩阵(类对角形矩阵)

此时常数列的值对应的就是x1 x2 x3的值(如图 x1 == 3 ,x2 == 1, x3 == 1)

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];

此处代码非常抽象,下面用图解的方法帮助理解

从下往上消,将阶梯型矩阵消成简化阶梯型矩阵(类对角型矩阵)

由于x3的系数是1,所以1上面的2(ai)就是要乘的倍数

ai -= ai * aj

以此类推

2.不完美的阶梯型矩阵

x的取值一共有两种情况:0或非0

在讨论之前先了解一个概念:

方程个数 < 未知数的个数 -> 此方程有无穷多组解

方程个数 = 未知数的个数 -> 此方程有唯一解

矩阵化简得过程出现了等式左右不相等的情况 -> 此方程无解

若x == 0 ,最后一行相当于是一个恒等式0 == 0 ,此时方程个数 < 未知数的个数,方程有无穷多组解

若x != 0, 最后一行出现了等式左右不相等的情况,此方程无解

if (r < n)//不完美的阶梯型矩阵

{
    for (int i = r; i < n; i++)
        if (fabs(a[i][n]) > eps)
            return 2; // 无解
    return 1; // 无穷多组解
}

//r == n 的情况 往上代入把矩阵化为简化阶梯型矩阵

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];
return 0;//唯一解

返回2代表无解,返回1代表无穷多组解,返回0代表唯一解

整体代码

include

include

include

using namespace std;
const int N = 110;
double aN;
int n;
const double eps = 1e-6;
int gauss()
{

int c, r;
for (c = 0, r = 0; c < n; c++)
{
    int t = r;
    for (int i = r + 1; i < n; i++)
        if (fabs(a[t][c]) < fabs(a[i][c]))
            t = i;
    if (fabs(a[t][c]) < eps) continue;
    for (int i = c; i <= n; i++)    swap(a[t][i], a[r][i]);
    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]) > eps)
            for (int j = n; j >= c; j--)
                a[i][j] -= a[i][c] * a[r][j];
    r++;
}
if (r < n)
{
    for (int i = r; i < n; i++)
        if (fabs(a[i][n]) > eps)
            return 2; // 无解
    return 1; // 无穷多组解
}
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];
return 0;

}
int main()
{

cin >> n;
for (int i = 0; i < n; i++)
    for (int j = 0; j < n + 1; j++)
        cin >> a[i][j];
int t = gauss();

if (t == 0)
{
    for (int i = 0; i < n; i++)
    {
        if (fabs(a[i][n]) < eps) a[i][n] = 0;
        printf("%.2lf\n", a[i][n]);
    }
}
else if (t == 1) puts("Infinite group solutions"); // 无穷多组解
else puts("No solution");//无解
return 0;

}

示例:

结果x1 == 1,x2 == -2,x3 == 3

相关文章
|
4月前
【数值分析】Jacobi、Seidel和Sor迭代法求解线性方程组(附matlab代码)
【数值分析】Jacobi、Seidel和Sor迭代法求解线性方程组(附matlab代码)
|
Linux API iOS开发
|
人工智能 开发者
求解拉格朗日乘子法 | 学习笔记
快速学习求解拉格朗日乘子法
133 0
求解拉格朗日乘子法 | 学习笔记
雅克比迭代法求解线性方程组
雅克比迭代法求解线性方程组
113 0
|
Python
使用Numpy求解方程组
使用Numpy求解方程组
147 0
使用Numpy求解方程组
|
算法
F#实现Runge–Kutta算法求解常微分方程
不少工程问题中涉及的微分方程,我们很难求出方程的解析解,或者说根本不存在精确的解析解。此时,我们需要利用电脑,结合数值分析的方法来近似求出微分方程的相关解,并研究其性质。通过求出多个自变量的值,并求出对应的解,那么可以绘制出图形来辅助研究方程的特征。本文将介绍F#实现Runge–Kutta算法求解微分方程。
842 0
F#实现Runge–Kutta算法求解常微分方程