高斯消去法(python实现)

高斯消去法

已知：$A\varphi =b$

向前消去：

第i行的系数减去第一行的系数乘以$\frac{a_{i1}}{a_{11}}$,消去第i行的${\varphi }_1$,重复操作，将A处理为一个上三角矩阵。伪代码如下：

for k = 1 to N-1
{
		for i = k + 1 to N  //从第二行开始计算
    {
      ratio = a_ik/a_kk
      {
        for j = k to N   // 从第k列开始
          	a_ij = a_ij - ratio * a_kj
      }
      b_i = b_i - ratio * b_k
    }
}

向后回带

消去完毕后的矩阵方程可以通过回带的方式求解${\varphi }_i$,如第N个方程只有一个未知量${\varphi }_N$,可以直接计算得到：
$${\varphi }_N=\frac{b_N^{N-1}}{a_{NN}^{N-1}}$$
${\varphi }_{N-1}$可以通过${\varphi }_N$的值计算得到：
$${\varphi }_{N-1}=\frac{b_{N-1}^{N-2}-a_{(N-1)N}^{N-2}{\varphi }_N}{a_{(N-1)(N-1)}^{N-2}}$$
以此类推：
$${\varphi }_i=\frac{b_i^{i-1}-{\sum }_{j=i+1}^N{a}_{ij}^{i-1}{\varphi }_j}{a_{ii}^{i-1}}$$
伪代码如下：

phi_N = b_N / a_NN
for i = N-1 to 1
{
		term = 0
		for j = i+1 to N
		{
				term = term + a_ij * phi_j
		}
		phi_i = (b_i - term) / a_ii
}

最终的python实现：

import numpy as np
if __name__ == '__main__':
    A = [[2, -1, 3, 2],
         [3, -3, 3, 2],
         [3, -1, -1, 2],
         [3, -1, 3, -1]]
    B = [6, 5, 3, 4]
    A = np.array(A).astype(float)
    B = np.array(B).astype(float)

    for k in range(A.shape[0]-1): # 0-n-1行
        for i in range(k+1, A.shape[0]): # k+1 - n 行
            ratio = A[i][k] / A[k][k]
            for j in range(k, A.shape[1]): # k - n 列
                A[i][j] = A[i][j] - ratio * A[k][j]
            B[i] = B[i] - ratio * B[k]
    print(A)
    print(B.size)
    n = B.size-1
    phi = np.zeros(B.size)
    phi[n] = B[n] / A[n][n]
    for i in range(n-1, -1, -1):
        term = 0
        for j in range(i+1, n+1): # i+1 - B.size
            term = term + A[i][j] * phi[j]
        phi[i] = (B[i]-term) / A[i][i]
    print(phi) # [1 1 1 1]