1変数n次代数方程式の根n個を求める

数学 プログラミング

WikipediaによればDurand-Kerner法という数値解法があって（他にも色々あるようだが）n個の根を同時に求められるらしい．

簡単そうだったのでC++で実装してみた：

プリプロセッサ・シンボル DEBUG を定義してコンパイルすると中間状態を吐き出す：

C:\home\nodakai\prog\solving\numerical>cl /EHsc /DDEBUG durand_kerner.cpp
Microsoft(R) 32-bit C/C++ Optimizing Compiler Version 15.00.21022.08 for 80x86
Copyright (C) Microsoft Corporation.  All rights reserved.

durand_kerner.cpp
Microsoft (R) Incremental Linker Version 9.00.21022.08
Copyright (C) Microsoft Corporation.  All rights reserved.

/out:durand_kerner.exe
durand_kerner.obj

C:\home\nodakai\prog\solving\numerical>durand_kerner.exe
> 2 1
2 x + 1
(1.83697e-016,3)
(-0.5,0)
(-0.5,0)
(-0.5,0)
> 4 4 1
4 x^2 + 4 x + 1
(3.53553,3.53553) (-3.53553,-3.53553)
(1.25009,1.78544) (-1.61249,-1.1885)
(0.180321,0.713706) (-0.922176,-0.445941)
(-0.239485,0.274456) (-0.661103,-0.169894)
(-0.400456,0.104936) (-0.561527,-0.0648694)
(-0.461975,0.0400879) (-0.523501,-0.0247765)
(-0.485476,0.0153125) (-0.508976,-0.00946372)
(-0.494452,0.00584889) (-0.503429,-0.00361482)
(-0.497881,0.00223408) (-0.50131,-0.00138074)
(-0.499191,0.000853342) (-0.5005,-0.000527394)
(-0.499691,0.000325948) (-0.500191,-0.000201447)
(-0.499882,0.000124501) (-0.500073,-7.69458e-005)
(-0.499955,4.75551e-005) (-0.500028,-2.93907e-005)
(-0.499983,1.81644e-005) (-0.500011,-1.12262e-005)
(-0.499983,1.81644e-005) (-0.500011,-1.12262e-005)
> 1 0 0 -8
x^3 - 8
(7.79423,4.5) (-7.79423,4.5) (-1.65327e-015,-9)
(5.21261,2.97149) (-4.49902,3.12369) (-0.286034,-4.69175)
(3.34777,1.52822) (-2.4886,2.03586) (-0.499805,-2.68915)
(2.24571,0.599639) (-1.36696,1.62039) (-0.814283,-1.80674)
(1.95759,0.0628997) (-0.998038,1.70837) (-0.999563,-1.72734)
(1.9999,-0.00049874) (-1,1.73208) (-1,-1.73205)
(2,-1.88872e-009) (-1,1.73205) (-1,-1.73205)
(2,-2.18582e-021) (-1,1.73205) (-1,-1.73205)
(2,-2.18582e-021) (-1,1.73205) (-1,-1.73205)

定義しないと答えだけを表示する．

C:\home\nodakai\prog\solving\numerical>cl /EHsc durand_kerner.cpp
Microsoft(R) 32-bit C/C++ Optimizing Compiler Version 15.00.21022.08 for 80x86
Copyright (C) Microsoft Corporation.  All rights reserved.

durand_kerner.cpp
Microsoft (R) Incremental Linker Version 9.00.21022.08
Copyright (C) Microsoft Corporation.  All rights reserved.

/out:durand_kerner.exe
durand_kerner.obj

C:\home\nodakai\prog\solving\numerical>durand_kerner.exe
> 1 -3000 2999998 -999996000 -2999999 999999000 0
x^6 - 3000 x^5 + 3e+006 x^4 - 9.99996e+008 x^3 - 3e+006 x^2 + 9.99999e+008 x
(-1,7.97467e-177) (1.99195e-180,-1.61087e-180) (999,1.17013e-127) (1001,3.33413e-128) (1000,1.01609e-133) (1,1.89929e-183)
> 1 -3000 3000000 -1000000000 0 0 0
x^6 - 3000 x^5 + 3e+006 x^4 - 1e+009 x^3
(-1.90797e-112,7.19348e-112) (-1.96178e-112,-1.35934e-111) (1000,0.00266052) (1000,-0.000145985) (999.998,-0.00293325) (1.06013e-111,1.85072e-112)

x^6 - 3000 x^5 + 2999998 x^4 - 999996000 x^3 - 2999999 x^2 + 999999000 x = (x+1)x(x-1)(x-999)(x-1000)(x-1001) であって，1ずつ離れた根の3つ子が1000離れて存在している．ランダムな初期値＆Newton法ではまず全てを見つけることはできないだろう．またx^6 - 3000 x^5 + 3000000 x^4 - 1000000000 x^3 = x^3(x-1000)^3 であって，3重根2つを持つ．このような場合でも問題ない．
誤差の範囲内でよく動いてくれているようだ．