1 Newton法

1.1 Newton法とは

ある方程式の解を数値的に求めるにはどうすればよいのだろうか。例えば

$$x^2-2=0$$ (1.1)

という $\pm\sqrt{2}=\pm 1.41421356\cdots$を解に持つ方程式を考える。ここで$y=x^2-2$のグラフを描き、解の近くの点、例えば$x=2$における接線を引いてみる(図1左)。

図 1: Newton法の考え方(左)と、接線による $x$ 軸との交点の変化(右)。

$y=x^2-2$の、点$(x_0,x_0^2-2)$における接線は

$$y-(x_0^2-2)=2x_0(x-x_0)$$

であるから、$x$軸との交点を$(x_1,0)$とすれば

$$-(x_0^2-2)=2x_0(x_1-x_0)$$

が成り立つので、$x_1$について解くと

$$x_1=x_0-\frac{x_0^2-2}{2x_0}=\frac{x_0^2+2}{2x_0}$$ (1.2)

となる。この式に$x_0=2$を代入すると$x_1=6/4=1.5$が得られ、まあまあ$\sqrt{2}$に近い。ここで再び$x=1.5$における接線を引き、$x_0=1.5$を代入して$x$軸との交点を求めると、今度は $x_1=4.25/3=1.41(66\cdots)$となり、$\sqrt{2}$に近くなっている。もう一度繰り返せばさらに$\sqrt{2}$に近い $1.41421(56\cdots)$が求まる。この操作を繰り返せば$\sqrt{2}$の良い近似値を求めることができそうである。

これを一般的に表せば次のようになる。$f(x)=0$の解を求めたいとすると、まずグラフ$y=f(x)$の$x_n$における接線の式は

$$y-f(x_n)=f'(x_n)(x-x_n)$$

となる。この接線と$x$軸との交点を$(x_{n+1},0)$とすれば

$$-f(x_n)=f'(x_n)(x_{n+1}-x_n)$$

という方程式が成り立つので、$x_{n+1}$は

$$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}$$ (1.3)

