optimize-newton-raphson1d.py 技術ドキュメント

プログラムの動作

optimize-newton-raphson1d.py は、与えられた1次元関数の極値（最小値または最大値）をニュートン・ラプソン法を用いて探索するPythonプログラムです。探索の過程はグラフとしてリアルタイムに可視化され、各反復での探索点の移動と収束状況が標準出力にログとして表示されます。

主な機能は以下の通りです。

1次元関数の極値探索: 指定された関数 $f(x)$ の1階微分 $f'(x)$ がゼロになる点、すなわち極値をニュートン・ラプソン法で探索します。
数値微分: 関数の1階微分および2階微分は、中心差分法を用いた数値微分によって計算されます。
リアルタイム可視化: Matplotlibを使用して、対象となる関数 $f(x)$ のグラフ上に、ニュートン・ラプソン法の各反復で計算された探索点がプロットされ、視覚的に収束過程を追うことができます。
収束判定: 連続する探索点間の変化量 dx が所定の閾値 eps を下回った場合に収束と判断し、探索を停止します。
初期値設定: コマンドライン引数で探索の初期推定値 $x_0$ を指定できます。

このプログラムは、数値最適化アルゴリズムの一つであるニュートン・ラプソン法の動作原理を理解し、その収束過程を視覚的に確認するための教育的ツールとして役立ちます。

原理

このプログラムは、1次元関数 $f(x)$ の極値（最小値または最大値）を探索するために、ニュートン・ラプソン法を適用しています。関数の極値は、その1階微分 $f'(x)$ がゼロとなる点 ($f'(x)=0$) で発生します。したがって、この問題は $g(x) = f'(x)$ の根を見つける問題に帰着されます。

ニュートン・ラプソン法の更新式は以下の通りです。 $$x_{n+1} = x_n - \frac{g(x_n)}{g'(x_n)}$$ ここで $g(x) = f'(x)$ なので、$g'(x) = f''(x)$ となります。よって、極値探索のためのニュートン・ラプソン法の更新式は次のようになります。 $$x_{n+1} = x_n - \frac{f'(x_n)}{f''(x_n)}$$

プログラムでは、以下の関数が使用されます。

対象関数 func(x): 探索対象となる関数 $f(x)$ です。現在の実装では次の多項式関数が定義されています。 $$f(x) = -3.0 + x - 30.0 x^2 + 1.5 x^3 + 3.0 x^4$$
1階微分 diff1(x): $f(x)$ の1階微分 $f'(x)$ を計算します。プログラムでは中心差分法による数値微分が用いられています。微小区間幅を $h$ とすると、以下の式で近似されます。 $$f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$$
2階微分 diff2(x): $f(x)$ の2階微分 $f''(x)$ を計算します。プログラムでは中心差分法による数値微分が用いられています。微小区間幅を $h$ とすると、以下の式で近似されます。 $$f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}$$

プログラム内で定義されている dump パラメータは、2階微分 $f''(x)$ がゼロに近づき、更新式の分母がゼロになることを防ぐための安定化項として機能する可能性があります。具体的には、$f''(x)$ が負の場合は dump を引き、正の場合は dump を加えることで、絶対値をわずかに増加させます。ただし、デフォルト値は 0.0 のため、現在の実行ではこの安定化機能は無効です。

アルゴリズムは以下のステップで進行します。

初期推定値 $x_0$ を設定します。
各反復 $n$ で、現在の $x_n$ における $f(x_n)$、$f'(x_n)$、$f''(x_n)$ を計算します。
ニュートン・ラプソン法の更新式を用いて次の探索点 $x_{n+1}$ を計算します。
$|x_{n+1} - x_n|$ が収束閾値 eps 未満になった場合、または最大反復回数 nmaxiter に達した場合に探索を終了します。
各反復で計算された探索点をグラフ上にプロットし、アニメーションとして表示します。

必要な非標準ライブラリとインストール方法

このプログラムは以下の非標準ライブラリを使用しています。

NumPy: 数値計算を効率的に行うためのライブラリ。
Matplotlib: グラフ描画のためのライブラリ。

これらのライブラリは、Pythonのパッケージマネージャ pip を使ってインストールできます。

pip install numpy matplotlib

なお、プログラムコードには import csv が含まれていますが、現在の実装では csv ライブラリは使用されていません。

必要な入力ファイル

