Entries

D・K・A! D・K・A!

※某漫画の話ではありません

○あらすじ
 普通、方程式を数値計算で解くのには、ニュートン法二分法がよく用いられます。ですが、これらの方法だと、解を1つづつ求めることしかできません。つまり、大雑把にでもグラフを描いて(増減を調べて)から計算する必要があるのです。……ですが、N次方程式(a0xn+a1xn-1+a2xn-2+...+an-1x+an=0)の場合、N個の解(複素数含む)をまとめて求める方法が存在します。今回は、そのアルゴリズムについての話。

○ベアストウ(Bairstow-Hitchcock)法
 高校の教科書にも書いてあることですが、f(a)=0となる時に限り、f(a)は(x-a)で割り切れます。要するに、
  1. 方程式に代入して0になるxを探す
  2. 割り算して1.に戻る
を延々繰り返せば、どんなN次方程式でも解くことが可能です……が、それが簡単に通用するのは高校までです。
 整数や分数ならともかく、実数係数の方程式ともなると、手作業で特定するのは困難を極めます(ましてや解が複素数の場合だと(ry)。そこでコンピュータの出番となるわけですが、二分法(試しに代入してみる)やニュートン法(接線を利用する)を使用するにしても、できれば実数計算だけで済ませた方が速く済みます。そこで、2次式でとりあえず割ってみるのがこのベアストウ法なのです。
 幸い、係数が全て実数である場合、a+biが解ならば、その共役複素数であるa-biも解です(iは虚数単位)。つまり、
  1. 方程式を割り切る二次式である、x2+px+qを探す
  2. x2+px+q=0を解く。割り算して商が1次式or2次式ならそのまま解いて終了
  3. 1.に戻る
とすれば、複素数の計算を避けて計算できます。……と書くと、どうやってその二次式求めるんだよという(当然の)疑問が浮かんでくると思うのですが、そこはニュートン法で収束させることで乗り切ります。

☆BASIC(PC-G850Vというポケコン用に書いたもの)で書いたコードを公開します(ファイル名はBAIRSTW.BAS)。

2012/07/23 19:19追記:BAIRSTW.BASに修正がありました。この一行を追加してください。

295 IF N=2 THEN GOTO 160



○DKA(Durand-Kerner-Aberth)法
 前項のような方法は、係数に複素数が含まれる場合には使えません。そんな場合にも使えるのがDKA法です。

 a0xn+a1xn-1+a2xn-2+...+an-1x+an=0 において、解がA1・A2・A3・……Anとします。
 すると、因数定理により、方程式はa0(x-A1)(x-A2)...(x-An)=0と因数分解されます。
 この左辺をF(x)とおいて微分(積の微分法)し、xにAjを代入(1≦j≦n)すると、F'(x)=a0Π(1≦i≦n、ただしi≠jAj-Ai)といった式になります。で、これをニュートン法の式に入れれば、後は反復計算で収束するというわけです。

 以上がDurand-Kerner法と呼ばれているものですが、ここでも初期値はどうするのかといった問題がつきまといます。前述のソースコードではPとQ、両方を1として初期値にしましたが、こちらの場合、収束する値を上手くバラけさせる方法が必要になります。その方法として有名なのがAberthの初期値と呼ばれるものです。式は略しますが、要するに「(複素平面上での)円の周囲に均等に配置すればいいんじゃね?」といった発想です。この辺でどうしても「複素数」を計算する必要が生じるため、プログラミングする際には、配列なり構造体なりクラスなり、適当な方法で実装しましょう。

 まず、所持するポケコン用に書いたものがDKA.BASです。関数定義とか構造体とかが無い点を補うため、コードがやけに長くなっています。
 次に、STLの<complex>全開で書いたC++のコードがDKA法.cppです。BASICよりも自然、かつ短いコードで済むのがC++の威力ですね。
※input.txtの例:

4 0.00001
1 -1
-3 5
6 -12
-6 14
4 -8


○その他
 実行してみれば分かります(ポケコン持ってる人は少ないと思うけど)が、それなりにサクサクと動きます。C++版に至っては、10次方程式だろうが100次方程式だろうが一瞬で解を出力してくれます。ただ、DKA法において反復する計算式は、もうちょっと改良のしようがあります。Ehrlish-Aberth法田辺法なんかがそうです(参考文献)。また、SOR法(証券でもオカルトでもない)の手法を取り入れることもできます。DKA法便利すぎワロタwwwwと感じたので、ブログに書いてみました。

2013/06/20追記:
リンク切れを修正。
関連記事
スポンサーサイト
この記事にトラックバックする(FC2ブログユーザー)
http://ysrken.blog.fc2.com/tb.php/21-e625e9e3

トラックバック

コメント

コメントの投稿

コメントの投稿
管理者にだけ表示を許可する

Appendix

プロフィール

YSR

Author:YSR
「YSR」「YSRKEN」「◆YSRKENkO6Y(~2013/08/25)」「◆YSRKEN.ceVZZ(2013/08/26~)」として活動しています。
プログラミングと艦これが趣味です。
プロフ画像はCrystalDiskInfoの水晶雫ちゃんです。
主な創作物についてはhttp://ysrken.blog.fc2.com/blog-entry-76.htmlをご覧ください。

カレンダー(月別)

08 ≪│2017/09│≫ 10
- - - - - 1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30

全記事表示リンク

全ての記事を表示する

QRコード

QR

総アクセス数

アクセス数

現在の閲覧者数: