【python3】【计算方法】高斯消元法求解线性方程组 [a]{x} = b

高斯消元法求解线性方程组 [a]{x} = b

求解过程

使用高斯消元法求解线性方程组 [a]{x} = b，其中 [a] 是一个 Vandermonde 矩阵，{x} 是未知向量，b 是常向量。

高斯消元法的求解过程可以表示为以下矩阵变换：

$$\left[\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,n}& b_1\\ a_{2,1}& a_{2,2}& \cdots & a_{2,n}& b_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{n,1}& a_{n,2}& \cdots & a_{n,n}& b_n\end{array}\right]⟶\left[\begin{array}{ccccc}1& a_{1,2}^{\prime }& \cdots & a_{1,n}^{\prime }& b_1^{\prime }\\ 0& 1& \cdots & a_{2,n}^{\prime }& b_2^{\prime }\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 1& b_n^{\prime }\end{array}\right]$$

其中， $a_{i,j}^{\prime }$ 表示经过消元得到的新系数， $b_i^{\prime }$ 表示经过消元得到的新常数项。
无序列表形式：

• 先将第一行除以 $a_{1,1}$。
• 将第一行的倍数加到后面的行上，使得第一列下面的元素都为零。
• 用同样的方法消去第二行、第三行，直到最后一行。
• 经过消元得到上三角矩阵后，可以通过回带求解出未知向量 {x} 的值。
• 从最后一行开始，依次求解每个未知量，直到求解出第一个未知量 {x}_1。

回带的过程可以表示为以下公式：

$$x_n=\frac{b_n^{\prime }}{a_{n,n}^{\prime }}$$

$$x_i=\frac{1}{a_{i,i}^{\prime }}\left(b_i^{\prime }-\sum _{j=i+1}^n{a}_{i,j}^{\prime }{x}_j\right)\text{for }i=n-1,n-2,\dots ,1$$

其中， $n$ 是未知向量的长度。

案例

代码

from numpy import zeros, float64, array, product, \
                  diagonal, matmul
from gaussElimin import *

def vandermode(v):
    n = len(v)
    a = zeros((n,n), dtype=float64)
    for j in range(n):
        a[:,j] = v**(n-j-1)
    return a

v = array([1.0, 1.2, 1.4, 1.6, 1.8, 2.0])
b = array([0.0, 1.0, 0.0, 1.0, 0.0, 1.0])
a = vandermode(v)
aOrig = a.copy()    # Save original matrix
bOrig = b.copy()    # and the constant vector
x = gaussElimin(a,b)
det = product(diagonal(a))
print('x =\n',x)
print('\ndet =',det)
print('\nCheck result: [a]{x} - b =\n', \
      matmul(aOrig,x) - bOrig)
input("\nPress return to exit")

解析

• 代码定义了一个名为 vandermode 的函数。
• 该函数接受一个一维 Numpy 数组 v 并返回一个 Vandermonde 矩阵 a，其元素是 v 中元素的幂。
• 代码使用此函数生成一个 Vandermonde 矩阵 a 和一个常向量 b。
• 它使用高斯消元法解出线性方程组 [a]{x} = b，其中 {x} 是解向量。
• 矩阵 a 的行列式也被计算了。
• 然后代码打印出解向量 x 和行列式 det。
• 最后，代码通过从原始矩阵 aOrig 和向量 bOrig 的乘积中减去 b 来检查结果，其中 aOrig 和 bOrig 是原始矩阵 a 和向量 b 的副本。
• 代码末尾的 input 语句等待用户按回车键退出程序。

算法

## module gaussElimin
''' x = gaussElimin(a,b).
    Solves [a]{b} = {x} by Gauss elimination.
'''
import numpy as np

def gaussElimin(a,b):
    n = len(b)
  # Elimination Phase
    for k in range(0,n-1):
        for i in range(k+1,n):
           if a[i,k] != 0.0:
               lam = a [i,k]/a[k,k]
               a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]
               b[i] = b[i] - lam*b[k]
  # Back substitution
    for k in range(n-1,-1,-1):
        b[k] = (b[k] - np.dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
    return b

• 这是一个名为 gaussElimin 的模块，其中包含一个函数 gaussElimin(a,b)，用于使用高斯消元法求解线性方程组 [a]{b} = {x}。

• 此函数接受两个参数 a 和 b，其中 a 是一个 Numpy 数组，表示系数矩阵，b 是一个 Numpy 数组，表示常数向量。

• 该函数的求解过程分为两个阶段：消元阶段和回带阶段。

• 在消元阶段，函数使用高斯消元法将系数矩阵 a 转化为上三角矩阵。

• 在回带阶段，函数使用回带法求解未知向量 {x}。

• 具体来说，消元阶段的过程是：对于每一行，找到第一个非零元素，将该元素所在的列消成上三角形式。

• 对于第 $k$ 行和第 $i$ 行（$i>k$），如果 $a_{i,k}\ne 0$，则计算消元因子 $\lambda =a_{i,k}/a_{k,k}$，然后对第 $i$ 行做如下操作：

  $$a_{i,k+1:n}=a_{i,k+1:n}-\lambda a_{k,k+1:n}$$

  同时，将常数向量 $b$ 的第 $i$ 个元素也更新：

  $$b_i=b_i-\lambda b_k$$

• 回带阶段的过程是：从最后一行开始，依次求解每个未知量，直到求解出第一个未知量 $x_1$。

• 对于第 $k$ 行，设 $x_{k+1:n}$ 表示未知向量的后 $n-k$ 个元素，使用如下公式求解 $x_k$：

  $$x_k=\frac{1}{a_{k,k}}\left(b_k-\sum _{j=k+1}^n{a}_{k,j}{x}_j\right)$$

• 最终，函数返回解向量 {x}。

案例计算结果

x =
 [   416.66666667  -3125.00000004   9250.00000012 -13500.00000017
   9709.33333345  -2751.00000003]

det = -1.1324620799859046e-06

Check result: [a]{x} - b =
 [ 0.00000000e+00  0.00000000e+00  1.81898940e-12  1.09139364e-11
 -1.81898940e-11  4.36557457e-11]

Press return to exit


Process finished with exit code 0

附：思维导图