APL/J言語:多項式(カーナー法による係数のルート)

APL/J言語:多項式(カーナー法による係数のルート)
カーナー法(Kerner's Method)という名前ではウィキペディアにはありませんが、「求根アルゴリズム」というページに若干の記述がありました。

      • -

ニュートン法は一度に一つのルートを対象とし、最初の近似が近いところにある必要があります。カーナー法では、すべてのルートを求める一般的な方法であり、リストから始め、残るf aの要素を対応するルートについての微分係数で割っていきます。この方法は最高次の多項式係数が1のものにしか対応しないので、まず最初に多項式を最後のもので割って正規化し、多項式が同じルートを持つようにします。この方式は当初の近似の最低限一つが複素数のとき、複素数に収束します。われわれは指数関数についてテイラー級数近似を用います。なぜなら対応する多項式複素数のルートを持つからです。

   ]d=: ^ t. i.6
1 1 0.5 0.166667 0.0416667 0.00833333

   ]c=: (norm=: % {:) d
120 120 60 20 5 1

   +. a=: (init=: r.@}.@i.@#) c |a
 0.540302  0.841471 1 1 1 1 1
_0.416147  0.909297
_0.989992  0.14112
_0.653644 _0.756802
 0.283662 _0.958924

   deriv=: [: */ 0&=@{.}@(-/~ ,: 1:)

   kerner=: 1 : '] - x&p. % deriv@]'

   r=: c kerner ^:_ a
   +.(/:|) r     NB.大きさによってソートされたルートの実数部分と虚数部分
_2.18061 4.57601e_31
 _1.6495     1.69393
 _1.6495    _1.69393
0.239806     3.12834
0.239806    _3.12834

   >./|c p. r
1.04488e_13

以上の結果は正規化しない場合の多項式係数dについての原始動詞p.(ピードット)の結果と比較できます。

   p. 2 4 2
 +-+-----+
 |2|_1 _1|
 +-+-----+

   ,.;}.p. d
 0.2398064j3.12834
0.2398064j_3.12834
   _1.6495j1.69393
  _1.6495j_1.69393
          _2.18061

ニュートン法複素数のルートに使うこともできます。

   +. d NEWTON ^:0 1 2 3 _ a=: 1j1
        1        1
0.0166065  0.99639
_0.990523 0.992532
 _1.95338  1.10685
  _1.6495   1.6939