Python でのピボットを使用したガウスの消去法
ガウスの消去法は、行削減法としても知られています。これは、線形問題を解決するために一般的に使用されるアルゴリズムです。アルゴリズムには、行列が階段形に縮小されるまで、線形方程式から引き出された係数の行列に対する一連の行演算が含まれます。
行列を階段形に変換するときに、次の 3つの操作が実行されます。これらには以下が含まれます：
- 1つの行にスカラーを乗算し、それを加算するか、別の行から減算します。
- 行を入れ替えます。
- 行にスカラーを掛けます。
各行の最初の要素（先頭のエントリとも呼ばれる）がゼロの場合、行列は行階段形になります。各行の先頭のエントリは、前の行の先頭のエントリの右側の 1 列です。さらに、要素がゼロの行は、行がゼロ以外の要素の下にある必要があります。
行演算は大きく 2つに分類されます。前方排除では、行列を階段形に縮小して、実行可能な解が存在するかどうか、およびそれが有限か無限かを判断します。一方、逆置換は、行列を行階段形にさらに縮小します。
Python でのピボットによるガウスの消去法
ピボットとは、適切なピボット要素を取得するための行と列の交換です。適切なピボット要素は、ゼロ以外で、他の行エントリと比較するとかなり大きいが小さい必要があります。
ピボットは、部分ピボットと完全ピボットに分類されます。部分ピボットアルゴリズムでは、丸め誤差を最小限に抑えるために、最大の要素がピボット要素と見なされます。
一方、完全なピボットには、行と列を交換して最適なピボット要素を取得することが含まれるため、精度が向上します。
指数が複数の変数がない場合、一連の方程式は線形と見なされます。ガウスの消去法には、一連の手順が含まれます。最初のステップは、係数行列を作成することです。
係数行列は、一連の線形方程式から引き出された係数の行列です。次のステップでは、拡大行列を作成します。この拡大行列は、階段形に縮小する一連の操作にかけられます。
ただし、ピボット要素がゼロまたは非常に小さい場合は、ピボット行を下の行と交換する必要があります。
ガウスの消去法を実装する際、配列インデックスに注意する必要があります。リストや配列などの 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 ライブラリが提供する配列関数と fabs 関数を使用して、行列を作成し、絶対値を読み取りました。Linalg も使用しました。固有値やベクトル、行列式の計算などの操作を実行するために使用される NumPy サブライブラリ。
from numpy import array, zeros, fabs, linalg
変数 a は変数の定数係数を表し、b は方程式に対応する定数項の行列を表します。一般に、行列とベクトルは、解く予定の 4つの線形方程式の行列表現から描画されます。
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)
除去の目的は、先頭の対角にゼロを持つように行列を変換することです。
ガウスの消去法を実行する場合、最初の行は固定されたままにする必要があります。ただし、最初の行のピボット要素がゼロに等しい場合は、最初の行を、先頭の要素としてゼロがない後続の行と交換する必要があります。
上記のコードでは、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
以下のコードでは、ループを使用して、最初の行に消去が適用されないようにします。
次に、固定行の下の各行に係数を掛けました。係数は、固定行に対応する対角要素を対応する行の最初の要素で割ったものに等しくなります。
次に、固定行から 2つの行を減算して、主対角線の下の 2 番目の列の要素を取得します。後続の行も同じ手順で減らすことができます。
前述したように、最初の行は常に変更されません。また、以下に示すように、各行の係数と 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 の値が最後の方程式から最初の方程式まで計算できるという事実を指します。この場合、x の係数を掛けた b スカラーベクトルの対応する値を利用します。次に、結果は主対角線の対応する要素で除算されます。
n から始まる x の値の計算と、n-1 から最初の値までの x の残りの値の計算を優先しました。一方、前の計算では x の値しかわからないため、合計は j +1 から開始する必要があります。さらに、合計の値は、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