ニュートン法

今与えられた課題は $f(x) = 0$ を満たす $x$ を求めることである。 その為に以下の手順で行う。

1. まず、$x$ 軸上のある点 $x_0$ を選ぶ。

2. $x = x_0$ での $f(x)$ の接線、つまり、点 $(x_0, f(x_0))$ での接線を求める。

3. その接線と $x$ 軸との交点を求める。 このときの $x$ の値を $x_1$ とする。

4. 同様にその点における接線を求め、$x$ 軸との交点を求め、… を繰り返す。

そうすると $x_n$ はどんどん求める解に近づいていく。

次のデモで試してみて欲しい。

基本的な考えかたはつかめただろうか。 次はこれを具体的なプログラムにする必要がある。

しかし、その前に次の問題を解いて欲しい。 数式は出てくるが、$_n$ とか $_{n-1}$ とかのせいでゴチャゴチャして見えるだけで、内容は非常に簡単 (接線と交点を求めるだけ) である。

(解答) $x = x_{n-1}$ で

1. $f(x)$ の値は $f(x_{n-1})$ であり、

2. $f(x)$ の接線の傾きは $f'(x_{n-1})$ である

から、求めるものは、「点 $(x_{n-1}, f(x_{n-1}))$ を通り、傾きが $f'(x_{n-1})$ であるような直線の式」となる。

求め方はいくつかあるが、一般に点 $(X, Y)$ を通る傾き $m$ の直線は $y - Y = m (x - X)$ と書けるから、 \[ y - f(x_{n-1}) = f'(x_{n-1}) (x - x_{n-1}) \] となる。

あるいは、$y$ 切片を $b$ とおいて、 \[ y = f'(x_{n-1}) x + b \] 点 $(x_{n-1}, f(x_{n-1}))$ を通ることから、 \[ f(x_{n-1}) = f'(x_{n-1}) x_{n-1} + b \] なので、$b$ を求めると、$b = f(x_{n-1}) - f'(x_{n-1}) x_{n-1}$。 従って求める式は、 \[ y = f'(x_{n-1}) x + f(x_{n-1}) - f'(x_{n-1}) x_{n-1} \] である、としても良い。

(解答) この接線と $x$ 軸との交点は、上で求めた式に $y = 0$ を代入すれば良いから、 \[ 0 - f(x_{n-1}) = f'(x_{n-1}) (x - x_{n-1}) \] を $x$ について解けば良く、 \[ x = x_{n-1} - f(x_{n-1}) / f'(x_{n-1}) \] となり、これが $x_n$ である。 \[ x_n = x_{n-1} - f(x_{n-1}) / f'(x_{n-1}) \]

もう一度原点に立ち返ろう。

考え方 で見たように $x_n$ で $n \rightarrow \infty$ とすると、$x_n$ は求める解に近づいていく。 ということは、まず $x_n$ の式を求めたい訳だが、それは直前ですでに見たように \[ x_n = x_{n-1} - f(x_{n-1}) / f'(x_{n-1}) \quad (n \ge 1) \] となる。 $n = 1$ のとき、 \[ x_1 = x_0 - f(x_0) / f'(x_0) \] であるから、$x_0$ を決めれば、$x_1$ が決まり、それによって $x_2$ が計算でき…、となる。 よってまずは $x_0$ を適当に決めればよい。 ($f(x)$ は与えられており、それにより $f'(x)$ が計算できるものと仮定していることに注意)。

上で接線の式を求めたが、それは $x_{n-1}$ から $x_n$ を求めるのに必要だっただけで、いったん上の式が求まってしまえば もう必要ないことにも注意しよう。

さて、ここで一つ大きな問題が出てくることに気付いたであろう。 数学的には $n \rightarrow \infty$ と一言で言えるが、実際に無限回計算する訳にはいかない。

アルゴリズムとは、堅苦しく言うと、「ある計算(あるいは処理)を実現するための手順」である。 これまでは、アルゴリズムを導き出す作業をしていた訳である。

アルゴリズムがアルゴリズムであるためには、対象の計算結果が正しく得られることも重要であるが、それに加えて「停止すること」が重要である。 これまでの議論では $n \rightarrow \infty$ とする必要があったが、これをプログラミングして実行しても永遠に終わらない。

実に当たり前のことを言っているように思うかもしれないが、この点は大変重要なので、ぜひ意識して欲しい。

さて、ニュートン法において、いつ停止させるか、という問題は実はそれほど簡単ではない。 その為の方法はいくつかあるが、一般的には $x_n$ と $x_{n-1}$ との差が十分小さくなったら、 あるいは、$f(x_n)$ の値が十分 0 に近づいたら、とすることが多い (ニュートン法が $f(x) = 0$ の解を求める手法であるということを考えると、後者の方が直観的であろう)。 「十分」というのも曖昧な表現であるが、これは必要な精度に依存する。

ここまでで材料はそろったはずである。 実際にプログラムを書いてみよう。

おおむね、以下のような感じになる (関数名等は変えてもよい):

def f(x):
    ...
    return ...

def df(x):
    ...
    return ...


print(f(2), df(2))

おおむね、以下のような感じになる (関数名等は変えてもよい):

def newton():
    ...
    return ...


root = newton()
print(root)
# or:
# print(newton())

ニュートン法は非常に優れた手法であるが、欠点がない訳ではない。

$f(x) = 0$ の解を求めるのではなく、ある関数 $g(x)$ の極値 (を与える $x$) を求めたい、という場合もよくある。 この場合は、求めるべきものは $g'(x) = 0$ の解であるから、ニュートン法を使うのであれば、 $$\begin{align*} f(x) & = g'(x) \\ f'(x) & = g''(x) \end{align*}$$ としてニュートン法を適用すれば良いことが分かる。