高斯消元法求解线性方程组（附python代码）

输入：a是m×n的系数矩阵，b是m×1的（列）向量。

输出：方程组的通解。

用高斯消元法（行化简法）解线性方程组步骤

1.构造方程组的增广矩阵
2.从最左边列往右，使用行化简算法把增广矩阵化为阶梯形，确定矩阵是否有解：
若最后一列为主元列（最后一行非零行形如 [0 0 0 5]），无解，返回无解。
3.继续行化简，把主元上面的所有的元素都化为0，把主元位置变成1.
4.把每个主元列对应的变量表示成非主元变量的线性组合.

import numpy as np

def arguemented_mat(a, b):return np.c_[a,b]# 行交换
def swap_row(a,i,j):m,n = a.shapeif i >= m or j >= m:print 'error: out of index ...'else:a[i] = a[i] ^ a[j]a[j] = a[i] ^ a[j]a[i] = a[i] ^ a[j]return a# 变成阶梯矩阵
def trape_mat(sigma):m,n = sigma.shape#保存主元所在的列数，一般来说，每行都有一个主元，除非某行全零main_factor = []main_col = 0while main_col < n and len(main_factor) < m:# 当行数多于列数的时候，出现所有的列已经处理完，结束if main_col == n:break# 逐列找主元,若该列全零（从第i行往下），则没有主元# 当前查找小矩阵首行所在原矩阵的行号first_row = len(main_factor)while main_col < n:new_col = sigma[first_row:, main_col]not_zeros = np.where(new_col > 0)[0]# 若该列没有非零值，则该列非主元列if len(not_zeros) == 0:main_col += 1break# 否则通过行变换找到主元else:main_factor.append(main_col)index = not_zeros[0]# 若首个元素不是主元,需要行变换if index != 0:sigma = swap_row(sigma, first_row, first_row + index)# 把该主元下面的元素全部变成0if first_row < m - 1:for k in xrange(first_row+1, m):times = float(sigma[k, main_col]) / sigma[first_row, main_col]sigma[k] = sigma[k] - times * sigma[first_row]# 处理完当前主元列之后继续main_col += 1breakreturn sigma, main_factor# 回代求解
def back_solve(sigma, main_factor):# 判断是否有解if len(main_factor) == 0:print 'wrong main_factor ...'return Nonem,n = sigma.shapeif main_factor[-1] == n-1:print 'no answer ...'return None# 把所有的主元元素上方的元素变成0for i in range(len(main_factor)-1,-1,-1):factor = sigma[i, main_factor[i]]sigma[i] = sigma[i] / float(factor)for j in xrange(i):times = sigma[j, main_factor[i]]sigma[j] = sigma[j] - float(times)*sigma[i]# 先看看结果对不对return sigma# 结果打印
def print_result(sigma, main_factor):if sigma is None:print 'no answer ...'return m,n = sigma.shaperesult = [''] * (n-1)main_factor = list(main_factor)for i in range(n-1):# 如果不是主元列，则为自由变量if i not in main_factor:result[i] = 'x_' + str(i+1) + '(free var)'# 否则是主元变量，从对应的行，将主元变量表示成非主元变量的线性组合else:# row_of_maini表示该主元所在的行row_of_maini = main_factor.index(i)result[i] = str(sigma[row_of_maini, -1])for j in range(i+1, n-1):ratio = sigma[row_of_maini, j]if ratio > 0:result[i] = result[i] + '-' + str(ratio) + 'x_' + str(j+1)if ratio < 0:result[i] = result[i] + '+' + str(-ratio) + 'x_' + str(j+1)print '方程的通解是：\n', for i in range(n-1):print 'x_' + str(i+1), '=', result[i]            # 得到简化的阶梯矩阵和主元列
def solve(a,b):sigma = arguemented_mat(a,b)print '增广矩阵为：'print sigmasigma, main_factor = trape_mat(sigma)sigma = back_solve(sigma, main_factor)print '方程的简化阶梯矩阵:'print sigmaprint '方程的主元列为：'print main_factorprint_result(sigma, main_factor)return sigma, main_factorif __name__ == '__main__':a = np.array([[1,6,2,-5,-2], [0,0,2,-8,-1], [0,0,0,0,1]])b = np.array([-4, 3, 7])sigma,mf  = solve(a,b)print '*' * 20a = np.array([[1,0,-5], [0,1,1], [0,0,0]])b = np.array([1,4,0])sigma,mf  = solve(a,b)print '*' * 20