このプログラムは、実行時に特定の入力ファイルを必要としません。探索の初期推定値 $x_0$ は、コマンドライン引数として与えるか、プログラム内のデフォルト値（x0 = -2.0）を使用します。

生成される出力ファイル

このプログラムは、実行時にファイルシステムに新しいファイルを生成しません。実行結果は以下の形式で出力されます。

標準出力: 各反復での探索点の値、移動量 dx、関数値 f(x)、および収束状況がテキスト形式で表示されます。
グラフィカルユーザーインターフェース (GUI): Matplotlibによって生成されるグラフィックウィンドウに、対象関数の曲線と、ニュートン・ラプソン法によって逐次計算される探索点がアニメーションで表示されます。このウィンドウはプログラムが終了するか、ユーザーがキー入力をするまで表示され続けます。

コマンドラインでの使用例 (Usage)

基本的な実行コマンドは以下の通りです。

python optimize-newton-raphson1d.py [x0_initial_guess]

python optimize-newton-raphson1d.py: 初期推定値 x0 にプログラム内で定義されたデフォルト値（-2.0）を使用します。
[x0_initial_guess]: 任意引数です。ここに浮動小数点数を指定することで、探索の初期推定値を変更できます。

コマンドラインでの具体的な使用例

1. デフォルトの初期値で実行する場合

python optimize-newton-raphson1d.py

実行結果の説明: このコマンドを実行すると、プログラムはデフォルトの初期推定値 $x_0 = -2.0$ からニュートン・ラプソン法を開始します。標準出力には以下のようなログが表示され、各反復での $x$ の値、次の $x$ の値、変化量 $dx$、および関数値 $f(x)$ が報告されます。

Find minimum / maximum point by Newton-Raphson method

Iter     0: x:  -2.000000000000 =>  -1.350849780092, dx =  0.6492   f=  29.000000000000
Iter     1: x:  -1.350849780092 =>  -1.071850785055, dx =  0.2790   f=   8.916894334057
Iter     2: x:  -1.071850785055 =>  -1.006935105221, dx =  0.06492  f=   3.116630090875
Iter     3: x:  -1.006935105221 =>  -1.000049445100, dx =  0.006886 f=   2.051919864273
Iter     4: x:  -1.000049445100 =>  -1.000000000000, dx =  4.945e-05 f=   2.000395560410
Iter     5: x:  -1.000000000000 =>  -1.000000000000, dx = -2.715e-14 f=   2.000000000000
  Success: Convergence reached: dx = -2.715105151042769e-14 < eps = 1e-10
Press enter to terminate:

同時にMatplotlibのグラフウィンドウが開かれ、関数 $f(x)$ の曲線上に、初期点から最終的な収束点（この例では $x = -1.0$ 付近の極小値）へと移動する青い丸（探索点）がアニメーションで表示されます。最終的な収束点における $f(x)$ の値は $2.0$ となります。

2. 別の初期値で実行する場合

python optimize-newton-raphson1d.py 3.0

実行結果の説明: このコマンドでは、初期推定値 $x_0 = 3.0$ を指定して実行します。プログラムは指定された初期値からニュートン・ラプソン法を開始し、標準出力には同様のログが表示されますが、探索経路は異なります。

Find minimum / maximum point by Newton-Raphson method

Iter     0: x:   3.000000000000 =>   1.932857142857, dx = -1.067   f= 111.000000000000
Iter     1: x:   1.932857142857 =>   1.196303273449, dx = -0.7366  f=  32.898688401389
Iter     2: x:   1.196303273449 =>   0.803738012170, dx = -0.3926  f=   7.828552697843
Iter     3: x:   0.803738012170 =>   0.720084511529, dx = -0.08365 f=   4.045053150530
Iter     4: x:   0.720084511529 =>   0.713837424151, dx = -0.006247 f=   3.999602058525
Iter     5: x:   0.713837424151 =>   0.713832049615, dx = -5.374e-06 f=   3.999554030737
Iter     6: x:   0.713832049615 =>   0.713832049580, dx = -3.493e-11 f=   3.999554030737
  Success: Convergence reached: dx = -3.492582852207901e-11 < eps = 1e-10
Press enter to terminate:

グラフウィンドウでは、初期点 $x=3.0$ から別の極小値（この例では $x \approx 0.71$ 付近）へと探索点が移動する様子がアニメーションで表示されます。最終的な収束点における $f(x)$ の値は $3.999...$ となります。このように、ニュートン・ラプソン法は初期値によって収束する極値が異なる場合があります。