python18：用python实现高斯消元法

"""
高斯消元法:
分为三部分去解决这个大问题，并且先用例子，再推广为一般式，再转化为python code
1. 高斯消元法在化学配平上的体现，如何转化为matrix形式
2. 考虑(Upper) triangular systems的解决办法
3. general化去解决配平问题
4. 整合（为什么要移位的原因在图片3）

Part1.

some chemistry：
Eg.
Unbalanced equation:
Fe2O3 + C -> Fe + CO2

count atoms:
ATOM Reactants products
Fe      2       1
O       3       2
C       1       1

ATOM Reactants products
Fe      4       4
O       6       6
C       3       3

so:
Balanced equation:
2Fe2O3 + 3C -> 2Fe + 3CO2

Let's balance equation algorithmically:
Eg.
General equation:
aFe2O3 + bC -> cFe + dCO2
so:
2a = c
3a = 2d
b = d

Reached a very general problem:
solve linear systems:
input: set of n linear equations in n variables
output: assignment of values to variables satisfying all equations

转化为matrix：
2a = c          2a - c = 0          2a + 0b - 1c + 0d = 0
3a = 2d         3a - 2d = 0         3a + 0b + 0c - 2d = 0
b = d           b - d = 0           0a + 1b + 0c - 1d = 0
a = 4           a = 4               1a + 0b + 0c + 0d = 4
so：转化为matrix：
[        *       [   =  [
2 0 -1 0        a       0
3 0 0 -2        b       0
0 1 0 -1        c       0
1 0 0 0         d       4
]               ]       ]
coefficient matrix A
            solution vector x
                        right hand side b
solve Linear systems:
input: n*n-matrixA(python table) of coefficients, n -dim vector b(python list)
output: n-dim vector x(python list) such that Ax = b

What are easy to solve linear systems?
Overall strategy: Transform and Conquer Paradigm
problem input ->(Transform) another problem/another representation/another input ->(Conquer) solution

Part2.solve (Upper) triangular systems

Eg.how to solve (Upper) triangular systems(because only coefficients in main diagonal and above are non-zero)
1 2 1   a   2
0 1 3   b   7
0 0 2   c   6

Back-substitution algorithm
2c = 6(c = 3)
b + 3c = 7(c is known)(b=-2)
a + 2a + 2 = 2(b,c is known)(a=3)

How to implement in python:
def solve_by_back_sub(u, b):
    n = len(b)
    x = n * [0]
    for i in range(n-1, -1, -1):
        # find value x[i]
        # 此处看图片的详细解释
        s = sum(x[j] * u[i][j] for j in range(n-1, i, -1))
        x[i] = (b[i] - s)/u[i][i]
    return x

Part3.消元

Gaussian Elimination：
general system ->(Transform) triangular system ->(Conquer) solution
所以要做的就是消元：
simplify column

详细解释看图
def triangular(a, b):
    n = len(a)
    for j in range（n）:
        for i in range（j + 1， n）
            # eliminate entry a[i][j]
            q = a[i][j] / a[j][j]
            # subtract q*a[j] from a[i]
            a[i] = [a[i][l] - q*a[j][l] for l in range(n)]
            # subtract q*b[j] from b[j]
            b[i] = b[i] - q*b[j]
    return a, b

Part4.整合

最后，整合成完整的代码并分析其复杂度为O（n^3）:
def solve_by_back_sub(u, b):
    n = len(b)
    x = n * [0]
    for i in range(n-1, -1, -1):
        # find value x[i]
        # 此处看图片的详细解释
        s = sum(x[j] * u[i][j] for j in range(n-1, i, -1))
        x[i] = (b[i] - s)/u[i][i]
    return x

def triangular(a, b):
    #用deepcopy不破坏input
    u, c = deepcopy(a), deepcopy(b) n = len(u)
    for j in range(n):----------O(n)
        k = pivot_index(u, j)
        # 移位，原因在图片
        u[j], u[k] = u[k], u[j]
        c[j], c[k] = c[k], c[j]
        for i in range(j + 1, len(a)):----------O(n^2)
            q = u[i][j] / u[j][j]
            u[i] = [u[i][l] - q*u[j][l] for l in range(n)] ----------O(n^3)
            c[i] = c[i] – q*c[j]
        #INV: u, c has same solution(s) as a, b
    return u, c

def solve_gauss_elim(a, b):
    u, c = triangular(a, b)
    return solve_backsub(u, c)
"""