理論と実装で理解するニュートン法|数値解析入門 #1

f:id:lib-arts:20200125192420p:plain

数値解析(Numerical Analysis)については所々取り扱っていると思いますが、数値解析として取り上げてはいなかったので数値解析の入門的な内容をまとめていければと思います。
#1では数値解析の基本的な手法であるニュートン法についてまとめつつ、Pythonを用いて簡単に実装を行えればと思います。

数値解析 - Wikipedia

ニュートン法 - Wikipedia

上記を参考に、数値解析やニュートン法について簡単に確認した後にニュートン法の実装を行ってみます。
以下、目次になります。
1. 数値解析について
2. ニュートン法の概要について
3. ニュートン法の実装
4. まとめ


1. 数値解析について
1節では数値解析についてWikipediaの記述を元に簡単に言葉の整理を行なっておきます。

f:id:lib-arts:20200125192805p:plain

f:id:lib-arts:20200125192743p:plain

数値解析 - Wikipedia

上記が冒頭部の記載ですが、数値解析(Numerical Analysis)は「代数学的な方法で解を得ることが不可能な解析学上の問題を、数値を用いて近似的に解く手法に関する学問」であるとされています。数値解析においては厳密な解を求めるというより、ある程度の誤差の範囲内の近似解を求めようとするとされています。

チェスなどのゲームの木の探索やセールスマン巡回問題、近年流行りのDeepLearningなどのアルゴリズムも現象を近似的にモデリングしようとしているので、多くの工学的なアプローチではここで述べた数値解析の枠組みを用いています。

数値解析については言葉の整理ができれば十分だと思われるので1節はここまでとします。


2. ニュートン法の概要について
2節では今回の本題であるニュートン法について確認します。

f:id:lib-arts:20200125193615p:plain

ニュートン法 - Wikipedia

上記がニュートン法(Newton's method)の概要で、「方程式系を数値計算によって解くための反復法による求根アルゴリズムの1つ」であるとされています。ニュートン法ニュートン・ラフソン法(Newton-Raphson method)とも呼ばれており、これはアイザックニュートンとジョセフラフソンに由来すると記載されています。

f:id:lib-arts:20200125192420p:plain

さて、具体的なアルゴリズムですが、上記の概念図のようにx_{n}をアップデートしていくことによって繰り返し的に方程式の解を求めます。数式としては、x_{n+1}=x_{n}-\frac{f(x)}{f'(x)}を用いてx_{n}を漸化式的に値を更新していきます。この数式の意味としては図と同様なのですが、(傾き)=(y方向の変化)/(x方向の変化)の数式をベースにしています。ここで傾きはf'(x)、y方向の変化はf(x)、x方向の変化はx_{n}-x_{n+1}にそれぞれ対応しており、これらを(傾き)=(y方向の変化)/(x方向の変化)の式に代入し変形することでx_{n+1}=x_{n}-\frac{f(x)}{f'(x)}を得ることができます。
ここで例としてf(x)=x^2-2があげられていますが、f(x)=x^2-2をおくとf'(x)=2xのため、x_{n+1}=x_{n}-\frac{x_{n}^2-2}{2x_{n}}=\frac{1}{2}(x_{n}+\frac{2}{x_{n}})となりこれを用いて繰り返し演算を行なっていくことになります。この辺の演算は実装を元に確認する方がわかりやすいので、計算については3節で取り扱います。


3. ニュートン法の実装
3節ではニュートン法を実際に実装で確認します。2節で取り扱ったf(x)=x^2-2の求根にあたってx_{0}=2とおいて計算してみます。

x = 2
for i in range(10):
    print(str(i)+":"+str(x))
    x = (x+2/x)/2

実行結果は下記のようになります。

f:id:lib-arts:20200125200113p:plain

x^2-2=0の解である\sqrt{2}=1.41421...が求まっていますが、初期値を変えて実験してみます。

x = 10
for i in range(10):
    print(str(i)+":"+str(x))
    x = (x+2/x)/2

f:id:lib-arts:20200125200144p:plain

x = 0.5
for i in range(10):
    print(str(i)+":"+str(x))
    x = (x+2/x)/2

f:id:lib-arts:20200125200225p:plain

x = -4
for i in range(10):
    print(str(i)+":"+str(x))
    x = (x+2/x)/2

f:id:lib-arts:20200125200303p:plain

上記のように初期値の設定によらずx^2-2=0の解に値が収束していることがわかります。また、この際に注意なのが、x^2-2=0二次方程式で二つの解を持つのですが、初期値から近い解に収束していることです。今回はy軸で対象な関数のx軸との関係性を方程式として考えているので、初期値が+の場合は\sqrt{2}、初期値が-の場合は-\sqrt{2}に収束するようになっています。


4. まとめ
#1ではニュートン法(ニュートンラフソン法)について取り扱いました。
#2ではニュートン法を元にした準ニュートン法について取り扱っていきます。