一変数多項式の解を近似で求める求根アルゴリズムについて(9.5割Newton法です)
※数学的記述はおきもちであって正確でない場合があります(指摘お待ちしております)

Newton法(Newton-Raphson法)

目標

一変数実多項式$f(x)=0$(の係数)が与えられたときに、$|f(x_0)|$が十分小さい$x_0$を求める

概要

$x_0'$と$x_0$を方程式$f(x)=0$の解に対する近似値とし、この方程式が区間$(x_0',x_0)$に唯一の解$ξ$を持つとします。$x_1$を、座標$(x_0,f(x_0))$における$f(x)$の接線と$x$座標の交点とします。このとき$x_1$は

$\begin{array}{c} 0 = f'(x_0)(x_1-x_0)+f(x_0) \\ \Leftrightarrow x_1 = x_0 - \frac{f(x_0)}{f'(x_0)} \end{array}$

なので、一般に$1$以上の整数$n$に対して

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

であることが分かります。グラフを見ると、これを繰り返していけば目標が達成できそうです。

誤差の評価

実際に誤差がどのような振る舞いをするのか気になるので、具体的にみるために$f(x)$を$x_{n-1}$まわりでTaylor展開することを考えます。
$f(x)$は$ξ$を含む開区間$(x_0',x_{n-1})$で$C^{∞}$級なのでこの展開は可能で、$2$次まで展開すると、

$f(x) = f(x_{n-1}) + f'(x_{n-1})(x-x_{n-1}) + \frac{1}{2}f''(x-x_{n-1})(x-x_{n-1})^2+O((x-x_{n-1})^3)$

これに$x=ξ$を代入すると、

$f(ξ) = f(x_{n-1}) + f'(x_{n-1})(ξ-x_{n-1}) + \frac{1}{2}f''(x_{n-1})(ξ-x_{n-1})^2+O((ξ-x_{n-1})^3)$

変形して、

$-\frac{f(x_{n-1})-f(ξ)}{f'(x_{n-1})}=ξ-x_{n-1}+\frac{1}{2}\frac{f''(x_{n-1})}{f'(x_{n-1})}(ξ-x_{n-1})^2+O((ξ-x_{n-1})^3)$

いま、

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

であるから(∵$f(ξ)=0$)、$2$つ前の式を代入すれば

$\begin{array}{c} x_n-ξ=x_{n-1}-ξ+ξ-x_{n-1}+\frac{1}{2}\frac{f''(x_{n-1})}{f'(x_{n-1})}(ξ-x_{n-1})^2+O((ξ-x_{n-1})^3) \\ = \frac{1}{2}\frac{f''(x_{n-1})}{f'(x_{n-1})}(x_{n-1}-ξ)^2+O((ξ-x_{n-1})^3) \end{array}$

よって、誤差$x_n-ξ$は前の誤差$x_{n-1}-ξ$の$2$乗に比例して小さくなることが分かります(二次収束)。
つまり、$n$が$1$増える毎に正確な桁数がだいたい$2$倍になるということです←すごい

※注意※

グラフを見ると、$x_0$ではなく$x_0'$から始めた場合は解の近似は得られません。また、$f'(x_{n})=0$の場合にも$x_{n+1}$が求まりません。適用できるのは、
「$f'(x)$と$f''(x)$が区間$(x_0',x_0)$で$0$とならず、$f(x_0')$と$f(x_0)$の符号が異なり、$f(x)$と$f''(x)$が同符号である端点から始める場合」です。このとき$(x_0',x_0)$に含まれる唯一の解$ξ$の近似が求まります。(力不足により証明略)

使用例

• $\sqrt{α}$の近似
  これは$x^2=\sqrt{α}\Leftrightarrow x^2-\sqrt{α}=0$を解く事と同じなので、これについてNewton法を適用します
  $f(x)=x^2-\sqrt{α}$とすれば$f''(x)=2>0$より初期値を適当な大きい数にすれば適用できる条件を満たせます
  C++での実装例:

double newton_sqrt(double x,double eps=1e-8){
  if(x < 0)return -1;
  if(!x)return 0;
  double res = x < 1 ? 1 : x; //初期値
  while(res*res - x > eps){
    res = res - (res*res - x) / (2*res);
  }
  return res;
}

mt19997で$10^7$回実行したら1700msかかりました　尚 std::sqrt 17ms …

• $\sqrt[3]{α}$の近似
  さっきみたいに$f(x)=x^3-α=0$を解きます、$f''(x)=6x$の符号は$x$の符号と必ず同じなので$α<0$なら十分小さい負数を、$α>0$なら十分大きい正数を初期値にします
  C++での実装例:

double newton_cbrt(double x,double eps=1e-8){
  if(!x)return 0;
  double res = abs(x) < 1 ? (x < 0 ? -1 : 1) : x ; //初期値
  while(abs(res*res*res - x) > eps){
    res = res - (res*res*res - x) / (3*res*res);
  }
  return res;
}

mt19997で$10^7$実行すると2670ms ( std::cbrt 750ms STLすごい)

• $e^α$とか$\log(α)$とか三角関数の近似
  Taylor展開である程度の次数まで展開してやれば、これらも有限次数の多項式で表せるので、Newton法が使えます

DKA法(Durand-Kerner-Aberth法)

Newton法を改良し、複素数の解も含んだ一般の代数方程式にも適用できるようにしたもの

謝罪

DKA法を読んでください(丸投げ)

まとめ

• Newton法は実多項式の近似解を素早く求められる(例外アリ)
• Newton法とTaylor展開で色んな初等関数の求根もできる
• DKA法はすごい
• Wolfram Alpha は最高

参考文献

A Tutorial of BNCpack[暫定版]
数値計算の常識-Newton法
DKA法