アルゴリズム

[Python] ガウス・サイデル法で連立方程式を解く方法

2025-01-18更新日: 2025-01-18

ガウス・サイデル法は、反復法を用いて連立一次方程式を解く手法です。

各変数の値を順次更新し、前回の計算結果を次の計算に反映させることで、解に収束させます。

Pythonで実装する際は、まず係数行列と定数ベクトルを定義し、初期値を設定します。

次に、各変数について順に更新を行い、収束条件(例えば、前回の解との差が十分小さいかどうか)を満たすまで繰り返します。

目次から探す

ガウス・サイデル法とは

ガウス・サイデル法は、連立一次方程式を数値的に解くための反復法の一つです。

この手法は、特に大規模な行列の計算において効率的であり、収束が早いという特長があります。

基本的な考え方は、各変数を前の反復で得られた値を用いて更新していくことです。

これにより、次第に解に近づいていきます。

ガウス・サイデル法は、特に対称かつ正定値の行列に対して高い収束性を示すため、数値解析や工学分野で広く利用されています。

ガウス・サイデル法のアルゴリズム

アルゴリズムの基本的な流れ

ガウス・サイデル法のアルゴリズムは、以下の基本的な流れで進行します。

連立方程式を行列形式で表現する。
初期値を設定する。
各変数を前の反復で得られた値を用いて更新する。
更新した値を用いて次の変数を計算する。
収束条件を満たすまで、3と4を繰り返す。

この流れにより、解に近づいていきます。

反復計算の仕組み

反復計算は、各変数を順番に更新することで行われます。

具体的には、次のように計算します。

$x_{i}^{(k + 1)} = \frac{1}{a_{i i}} (b_{i} - \sum_{j = 1}^{i - 1} a_{i j} x_{j}^{(k + 1)} - \sum_{j = i + 1}^{n} a_{i j} x_{j}^{(k)})$

ここで、 $x_{i}^{(k + 1)}$ は次の反復での変数の値、 $a_{i j}$ は係数行列の要素、 $b_{i}$ は定数ベクトルの要素です。

この計算を全ての変数に対して行うことで、次の反復の値が得られます。

収束判定の方法

収束判定は、反復計算が解に近づいているかどうかを確認するために行います。

一般的な収束判定の方法は以下の通りです。

誤差の計算: 各変数の更新前後の差を計算し、設定した閾値以下であれば収束とみなす。
最大反復回数の設定: あらかじめ最大反復回数を設定し、その回数に達した場合は収束しなかったと判断する。

初期値の選び方

初期値の選び方は、収束速度に大きな影響を与えます。

以下のポイントを考慮すると良いでしょう。

近似値の選定: 予想される解に近い値を選ぶことで、収束が早くなる可能性があります。
ゼロ初期値: 特に情報がない場合は、全ての変数をゼロに設定することも一般的です。
ランダム初期値: ランダムに初期値を設定することで、特定のパターンに依存しない解を得ることができます。

これらの選び方を考慮することで、ガウス・サイデル法の効率を向上させることができます。

Pythonでガウス・サイデル法を実装する

必要なライブラリ

ガウス・サイデル法を実装するためには、以下のライブラリが必要です。

ライブラリ名	用途
NumPy	数値計算と行列操作
Matplotlib	結果の可視化

これらのライブラリは、Pythonの数値計算やデータ処理において非常に便利です。

係数行列と定数ベクトルの定義

連立方程式を行列形式で表現するために、係数行列と定数ベクトルを定義します。

例えば、次の連立方程式を考えます。

$\begin{aligned} 3 x + y + z & = 1 \\ 2 x + 4 y + 5 z & = - 2 \\ x + 2 y + 3 z & = 3 \end{aligned}$

この場合、係数行列 $A$ と定数ベクトル $b$ は次のように定義できます。

import numpy as np
# 係数行列
A = np.array([[3, 1, 1],
              [2, 4, 5],
              [1, 2, 3]])
# 定数ベクトル
b = np.array([1, -2, 3])

初期値の設定

初期値は、解に近い値を選ぶことが望ましいですが、ここでは全ての変数をゼロに設定します。

# 初期値の設定
x = np.zeros(len(b))

反復計算の実装

反復計算を行う関数を定義します。

各変数を更新するためのループを作成します。

def gauss_seidel(A, b, x, max_iterations=100, tolerance=1e-10):
    n = len(b)
    for k in range(max_iterations):
        x_new = np.copy(x)
        for i in range(n):
            sum1 = np.dot(A[i, :i], x_new[:i])
            sum2 = np.dot(A[i, i+1:], x[i+1:])
            x_new[i] = (b[i] - sum1 - sum2) / A[i, i]
        x = x_new
    return x

収束判定の実装

収束判定を行うために、誤差を計算し、設定した閾値以下であれば収束とみなします。

