feature image

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

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

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

Newton法(Newton-Raphson法)

newton2

目標

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

概要

を方程式の解に対する近似値とし、この方程式が区間に唯一の解を持つとします。を、座標におけるの接線と座標の交点とします。このとき

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

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

誤差の評価

実際に誤差がどのような振る舞いをするのか気になるので、具体的にみるためにまわりでTaylor展開することを考えます。
を含む開区間級なのでこの展開は可能で、次まで展開すると、

これにを代入すると、

変形して、

いま、

であるから(∵)、つ前の式を代入すれば

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

※注意※

グラフを見ると、ではなくから始めた場合は解の近似は得られません。また、の場合にもが求まりません。適用できるのは、
が区間とならず、の符号が異なり、が同符号である端点から始める場合」です。このときに含まれる唯一の解の近似が求まります。(力不足により証明略)

使用例

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

DKA法(Durand-Kerner-Aberth法)

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

謝罪

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

まとめ

参考文献

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

shisei_youki icon
この記事を書いた人
shisei_youki

この記事をシェア

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

関連する記事

2019年12月21日
モデリングを始めてみたい君へ、MagicaVoxelのススメ
isak icon isak
2019年12月13日
ゲーム紹介「League of Legends」【AdC2019 44日目】
Yataka_ML icon Yataka_ML
2019年12月25日
TensorFlow.jsでwasmを使ってみるためにコントリビュートした【AdC2019 56日目】
sappi_red icon sappi_red
2019年12月25日
無料でDTM環境構築 電子音系編
liquid1224 icon liquid1224
2019年12月4日
部内製チャットサービス「traQ」UIのこれまで 【AdC2019 35日目】
spa icon spa
2019年12月23日
無料でDTM環境構築 生音系編
kashiwade icon kashiwade
記事一覧 タグ一覧 Google アナリティクスについて 特定商取引法に基づく表記