在 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