增广矩阵为：
[[ 1  6  2 -5 -2 -4][ 0  0  2 -8 -1  3][ 0  0  0  0  1  7]]
方程的简化阶梯矩阵:
[[ 1  6  0  3  0  0][ 0  0  1 -4  0  5][ 0  0  0  0  1  7]]
方程的主元列为：
[0, 2, 4]
方程的通解是：
x_1 = 0-6x_2-3x_4
x_2 = x_2(free var)
x_3 = 5+4x_4
x_4 = x_4(free var)
x_5 = 7
********************
增广矩阵为：
[[ 1  0 -5  1][ 0  1  1  4][ 0  0  0  0]]
方程的简化阶梯矩阵:
[[ 1  0 -5  1][ 0  1  1  4][ 0  0  0  0]]
方程的主元列为：
[0, 1]
方程的通解是：
x_1 = 1+5x_3
x_2 = 4-1x_3
x_3 = x_3(free var)
********************

高斯消元法求解线性方程组（附python代码）相关推荐

未知数数量大于方程数量，如何求解，附Python 代码
未知数如果大于方程数量,意味着限制少于自由度,方程要么无解,要么有无穷个解. 怎么办?我就要部分答案,只要满足方程限制就行? 示例如下: import sympy as sp x, y, z = sp ...
python重点知识归纳_一文了解机器学习知识点及其算法（附python代码）
一文了解机器学习知识点及其算法(附python代码) 来源:数据城堡时间:2016-09-09 14:05:50 作者: 机器学习发展到现在,已经形成较为完善的知识体系,同时大量的数据科学家的研究成 ...
教程 | 理解和实现自然语言处理终极指南（附Python代码）
教程 | 理解和实现自然语言处理终极指南(附Python代码) 时间 2017-02-16 14:41:39 机器之心原文 http://www.jiqizhixin.com/article ...
基于TextRank算法的文本摘要（附Python代码）
基于TextRank算法的文本摘要(附Python代码): https://www.jiqizhixin.com/articles/2018-12-28-18
python 合并工作簿_Excel:快速合并多张表格或多个文件（工作簿）的数据（附Python代码）...
Excel:快速合并多张表格或多个文件(工作簿)的数据(附Python代码) 现实工作中经常遇到将零散的原始数据合并统计的工作要求,如月度统计或年度统计等.原始数据的收集大多是按时间(如日期或小时)进 ...
XGBoost参数调优完全指南（附Python代码）
XGBoost参数调优完全指南(附Python代码) 译注:文内提供的代码和运行结果有一定差异,可以从这里下载完整代码对照参考.另外,我自己跟着教程做的时候,发现我的库无法解析字符串类型的特征,所以只 ...
python随机森林变量重要性_推荐：一文读懂随机森林的解释和实现（附python代码）...
原标题:推荐 :一文读懂随机森林的解释和实现(附python代码) 作者:WilliamKoehrsen:翻译:和中华:校对:李润嘉本文约6000字,建议阅读15分钟. 本文从单棵决策树讲起,然后逐 ...
材料力学求解器-刚架与桁架杆系的计算机求解（附matlab代码）
材料力学求解器-刚架与桁架杆系的计算机求解(附matlab代码) 1 刚架的计算机求解 1.1位移法与刚度矩阵 1.2 matlab程序 2 桁架的计算机求解材料力学是一门非常成熟的学科,里面有大量 ...
Excel:快速合并多张表格或多个文件（工作簿）的数据（附Python代码）
Excel:快速合并多张表格或多个文件(工作簿)的数据(附Python代码) 现实工作中经常遇到将零散的原始数据合并统计的工作要求,如月度统计或年度统计等.原始数据的收集大多是按时间(如日期或小时)进 ...
利用OpenSearch API检索和下载数据附Python代码实例
利用OpenSearch API检索和下载数据附Python代码实例在数据下载过程中,我们常常会需要下载非常多的数据文件,这时我们可以利用wget等软件或者编写数据下载脚本来实现数据下载的批处理. ...

高斯消元法求解线性方程组（附python代码）

用高斯消元法（行化简法）解线性方程组步骤

高斯消元法求解线性方程组（附python代码）相关推荐

最新文章

热门文章