在 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 function 和 fabs 函式來建立一個矩陣並讀取絕對值。我們也使用了 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 表示固定行的索引，範圍從 1 到 n-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]
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