在 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