で与えられる。この漸化式を用いて$x_n$を次々に計算すれば、$f(x)=0$の解が求まりそうである。この方法をNewton法(Newton's Method)という。また

$$F(x)=x-\frac{f(x)}{f'(x)}$$ (1.4)

を用いて $x_{n+1}=F(x_n)$という操作を繰り返す事になるので、この$F(x)$を$f(x)$の反復関数(iteration function)と呼ぶ。具体的に$f(x)=x^2-c$を代入すれば

$$F(x)=x-\frac{x^2-c}{2x}=\frac{x^2+c}{2x}$$

が得られる。

求まりそうである、と言ったのは求まる事を証明していないからだが、実は十分に良い近似値から反復を繰り返すことで必ず解が求まることが知られている(次の章で証明する)。例えば$x^2+1=0$のような実数解を持たない方程式でも、$\pm i$の付近から反復を始めるとちゃんと$\pm i$が求まる。しかしこの「十分に」というのが曖昧で、どの程度近ければよいのかはそれぞれの$p(x)$に依存する。

また初期値がどの解に近いかによって当然求まる解も違ってくる。先の$f(x)=x^2-2$だと、図1右から考えて、$x_0>0$ならば$+\sqrt{2}$に、$x_0<0$ならば$-\sqrt{2}$に収束するようである。$x_0=0$だと、$x_1$以降は常に無限大となるので、全く意味がない。

1.2 収束性と収束の速さ

実際に数値計算を行う際には、どのくらい少ない回数の反復で必要な精度の近似値を求められるかが重要になってくる。これは近似値を真値との誤差で表すことで調べられる。

具体的に平方根を求める反復関数

$$F(x)=\frac{x^2+c}{2x}$$

の収束の様子を調べてみよう。$x_k$を

$$x_k=(1+\varepsilon )\sqrt{c}$$ (1.5)

と表すと$\varepsilon$が誤差であり、$F(x_k)$は

$$F(x_k)=\frac{((1+\varepsilon )\sqrt{c})^2+c}{2(1+\varepsilon... ...left(1+\frac{\varepsilon ^2}{2(1+\varepsilon )}\right)\sqrt{c}$$ (1.6)

となるので、誤差は$\varepsilon$から $\varepsilon ^2/2(1+\varepsilon )$になる。図2右

図 2: 平方根を求める反復関数(左、$c=2$)と、誤差の変化(右)。同じ$x$ $(\varepsilon )$座標で見たとき、1回の反復によって実線の値が破線の値になる。

からも分かるように、$\varepsilon$が正ならば1回の反復で誤差は約半分になる。$\varepsilon$が $-1$に近い場合には誤差が大きくなってしまうが、次の反復からは半分程度に減っていくので収束はする。

ここで重要なのは、 $\varepsilon \approx 0$の場合に

$$\frac{\varepsilon ^2}{2(1+\varepsilon )} \approx \frac{\varepsilon ^2}{2}$$ (1.7)

という近似式が成り立ち、$\varepsilon$の次数が2であるから1回の反復で正しい桁数が約2倍になる事である。具体的に$c=4,x_0=2.02$ $(\varepsilon =0.01)$とすれば、 $F(x_0) \approx 2.0001$ $(\varepsilon \approx 0.00005)$となり、正しい桁数が2倍に増えている。このような収束を2次収束という。

一般に、$\alpha$に収束する数列$\{a_n\}$において、ある定数$n_0$よりも大きい全ての$n$に対して

$$\vert a_{n+1}-\alpha\vert \le c\vert a_n-\alpha\vert^r$$ (1.8)

となるような定数$c$が存在するとき、この数列は$\alpha$に$r$次収束すると言う。また誤差 $\varepsilon _k$を、$\alpha=0$の場合は $\varepsilon _n=a_n$、それ以外の場合は

$$\varepsilon _n=\frac{a_n}{\alpha}-1 \quad もしくは \quad \varepsilon _n=1-\frac{a_n}{\alpha}$$ (1.9)

と定義する¹。このとき$\alpha \ne 0$ならば、$a_n$は誤差 $\varepsilon _n$を用いて

$$a_n=(1+\varepsilon _n)\alpha \quad もしくは \quad a_n=(1-\varepsilon _n)\alpha$$ (1.10)

の様に表される。$\varepsilon$を用いると、(1.8)は

$$\vert\varepsilon _{n+1}\vert \le c\vert\varepsilon _n\vert^r$$ (1.11)

となる。また先の例からも分かるように、$a_n$の正しい桁数は $\varepsilon _n$の上位の桁を埋める$0$の数程度である。

1.3 逆数

加減乗算のみで逆数を求めるには、方程式

$$x^{-1}-c=0$$ (1.12)

の解をNewton法で求めるとよい。この場合反復関数は

$$F(x)=x-\frac{x^{-1}-c}{-x^{-2}}=2x-cx^2$$ (1.13)

となる。誤差の評価をすれば

$$F\left(\frac{1-\varepsilon }{c}\right) =2\left(\frac{1-\vare... ...arepsilon )-(1-\varepsilon )^2}{c} =\frac{1-\varepsilon ^2}{c}$$ (1.14)

となるので2次収束である。

一般に逆数を$r$次収束で求める反復関数は、 $\varepsilon =1-cx$を用いて

$$F(x)=x(1+\varepsilon +\varepsilon ^2+ \ldots +\varepsilon ^{r-1})$$ (1.15)

となる。これは $x=(1-\varepsilon )/c$を実際に代入すれば

$$F(x)=\frac{1-\varepsilon }{c}(1+\varepsilon +\varepsilon ^2+ \ldots +\varepsilon ^{r-1}) =\frac{1-\varepsilon ^r}{c}$$ (1.16)

のように誤差の項が消えてしまう事からすぐに分かる。

逆数を求める場合には初期値に注意しなければならない。どんな場合でも誤差が$r$乗になるため、図3

図 3: 逆数を求める反復関数(左、$c=2$)と、誤差の変化(右)。2次収束から4次収束の場合を示した。

からも分かるように、 $\vert\varepsilon _0\vert=\vert 1-cx_0\vert<1 \Longleftrightarrow \vert x\vert<2/c)$となるようにしなければならない。

1.4 平方根

何度か出てきたように、平方根は反復関数

$$F(x)=\frac{x^2+c}{2x}$$

を用いて2次収束で求めることができる。ここで方程式$x-cx^{-1}=0$の解も$\pm\sqrt{c}$であるから、 $g(x)=x-cx^{-1}$の反復関数を求めると

$$G(x)=x-\frac{x-cx^{-1}}{1+cx^{-2}} =\frac{2cx}{x^2+c}$$ (1.17)

となり、別の反復関数が得られる。この場合誤差は

$$G((1-\varepsilon )\sqrt{c}) =\frac{2(1-\varepsilon )}{2-2\va... ...\varepsilon ^2}{2-2\varepsilon +\varepsilon ^2}\right)\sqrt{c}$$

となるから2次収束である。また $h(x)=1-cx^{-2}$の反復関数は

$$H(x) =x-\frac{1-cx^{-2}}{2cx^{-3}} =x+\frac{cx-x^3}{2c} =\frac{3cx-x^3}{2c}$$ (1.18)

となり、さらに別の反復関数が得られる。誤差は

$$H((1-\varepsilon )\sqrt{c}) =\frac{3(1-\varepsilon )-(1-\var... ...c{3}{2}\varepsilon ^2+\frac{1}{2}\varepsilon ^3\right)\sqrt{c}$$ (1.19)

となるから2次収束である。

これらの反復関数のうち、$G(x)$は役に立たない。何故なら$G(x)$の評価には乗算2回と除算1回が必要であるが、$F(x)$は乗算1回と除算1回で済むからである。$H(x)$の評価には乗算2回と除算1回が必要に見えるが、$c$は定数であるから、除算は反復の最初に1回行うだけでよい。そうすれば後は乗算3回で$H(x)$を評価できる。$F(x)$と$H(x)$のどちらを使ったほうが速いかは環境に依存する。

多倍長演算を行う場合には、もし$2c$が1桁なら$H(x)$には多倍長の除算がなくなり、常に$F(x)$より速く求まる。あるいは$H(x)$において$c$を$1/c$で置き換えると、

$$H(x)=x+\frac{x-cx^3}{2}$$ (1.20)

となって多倍長の除算がなくなる。この反復関数を使えば$1/\sqrt{c}$が得られるから、求まった解に$c$を掛ければ全く除算を用いずに$\sqrt{c}$が得られる。さらに、最後の反復において$y_n=cx_n$とおき、

$$y_{n+1}=y_n+\frac{x_n}{2}(c-y_n^2)$$ (1.21)

とすれば乗算の回数を減らすことができる。