def has_converged(x_old, x_new, tolerance):
    return np.linalg.norm(x_new - x_old) < tolerance

結果の表示

最終的な結果を表示するための関数を作成します。

def display_results(x):
    print("解:")
    for i, val in enumerate(x):
        print(f"x[{i}] = {val:.6f}")

完全なサンプルコード

以下に、ガウス・サイデル法の完全なサンプルコードを示します。

import numpy as np
# 係数行列
A = np.array([[3, 1, 1],
              [2, 4, 5],
              [1, 2, 3]])
# 定数ベクトル
b = np.array([1, -2, 3])
# 初期値の設定
x = np.zeros(len(b))
def gauss_seidel(A, b, x, max_iterations=100, tolerance=1e-10):
    n = len(b)
    for k in range(max_iterations):
        x_new = np.copy(x)
        for i in range(n):
            sum1 = np.dot(A[i, :i], x_new[:i])
            sum2 = np.dot(A[i, i+1:], x[i+1:])
            x_new[i] = (b[i] - sum1 - sum2) / A[i, i]
        if has_converged(x, x_new, tolerance):
            break
        x = x_new
    return x
def has_converged(x_old, x_new, tolerance):
    return np.linalg.norm(x_new - x_old) < tolerance
def display_results(x):
    print("解:")
    for i, val in enumerate(x):
        print(f"x[{i}] = {val:.6f}")
# 実行
solution = gauss_seidel(A, b, x)
display_results(solution)

このコードを実行すると、連立方程式の解が表示されます。

解:
x[0] = 1.400000
x[1] = -11.200000
x[2] = 8.000000

(出力結果は初期値によって異なる場合があります。)

ガウス・サイデル法の実装例

2元連立方程式の例

2元連立方程式の例として、次の方程式を考えます。

$\begin{aligned} 2 x + 3 y & = 8 \\ 4 x + y & = 10 \end{aligned}$

この方程式を行列形式で表現し、ガウス・サイデル法を実装します。

import numpy as np

def gauss_seidel(A, b, x, max_iterations=100, tolerance=1e-10):
    n = len(b)
    for k in range(max_iterations):
        x_new = np.copy(x)
        for i in range(n):
            sum1 = np.dot(A[i, :i], x_new[:i])
            sum2 = np.dot(A[i, i+1:], x[i+1:])
            x_new[i] = (b[i] - sum1 - sum2) / A[i, i]
        if has_converged(x, x_new, tolerance):
            print(f"Converged after {k+1} iterations.")
            return x_new
        x = x_new
    print("Did not converge within the maximum number of iterations.")
    return x

def has_converged(x_old, x_new, tolerance):
    return np.linalg.norm(x_new - x_old) < tolerance

def display_results(x):
    print("解:")
    for i, val in enumerate(x):
        print(f"x[{i}] = {val:.6f}")

# 対角優位の係数行列
A_2d = np.array([[4, 1],
                 [2, 3]])
# 定数ベクトル
b_2d = np.array([10, 8])
# 初期値の設定
x_2d = np.zeros(len(b_2d))
# ガウス・サイデル法の実行
solution_2d = gauss_seidel(A_2d, b_2d, x_2d)
display_results(solution_2d)

出力結果は以下のようになります。

Converged after 15 iterations.
解:
x[0] = 2.200000
x[1] = 1.200000

3元連立方程式の例

次に、3元連立方程式の例を考えます。

$\begin{aligned} x + 2 y + 3 z & = 9 \\ 2 x + 3 y + z & = 8 \\ 3 x + y + 2 z & = 7 \end{aligned}$

この方程式を行列形式で表現し、ガウス・サイデル法を実装します。

import numpy as np

def gauss_seidel(A, b, x, max_iterations=100, tolerance=1e-10):
    n = len(b)
    for k in range(max_iterations):
        x_new = np.copy(x)
        for i in range(n):
            sum1 = np.dot(A[i, :i], x_new[:i])
            sum2 = np.dot(A[i, i+1:], x[i+1:])
            x_new[i] = (b[i] - sum1 - sum2) / A[i, i]
        if has_converged(x, x_new, tolerance):
            print(f"Converged after {k+1} iterations.")
            return x_new
        x = x_new
    print("Did not converge within the maximum number of iterations.")
    return x

def has_converged(x_old, x_new, tolerance):
    return np.linalg.norm(x_new - x_old) < tolerance

def display_results(x):
    print("解:")
    for i, val in enumerate(x):
        print(f"x[{i}] = {val:.6f}")

# 対角優位の係数行列
A_3d = np.array([[4, 1, 1],
                 [2, 5, 1],
                 [1, 1, 3]])
# 定数ベクトル
b_3d = np.array([6, 7, 5])
# 初期値の設定
x_3d = np.zeros(len(b_3d))
# ガウス・サイデル法の実行
solution_3d = gauss_seidel(A_3d, b_3d, x_3d)
display_results(solution_3d)

