秋分到了，深圳还是那么热。“秋天的第一杯奶茶”不知为何在朋友圈流传起来。

做为一只抽了三次工程硕士数学都没抽中的非酋菜鸡，不能和同龄人一起学习，只能默默自学。

废话不多说，今天先自学线性方程组的直接求解法。公式和推导相关参考书都有，这里主要关注怎么用Python实现这些数值计算方法。能力有限，如有误请反馈于我及时修改。另外，以下代码没有考虑是否达到时间复杂度最优，只是可以运行得到理想的结果。

一、线性方程组的直接求解法

1.高斯消元法

算法实现思路：

1.第一步，将增广矩阵消元形成上三角矩阵；

2.第二步，从下向上回带进行解方程。

参考《数值分析》(李庆扬等)的第145页的推导公式。对代码的解释放在了代码的注释当中。只使用Numpy工具包。

以下是实现消元成上三角矩阵的代码。

def getTriangularMatrix(matA):

'''

matA:增广矩阵；array格式，且数值类型必须为浮点数，否则可能出现失真。

该函数为输入一个增广矩阵，输出为一个消元后的上三角矩阵

'''

if matA[0,0]==0: # 第一个元素为0的情况

print('不符合要求！需要换行操作。')

else:

numRows = matA.shape[0] # 获得总行数

for row in np.arange(0,numRows-1): # 前n-1行，row代表第几行

for k in range(row+1,numRows): # 在第row行的情况下，遍历剩下的行

if matA[k,row] != 0:

#第k(k小于row)行的新值等于该行减去第row行的值乘于一个系数。这个系数存在目的就是消元

matA[k,:]=matA[k,:]-(matA[k,row]/matA[row,row])*matA[row,:]

return matA # 此时输入的矩阵也被改变

以下是回代求解的代码

# 回代求解

def Gauss_solve(matA):

numRows = matA.shape[0] # 获得总行数

X = np.zeros((numRows,1)) # 建一个n*1的0列向量，来存结果

A = matA[0:numRows,0:numRows] # 为了方便表示，取出A

b = matA[0:numRows,-1] # 为了方便表示，取出b

if A[-1,-1] != 0: # 简单判断一下

X[numRows-1,0]=b[-1]/A[-1,-1] # 这里先求最后的那个解，目前我没想到不先求这个解的代码写法

for row_ in np.arange(numRows-2,-1,-1): # 着手求其它解

if A[row_,row_] != 0: # 简单判断一下

temp = np.dot(A[row_,row_+1:numRows],X[row_+1:numRows]) # 这一步花了很多时间去想。其思想在后面的图中。

X[row_,0] = (b[row_]-np.sum(temp))/A[row_,row_] # 有了上一行代码的基础，这一行就是简单的推导公式

return X

来两个例子测试一下效果。首先，先生成两个不同的测试矩阵。

import numpy as np

test_mat1 = np.array([[2,0,6,20],[1,4,3,18],[3,2,1,10]],dtype=float)

test_mat1 = array([[ 2., 0., 6., 20.],

[ 1., 4., 3., 18.],

[ 3., 2., 1., 10.]])

test_mat2 = np.array([[2,1,0,3,16],[0,2,0,4,20],[5,0,1,2,16],[1,1,2,2,17]],dtype=float)

test_mat2 = array([[ 2., 1., 0., 3., 16.],

[ 0., 2., 0., 4., 20.],

[ 5., 0., 1., 2., 16.],

[ 1., 1., 2., 2., 17.]])

开始运行函数计算结果，代码如下：

X1 = Gauss_solve(getTriangularMatrix(test_mat1))

X1 = array([[1.],

[2.],

[3.]])

X2 = Gauss_solve(getTriangularMatrix(test_mat2))

X2 = array([[1.],

[2.],

[3.],

[4.]])

得到两组结果，你可以代入验算，结果是正确的。暂时没有再做更多算例测试。

这里再展示一个Python的scipy库中直接求解线性方程组的函数。

from scipy import linalg

A = test_mat2[0:4,0:4] # 取出第二个算例的A

b = test_mat2[0:4,-1] # 取出第二个算例的b

x = inalg.solve(A,b)

# 解得

x = array([1., 2., 3., 4.])

2.列主元消去法

夜已深，明天继续........

今日结语：

虽然我知道我花这么多时间写这些对我的发展并没有什么用，即拿不到学分，也对之后的科研没有帮助。但是写代码时进行的逻辑思考的推理验证使我快乐。我并不打算做一名程序员，因为我知道写代码当成爱好是快乐，当成工作就只能是折磨了哈哈哈。