秋分到了，深圳还是那么热。“秋天的第一杯奶茶”不知为何在朋友圈流传起来。
做为一只抽了三次工程硕士数学都没抽中的非酋菜鸡，不能和同龄人一起学习，只能默默自学。
废话不多说，今天先自学线性方程组的直接求解法。公式和推导相关参考书都有，这里主要关注怎么用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.列主元消去法
夜已深，明天继续........
今日结语：
虽然我知道我花这么多时间写这些对我的发展并没有什么用，即拿不到学分，也对之后的科研没有帮助。但是写代码时进行的逻辑思考的推理验证使我快乐。我并不打算做一名程序员，因为我知道写代码当成爱好是快乐，当成工作就只能是折磨了哈哈哈。
