SOR法のPython実装

参考文献

数値計算の基礎と応用[新訂版]
数値解析学への入門
杉浦　洋(南山大学教授)　著
発行日 2009/12/10

参考ページ

準備

オンラインコンパイラを使用します。

ソースコード

sample.py

import numpy as np

# 関数vの定義
def v(x, y):
    return np.exp(x) * np.cos(y)

# 関数gの定義
def g(x, y):
    return v(x, y)

# メイン関数
def main():
    n = 20
    d = 1.0 / (n + 1)
    Pi = np.pi
    w = 2.0 / (1 + np.sin(Pi / (n + 1)))
    u = np.zeros((n+2, n+2))
    mu = 0

    # 境界条件の設定
    for i in range(n+2):
        u[i][0] = g(i*d, 0)
        u[i][n+1] = g(i*d, 1)
        u[0][i] = g(0, i*d)
        u[n+1][i] = g(1, i*d)
        mu += u[i][0] + u[i][n+1] + u[0][i] + u[n+1][i]
    mu /= 4*n

    # 初期値の設定
    for i in range(1, n+1):
        for j in range(1, n+1):
            u[i][j] = mu

    kmax = 300
    eps = 1.0e-8
    for k in range(1, kmax+1):
        rmax = 0
        for i in range(1, n+1):
            for j in range(1, n+1):
                r = (u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1]) / 4 - u[i][j]
                rmax = max(rmax, abs(r))
                u[i][j] += w * r
        if rmax <= eps:
            break

    print(f"n={n}, eps={eps:.2e}, rmax={rmax:.2e}")
    print(f"number of iterations k={k}, (kmax={kmax})")
    ermax = 0
    for i in range(1, n+1):
        for j in range(1, n+1):
            ermax = max(ermax, abs(u[i][j] - v(i*d, j*d)))
    print(f"Maximum Error={ermax:.2e}")

main()

実行結果

console

n=20, eps=1.00e-08, rmax=9.35e-09
number of iterations k=67, (kmax=300)
Maximum Error=4.12e-05