数値解析 第3回 (4) ニュートン法の誤差解析

ニュートン法の誤差解析

 テイラー展開を用いてニュートン法の誤差解析をしましょう。 初期値 $x_0$ が解 $\alpha$ に十分近いと仮定して議論を始めます。
  1. まず $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)$   $\cdots\cdots\ $ $(1)$
    ただし $O(t^3)$ は $t^3$ 以上の項を表す記号です。前回注意したように $O(t^3)$ はとても小さい値です。
  2. $f'(x)$ も $x=\alpha$ を中心にしてテイラー展開します。
    $f'(x)=f'(\alpha)+f''(\alpha)t+O(t^2)$   $\cdots\cdots\ $ $(2)$
  3. $\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 \quad \quad \quad \cdots\cdots\ (3) \\ \end{align} ただし $(\ast)$ では前回の教材で確認したべき級数展開
    $\dps{\frac{1}{1+X} = 1 - X + X^2 - X^3 + X^4 - \cdots}$
    を $X=\beta\,t+O(t^2)$ で用いています。
  4. $x_n$ の絶対誤差は $e_n = x_n - \alpha$ と表され、3° より \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{ ≒ }\ e_n - \left(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{ ≒ } \left(\frac{m-1}{m}\right)e_n$
となります。例えば3重解なら
$e_{n+1} \mbox{ ≒ } \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}}$
に収束しました。