"""
自己無撞着法（不動点反復法）を用いて方程式の解を求めるスクリプト。
このモジュールは、`x = f(x)` の形式で表現された方程式を、反復式 `x_{n+1} = f(x_n)` を用いて解きます。
初期値から開始し、指定された収束条件または最大反復回数に達するまで計算を繰り返します。
混合比を用いて反復を安定化させることも可能です。
:doc:`equation-selfconsistent_usage`
"""
import csv
import numpy as np
from numpy import sin, cos, tan, pi
import sys
x0 = 0.0
kmix = 1.0
eps = 1.0e-10
nmaxiter = 50
[ドキュメント]
def func(x):
"""
自己無撞着法で解く方程式の右辺を定義します。
この関数は、`x = f(x)` 形式の方程式における `f(x)` に相当します。
具体的には、`f(x) = 0.25 * (-x^3 + x^2 - 2.0)` を計算します。
:param x: float 独立変数。
:returns: float 計算された `f(x)` の値。
"""
return 0.25 * (-x*x*x + x*x - 2.0)
[ドキュメント]
def main():
"""
自己無撞着法を実行し、方程式の解を求めます。
グローバル変数 `x0` (初期値)、`kmix` (混合比)、`eps` (収束判定条件)、
`nmaxiter` (最大反復回数) を使用して、反復計算を行います。
計算の過程と結果は標準出力に表示されます。
:returns: int 収束した場合は1、収束しなかった場合は0。
"""
global x0, kmix, eps, nmaxiter
print("Solution of an equation by self-consistent method")
x = x0
for i in range(nmaxiter):
xnext = func(x)
dx = xnext - x
print("Iter {:5d}: x: {:>16.12f} => {:>16.12f}, dx = {:>10.4g}".format(i, x, xnext, dx))
if abs(dx) < eps:
print(" Success: Convergence reached: dx = {} < eps = {}".format(dx, eps))
return 1
x = (1.0 - kmix) * x + kmix * xnext
print(" Failed: Convergence did not reach: dx = {} > eps = {}".format(dx, eps))
return 0
if __name__ == "__main__":
main()