2019年11月28日 | ブログ記事

面倒な多項式はCPUに任せよう

shisei_youki

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

Newton法(Newton-Raphson法)

newton2

目標

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

概要

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

0=f(x0)(x1x0)+f(x0)x1=x0f(x0)f(x0)\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}

なので、一般に11以上の整数nnに対して

xn=xn1f(xn1)f(xn1)x_n = x_{n-1} - \frac{f(x_{n-1})}{f'(x_{n-1})}

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

誤差の評価

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

f(x)=f(xn1)+f(xn1)(xxn1)+12f(xxn1)(xxn1)2+O((xxn1)3)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=ξx=ξを代入すると、

f(ξ)=f(xn1)+f(xn1)(ξxn1)+12f(xn1)(ξxn1)2+O((ξxn1)3)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)

変形して、

f(xn1)f(ξ)f(xn1)=ξxn1+12f(xn1)f(xn1)(ξxn1)2+O((ξxn1)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)

いま、

xnξ=xn1ξf(xn1)f(xn1)=xn1ξf(xn1)f(ξ)f(xn1)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(ξ)=0f(ξ)=0)、22つ前の式を代入すれば

xnξ=xn1ξ+ξxn1+12f(xn1)f(xn1)(ξxn1)2+O((ξxn1)3)=12f(xn1)f(xn1)(xn1ξ)2+O((ξxn1)3)\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}

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

※注意※

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

使用例

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で10710^7回実行したら1700msかかりました 尚 std::sqrt 17ms …

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で10710^7実行すると2670ms ( std::cbrt 750ms STLすごい)

DKA法(Durand-Kerner-Aberth法)

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

謝罪

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

まとめ

参考文献

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

この記事を書いた人
shisei_youki

この記事をシェア

このエントリーをはてなブックマークに追加

関連する記事

2019年12月11日
円周率が無理数であることの証明【AdC2019 42日目】
Tarara
2019年12月11日
ぷよぷよってナンですか?【AdC2019 42日目】
arahi10
2019年12月10日
ブログ投稿ツイートを無理やり自動化した話【AdC2019 41日目】
xecua
2019年12月10日
理工系科目英語クラスについて【AdC2019 41日目】
kano
2019年12月9日
大学生活を豊かにする自己管理ツール集【AdC2019 40日目】
Deka
2019年12月8日
ハッカソンに参加して素材を作った話
NABE

活動の紹介

カテゴリ

タグ