高斯消去法(python实现)

303 阅读1分钟

高斯消去法

在这里插入图片描述

已知:Aϕ=bA\phi=b

向前消去:

第i行的系数减去第一行的系数乘以ai1a11\frac{a_{i1}}{a_{11}},消去第i行的ϕ1\phi_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
    }
}

在这里插入图片描述

向后回带

消去完毕后的矩阵方程可以通过回带的方式求解ϕi\phi_i,如第N个方程只有一个未知量ϕN\phi_N,可以直接计算得到:

ϕN=bNN1aNNN1\phi_N=\frac{b_N^{N-1}}{a_{NN}^{N-1}}

ϕN1\phi_{N-1}可以通过ϕN\phi_N的值计算得到:

ϕN1=bN1N2a(N1)NN2ϕNa(N1)(N1)N2\phi_{N-1}=\frac{b_{N-1}^{N-2}-a_{(N-1)N}^{N-2}\phi_N}{a_{(N-1)(N-1)}^{N-2}}

以此类推:

ϕi=bii1j=i+1Naiji1ϕjaiii1\phi_i = \frac{b_i^{i-1}-\sum_{j=i+1}^Na_{ij}^{i-1}\phi_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]