在 Python 中使用 Pivoting 进行高斯消除

Isaac Tony 2022年5月17日
在 Python 中使用 Pivoting 进行高斯消除

高斯消元法也称为行约简法。它是一种常用于解决线性问题的算法。该算法涉及对从线性方程中提取的系数矩阵进行一系列行操作,直到矩阵简化为梯形。

将矩阵转换为梯形形式时,会执行以下三个操作。这些包括:

  • 将一行与一个标量相乘,然后将其与另一行相加或减去。
  • 交换行。
  • 将行与标量相乘。

如果每行中的第一个元素(也称为前导条目)为零,则矩阵为行梯形;每行中的前导条目是前一行中前导条目右侧的一列。此外,具有零元素的行应低于具有非零行的元素。

行操作主要分为两种。前向消除涉及将矩阵简化为梯形,以确定是否存在可行的解决方案以及它是有限的还是无限的。另一方面,反向替换进一步将矩阵简化为行简化梯形。

在 Python 中使用枢轴进行高斯消除

枢轴是行和列的交换以获得合适的枢轴元素。与其他行条目相比,合适的枢轴元素应该既非零又显着大但更小。

旋转分为部分旋转和完全旋转。在部分枢轴算法下,最大元素被认为是枢轴元素以最小化舍入误差。

另一方面,完整的透视包括行和列的交换以获得最佳的透视元素,从而提高准确性。

如果没有变量的指数大于 1,则认为一组方程是线性的。高斯消元涉及一系列步骤;第一步包括创建一个系数矩阵。

系数矩阵只是从一组线性方程中得出的系数矩阵。下一步涉及创建一个增强矩阵,然后对其进行一系列操作,将其简化为梯形。

但是,在枢轴元素为零或非常小的情况下,我们必须将枢轴行与较低的行互换。

在实现高斯消元法时,我们需要注意数组索引。Python 可迭代对象(例如列表和数组)通常从索引 0 开始,到索引 n-1 结束。

然后我们可以读取增广矩阵的内容,应用消去法,反向代入,最后显示答案。

from numpy import array, zeros, fabs, linalg

a = array([[0, 6, -1, 2, 2],
           [0, 3, 4, 1, 7],
           [5, 1, 0, 3, -1],
           [3, 1, 3, 0, 2],
           [4, 4, 1, -2, 1]], float)
#the b matrix constant terms of the equations 
b = array([5, 7, 2, 3, 4], float)

print("Solution by NumPy:")


print(linalg.solve(a, b))

n = len(b)
x = zeros(n, float)

#first loop specifys the fixed row
for k in range(n-1):
    if fabs(a[k,k]) < 1.0e-12:
        
        for i in range(k+1, n):
            if fabs(a[i,k]) > fabs(a[k,k]):
                a[[k,i]] = a[[i,k]]
                b[[k,i]] = b[[i,k]]
                break

 #applies the elimination below the fixed row

    for i in range(k+1,n):
        if a[i,k] == 0:continue

        factor = a[k,k]/a[i,k]
        for j in range(k,n):
            a[i,j] = a[k,j] - a[i,j]*factor
            #we also calculate the b vector of each row
        b[i] = b[k] - b[i]*factor
print(a)
print(b)


x[n-1] = b[n-1] / a[n-1, n-1]
for i in range(n-2, -1, -1):
    sum_ax = 0
  
    for j in range(i+1, n):
        sum_ax += a[i,j] * x[j]
        
    x[i] = (b[i] - sum_ax) / a[i,i]

print("The solution of the system is:")
print(x)

输出:

Solution by NumPy:
[ 0.38947368  0.49473684 -0.10877193  0.12982456  0.83157895]
[[ 5.00000000e+00  1.00000000e+00  0.00000000e+00  3.00000000e+00  -1.00000000e+00]
 [ 0.00000000e+00  3.00000000e+00  4.00000000e+00  1.00000000e+00  7.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  4.50000000e+00  0.00000000e+00  6.00000000e+00]
 [ 0.00000000e+00  4.44089210e-16  0.00000000e+00  3.52702703e+00  2.95945946e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00  2.11354647e+00]]
[2.         7.         4.5        2.91891892 1.75758075]
The solution of the system is:
[ 0.38947368  0.49473684 -0.10877193  0.12982456  0.83157895]

