数値解析 第3回 (4) ニュートン法の誤差解析
ニュートン法の誤差解析
テイラー展開を用いてニュートン法の誤差解析をしましょう。
初期値 $x_0$ が解 $\alpha$ に十分近いと仮定して議論を始めます。
- まず $f(x)$ を $x=\alpha$ を中心にしてテイラー展開します。$x=\alpha+t$ とすると、$f(\alpha)=0$ より
$$f(x)=f(\alpha+t)=f'(\alpha)t+\frac{1}{2}f''(\alpha)t^2+O(t^3) \tag{1}$$
ただし $O(t^3)$ は $t^3$ 以上の項を表す記号です。前回注意したように $O(t^3)$ はとても小さい値です。
- $f'(x)$ も $x=\alpha$ を中心にしてテイラー展開します。
$$f'(x)=f'(\alpha)+f''(\alpha)t+O(t^2) \tag{2}$$
- $\dps{\beta=\frac{f''(\alpha)}{f'(\alpha)}}$ とおくと、$(1)$ ÷ $(2)$ から
\begin{align}
\frac{f(x)}{f'(x)}
&= \frac{f'(\alpha)t+\frac{1}{2}f''(\alpha)t^2+O(t^3)}{f'(\alpha)+f''(\alpha)t+O(t^2)} \\
&= \frac{t+\frac{1}{2}\beta\,t^2+O(t^3)}{1+\beta\,t+O(t^2)} \\
&= \left(t+\frac{1}{2}\beta\,t^2+O(t^3)\right) \times \left(1-\beta\,t+O(t^2)\right)\quad \cdots\cdots\ (\ast) \\
&= t-\frac{1}{2}\beta\,t^2+O(t^3)
\ \mbox{ ≒ }\ t - \frac{1}{2}\beta\,t^2 \tag{3} \\
\end{align}
ただし $(\ast)$ では前回の教材で確認したべき級数展開
$\dps{\frac{1}{1+X} = 1 - X + X^2 - X^3 + X^4 - \cdots}$
を $X=\beta\,t+O(t^2)$ で用いています。
- $x_n$ の絶対誤差は $e_n = x_n - \alpha$ と表されますので、(3) に $x = x_n= \alpha + e_n$ を代入することにより
\begin{align}
e_{n+1} = x_{n+1} - \alpha
&= \left( x_n - \frac{f(x_n)}{f'(x_n)}\right) - \alpha \\
&= (x_n - \alpha) - \frac{f(\alpha+e_n)}{f'(\alpha+e_n)} \\
&\mbox{ ≒ }\ \require{cancel}\cancel{e_n} - \left(\cancel{e_n} - \frac{1}{2}\beta\, {e_n}^2\right) = \frac{1}{2}\beta\, {e_n}^2\\
\end{align}
従って
$\dps{e_{n+1} \mbox{ ≒ } \frac{1}{2}\beta\, {e_n}^2}$
まとめると
Th.9 $x_0$ が解 $\alpha$ に十分近いとき、
$\dps{\beta=\frac{f''(\alpha)}{f'(\alpha)}}$ とおくと、
ニュートン法の数列 $x_n$ の絶対誤差 $e_n = x_n - \alpha$ は
$\dps{e_{n+1} \mbox{ ≒ } \frac{1}{2}\beta\, {e_n}^2}$
を満たす。
この式は、精度の「桁数」が倍々で増えてゆくことを示しています。
例えば $|\,\frac{1}{2}\beta\,| \leqq 1$ であるような状況で、
初期値の絶対誤差が $|\,e_0\,| \leqq 10^{-1}$ ぐらい小さければ、おおよそ
\begin{align}
|\,e_1\,| &\leqq 1 \times (10^{-1})^2 \leqq 10^{-2} \\
|\,e_2\,| &\leqq 1 \times (10^{-2})^2 \leqq 10^{-4} \\
|\,e_3\,| &\leqq 1 \times (10^{-4})^2 \leqq 10^{-8} \\
|\,e_4\,| &\leqq 1 \times (10^{-8})^2 \leqq 10^{-16} \\
\end{align}
のように誤差が減少してゆきます。
ニュートン法についての補足
Rem.10 上の議論は「初期値 $x_0$ が解 $\alpha$ に十分近い」ことが前提です。
しかし $\alpha$ がわからない状況で $x_0$ を適当に与えてしまうと
- $x_0$ に最も近い解が求まる保証は無いし、
- そもそも収束する保証もありません。
そこで、ループ回数に上限を決めることが必要になります。
Rem.11 解 $\alpha$ が重解のときは先に述べたように $f'(\alpha)=0$ となるので、
$\dps{\beta=\frac{f''(\alpha)}{f'(\alpha)}}$ が考えられず、上の議論が成り立ちません。
一般に $\alpha$ を $m$ 重解とすると、絶対誤差 $e_n$ の評価式は
$e_{n+1} \mbox{ ≒ } \dps{\left(\frac{m-1}{m}\right)e_n}$
となります。例えば3重解なら
$e_{n+1} \mbox{ ≒ } \dps{\left(\frac{2}{3}\right)e_n}$
となって二分法より収束が遅くなります。
2重解とわかっているのなら「 $f'(x)=0$ の解を求めて、その中から $f(x)=0$ の解を探す 」などの工夫が必要でしょう。
Rem.12 ニュートン法の収束原理はテイラー展開にあるので、
$f(x)$ が複素関数でも成り立ちますし、
初期値 $x_0$ を複素数にすると $f(x)=0$ の複素数解を求めることもできます。
Ex.13 $f(x)=x^3-1$ のとき、$x_0=i$ ( 虚数単位 ), 許容誤差 $=10^{-12}$ として実行すると
x0 = 1i
x1 = (-0.3333333333333333+0.6666666666666667i)
x2 = (-0.5822222222222222+0.9244444444444444i)
x3 = (-0.5087908032893191+0.8681655118873491i)
x4 = (-0.5000687390673927+0.8659822186925401i)
x5 = (-0.49999999628902975+0.8660253983385868i)
x6 = (-0.5+0.8660254037844387i)
x7 = (-0.5+0.8660254037844386i)
となり、ループ7回で虚数解
$\dps{\exp\left({\frac{2\pi}{3}i}\right)=\cos\left(\frac{2\pi}{3}\right)+i\sin\left(\frac{2\pi}{3}\right)=\frac{-1+\sqrt{3}\,i}{2}}$
に収束しました。