高斯消去法
已知:
向前消去:
第i行的系数减去第一行的系数乘以,消去第i行的,重复操作,将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
}
}
向后回带
消去完毕后的矩阵方程可以通过回带的方式求解,如第N个方程只有一个未知量,可以直接计算得到:
可以通过的值计算得到:
以此类推:
伪代码如下:
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]