在上面的代码中,我们使用了 NumPy 库提供的 array functionfabs 函数来创建一个矩阵并读取绝对值。我们也使用了 Linalg;一个 NumPy 子库,用于执行计算特征值、向量和行列式等操作。

from numpy import array, zeros, fabs, linalg

变量 a 表示变量的常数系数,而 b 表示对应于方程的常数项矩阵。通常,矩阵和向量是从我们打算求解的四个线性方程的矩阵表示中得出的。

a = array([[0, 6, -1, 2, 2],
           [0, 3, 4, 1, 7],
           [5, 1, 0, 3, -1],
           [3, 1, 3, 0, 2],
           [4, 4, 1, -2, 1]], float)
#the b matrix constant terms of the equations
b = array([5, 7, 2, 3, 4], float)

消除的目标是将矩阵变换为前导对角线为零。

执行高斯消除时,第一行应保持固定。但是,如果第一行中的枢轴元素等于 0,我们需要将第一行与任何后续行交换,而不是零作为前导元素。

我们在上面的代码中使用了 k 循环来指定固定的行。在这种情况下,k 表示固定行的索引,范围从 1n-1

随后的 if 语句检查枢轴元素是否为零并将其与适当的行交换。使用的 i 表示固定行下方的行的索引。

for k in range(n-1):
    if fabs(a[k,k]) < 1.0e-12:
       
        for i in range(k+1, n):
            if fabs(a[i,k]) > fabs(a[k,k]):
                a[[k,i]] = a[[i,k]]
                b[[k,i]] = b[[i,k]]
                break

在下面的代码中,我们使用循环来确保第一行不应用消除。

然后,我们将固定行下方的每一行乘以一个因子。因子等于对应于固定行的对角元素除以对应行中的第一个元素。

接下来,我们从固定行中减去两行,得到主对角线下第二列的元素。随后的行也可以在相同的过程中减少。

正如我们所提到的,第一行将始终保持不变。我们还计算了每行的因子和 b 向量,如下所示。

 for i in range(k+1,n):
       #first row should remain intact
        if a[i,k] == 0:continue
       #calculating the factor
        factor = a[k,k]/a[i,k]
        for j in range(k,n):
            a[i,j] = a[k,j] - a[i,j]*factor
            #calculating the b vector
        b[i] = b[k] - b[i]*factor
print(a)
print(b)

术语回代是指 x 的值可以从最后一个方程计算到第一个方程的事实。在这种情况下,我们利用 b 标量向量的相应值乘以 x 的系数。然后将结果除以主对角线中的相应元素。

我们优先计算从 n 开始的 x 值和从 n-1 到第一个的 x 的剩余值。另一方面,求和应该从 j+1 开始,因为从之前的计算中我们只知道 x 的值。此外,和的值是 x 的值乘以 a 的值的结果。我们现在可以计算 x 的值,因为我们有总和的值。

#calculating the value of x beginning from n-1
x[n-1] = b[n-1] / a[n-1, n-1]
for i in range(n-2, -1, -1):
    sum_ax = 0
#calculating the value of the sum
    for j in range(i+1, n):
        sum_ax += a[i,j] * x[j]
#calculating the value of x
    x[i] = (b[i] - sum_ax) / a[i,i]

print(x)

输出:

[[ 2.00000000e+00  3.00000000e+00  4.00000000e+00  1.00000000e+00
   7.00000000e+00]
 [ 0.00000000e+00  7.00000000e+00 -1.00000000e+00  3.00000000e+00
   1.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00 -1.30000000e+01  2.00000000e+00
  -2.10000000e+01]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  2.81250000e+00
   5.81250000e+00]
 [ 0.00000000e+00  8.88178420e-16  0.00000000e+00 -4.44089210e-16
   4.95734797e+00]]
[  7.           5.         -14.           0.625        0.15371622]
The solution of the system is:
[0.02170543 0.79224806 1.05116279 0.15813953 0.03100775]

下面的 NumPy 函数也应该得到相同的答案。

print("The numpy solution:")
print(linalg.solve(a, b))

输出:

Solution by Nuumpy:
[0.02170543 0.79224806 1.05116279 0.15813953 0.03100775]
Author: Isaac Tony
Isaac Tony avatar Isaac Tony avatar

Isaac Tony is a professional software developer and technical writer fascinated by Tech and productivity. He helps large technical organizations communicate their message clearly through writing.

LinkedIn