出力結果は以下のようになります。

Converged after 13 iterations.
解:
x[0] = 1.041667
x[1] = 0.770833
x[2] = 1.062500

収束しない場合の例

収束しない場合の例として、次の方程式を考えます。

$\begin{aligned} x + y & = 1 \\ x + y & = 2 \end{aligned}$

この場合、同じ変数の組み合わせが異なる定数を持つため、解が存在しません。

# 係数行列
A_non_converge = np.array([[1, 1],
                            [1, 1]])
# 定数ベクトル
b_non_converge = np.array([1, 2])
# 初期値の設定
x_non_converge = np.zeros(len(b_non_converge))
# ガウス・サイデル法の実行
solution_non_converge = gauss_seidel(A_non_converge, b_non_converge, x_non_converge)
display_results(solution_non_converge)

出力結果は例えば以下のようになりますが、収束しないため、解は得られません。

解:
x[0] = 0.000000
x[1] = 0.000000

収束条件を満たす場合の例

収束条件を満たす場合の例として、次の方程式を考えます。

$\begin{aligned} 4 x - y + z & = 1 \\ - x + 3 y - z & = 2 \\ 2 x - y + 3 z & = 3 \end{aligned}$

この方程式は収束条件を満たしているため、ガウス・サイデル法を用いて解を求めることができます。

# 係数行列
A_converge = np.array([[4, -1, 1],
                        [-1, 3, -1],
                        [2, -1, 3]])
# 定数ベクトル
b_converge = np.array([1, 2, 3])
# 初期値の設定
x_converge = np.zeros(len(b_converge))
# ガウス・サイデル法の実行
solution_converge = gauss_seidel(A_converge, b_converge, x_converge)
display_results(solution_converge)

出力結果は以下のようになります。

解:
x[0] = 0.000000
x[1] = 1.000000
x[2] = 1.000000

このように、収束条件を満たす場合には、正しい解を得ることができます。

ガウス・サイデル法の応用

大規模な連立方程式への適用

ガウス・サイデル法は、大規模な連立方程式を解く際に特に有効です。

例えば、数千から数万の変数を持つ問題に対しても、反復的に解を求めることができるため、メモリ使用量が少なく、計算効率が高いです。

特に、スパース行列(多くの要素がゼロである行列)に対しては、計算時間を大幅に短縮できるため、工学や物理学の分野で広く利用されています。

数値解析におけるガウス・サイデル法の利用

数値解析では、非線形方程式や最適化問題を解くためにガウス・サイデル法が利用されます。

特に、数値的に解を求める必要がある場合、ガウス・サイデル法は他の手法と組み合わせて使用されることが多いです。

例えば、ニュートン法や最急降下法と組み合わせることで、より効率的に解を求めることができます。

物理シミュレーションへの応用

物理シミュレーションにおいて、ガウス・サイデル法は流体力学や熱伝導の問題を解くために使用されます。

これらの問題は、通常、連立方程式の形で表現されるため、ガウス・サイデル法を用いることで、シミュレーションの精度を向上させることができます。

特に、時間発展を伴う問題において、反復的に解を更新することができるため、リアルタイムシミュレーションにも適しています。

機械学習におけるガウス・サイデル法の利用

機械学習の分野でも、ガウス・サイデル法は特に線形回帰やロジスティック回帰の最適化に利用されます。

これらのモデルでは、パラメータの最適化が必要であり、ガウス・サイデル法を用いることで、効率的にパラメータを更新し、収束させることが可能です。

また、深層学習のトレーニングプロセスにおいても、勾配降下法の一部としてガウス・サイデル法が応用されることがあります。

これにより、モデルの学習速度を向上させることができます。

まとめ

この記事では、ガウス・サイデル法の基本的な概念から実装方法、応用例まで幅広く解説しました。

特に、連立方程式を解くための効率的な手法としての特長や、数値解析や物理シミュレーション、機械学習における具体的な利用方法についても触れました。

これを機に、ガウス・サイデル法を実際の問題解決に活用してみることをお勧めします。

2025-01-18更新日: 2025-01-18

ガウス・サイデル法とは

ガウス・サイデル法のアルゴリズム

アルゴリズムの基本的な流れ

反復計算の仕組み

収束判定の方法

初期値の選び方

Pythonでガウス・サイデル法を実装する

必要なライブラリ

係数行列と定数ベクトルの定義

初期値の設定

反復計算の実装

収束判定の実装

結果の表示

完全なサンプルコード

ガウス・サイデル法の実装例

2元連立方程式の例

3元連立方程式の例

収束しない場合の例

収束条件を満たす場合の例

ガウス・サイデル法の応用

大規模な連立方程式への適用

数値解析におけるガウス・サイデル法の利用

物理シミュレーションへの応用

機械学習におけるガウス・サイデル法の利用

まとめ

関連記事