Programmierkurs
für Naturwissenschaftler/innen

ニュートン法

ニュートン法 (ニュートン・ラフソン (Newton-Raphson) 法) はある関数 $f(x)$ と $x$ 軸との交点を求める手法の一つである。 $f(x) = 0$ の解を求める、と言ってもよい。

例えば今 $\sqrt{2}$ の値を求めたいとする。 $\sqrt{2}$ は $$ x^2 = 2 $$ の解であるから、$f(x) = x^2 - 2$ としたときの $f(x) = 0$ の解を求めれば良いということになる1

考え方

今与えられた課題は $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}$ における $y = f(x)$ の接線の式を求めよ

(解答) $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$ 軸との交点 $x_n$ を求めよ

(解答) この接線と $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$ の解を求める手法であるということを考えると、後者の方が直観的であろう)。 「十分」というのも曖昧な表現であるが、これは必要な精度に依存する。

ニュートン法のアルゴリズム

では、具体的にニュートン法のアルゴリズムを導き出してみよう。 要は上で書いたことを、プログラムしやすいように言い替えるだけである。

ただし、ここで一つ注意すべき点がある。 これまで $x_n$ などと書いてきたので $x$ を配列に格納したくなるかもしれない。 しかし、$x_n$ を求めるのに必要なのは $x_{n-1}$ だけなので、古い $x$ は不要である。 すなわち、直前の $x$ だけ残しておけばよい。 以下では、直前の $x$ を覚えておくための変数として $x_\mathrm{prev}$ を使う。 $f(x)$, $f'(x)$ は定義されているものとする。

  1. $x_0$ を決める (与える)。これを (最初の) $x_\mathrm{prev}$ とする

  2. $x$ を $x_\mathrm{prev}$ から求める。すなわち $x = x_\mathrm{prev} - f(x_\mathrm{prev}) / f'(x_\mathrm{prev})$ を計算する

  3. 終了するための条件を満たしたか判定 (収束判定)

  4. まだ収束していない (と判断した) 場合、$x_\mathrm{prev} \leftarrow x$ として、2. に戻る。さもなくば次へ進む

  5. 終了

このアルゴリズムのフローチャートを描いてみよ

収束しなかった場合の処理

ここまでは収束することを前提に話を進めてきたが、実はニュートン法は必ず収束するとは限らない。 そのため、最悪の場合は停止しないことになる。 このような場合のために、最大の繰り返し回数を事前に決めておき、それを越えたら諦める、とすることが多い。

実習

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

$f(x) = x^2 - 2$ について、

  1. $f(x)$ 及び $f'(x)$ を求める関数を書け

  2. その関数を用いて、$f(2)$, $f'(2)$ を求めよ

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

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

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


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

上で定義した関数を用いて、ニュートン法で $\sqrt{2}$ を求める関数を書け。 その関数を用いて $\sqrt{2}$ を求めよ。

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

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


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

: 発展問題

上のプログラムは $\sqrt{2}$ しか求められないが、せっかくなので任意の $a$ について $\sqrt{a}$ を求めたいと思うであろう。 そこで、一般の $a$ ($a$ は 0 以上の実数) について $\sqrt{a}$ を求められるように上のプログラムを修正せよ。

ニュートン法の問題点

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

収束

既に述べたように、ニュートン法での計算は収束するとは限らない。 これは関数の性質 (及び初期値の取り方) に大きく依存する。

初期値

もう一つ重要な点は、解が複数ある場合でも解が一つしか求まらないこと、その解は初期値に依存すること、である (これが最初の脚注で述べた「大きな問題」である)。

例えば平方根を求めるプログラムで初期値を $-1$ などにすれば、負の値 ($-\sqrt{a}$) が求まるであろう。 これは正しい解であるが、今欲しいものではない。 平方根の場合は関数の性質が良く分かっているので、正の解を得るには初期値を正に取っておけば問題なさそうであることが分かる。

しかし、一般の関数では話はそう単純ではない。 解の個数が未知の場合はもっと話が複雑になる。 その関数がいくつ解を持つかはニュートン法では分からない。 そのため、解の個数が未知の場合、さまざまな初期値で試してみるなどの方法を取る必要がある。

一般的に、初期値は解に十分近い点を選ぶことが必要である。

導関数が必要

ニュートン法を適用するには、対象となる $f(x)$ だけでなく $f'(x)$ も必要であり、これが問題になる場合もある。 $f'(x)$ が簡単に表せないような場合では工夫が必要となる。

補足

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


1

気付いている人もいるであろうが、一つ大きな問題がある。これについては後で述べる。