ベアストウ法

過去ブログの転載です。

高次方程式の数値解を求めるうまい手法に、ベアストウ法というのがあります。

数値解を得るのは一般的にニュートン法で十分じゃないかと思われますが、虚数解を持つ高次方程式のすべての解をニュートン法で得るのは面倒なことが多くてあんまりやりたくなるなるものです。ベアストウ法は虚数解を含めてすべての解を取得するまでのアルゴリズムがサポートされており、実装がラクだという特徴があります。

n次方程式が次のように表されているものとします。
x^n＋a[1]x^(n-1)＋…＋a[n-1]x＋a[n]＝0
a[1]～a[n]は何か係数です。この左辺を2次式x^2＋px＋qで割ると、
(x^2＋px＋q)(x^(n-2)＋b[1]x^(n-3)＋…＋b[n-3]x＋b[n-2])＋Rx＋S
という形になります。
大きい括弧の部分が商で、Rx＋Sは余りですね。
もし余りが0、つまりRx＋S＝0になりましたら、これは即ち因数分解に成功したことを意味します。
そうすると、少なくともx^2＋px＋q＝0ですから、これは普通の2次方程式なのでx＝－p±√(p^2－4q)／2として解を得ることができます。場合分けすれば虚数解でも代数演算で普通に求められますよね。

多分、この手法の狙いはここだと思います。これにより虚数解をラクに取得できます。で、残った商の部分でも同じように繰り返せば、すべての解が得られるという仕組みです。

問題なのは、どうやって余りを0にするかですが。

記号を使った一般的なお話は別のサイトで見ていただくことにしてここでは何かてきとうな例を挙げてそれをベアストウ法で解いてみることをします。

x^3－5x^2＋11x－15＝0 の解を複素数の範囲ですべて求めよ

高校生なら、解けなくちゃまずい問題ですけどこれをコンピュータに解いてもらうことを考えましょう。

まず、x^3－5x^2＋11x－15を2次式x^2＋px＋qで割って、
(x^2＋px＋q)(x＋b)＋Rx＋S＝0
という形にできたとしましょう。
R＝S＝0にしたいわけですが、RとSはpとqに依存していますので、それぞれR(p, q)、S(p, q)という二変数関数で表すことができます。

ベアストウ法では、このRとSをまず全微分します。
dR＝∂R／∂p dp＋∂R／∂q dq
dS＝∂S／∂p dp＋∂S／∂q dq
微分の定義式から、p, qの微小な増分Δp, Δqとすると、上式は
R(p＋Δp, q＋Δq)－R(p, q)≒∂R／∂p Δp＋∂R／∂q Δq
S(p＋Δp, q＋Δq)－S(p, q)≒∂S／∂p Δp＋∂S／∂q Δq
と書くことができますね。この際≒は面倒なので以後＝とします。

で、何がしたかったかというとR＝S＝0にしたかったわけですから、もとのR(p, q)、S(p, q)から増分Δp, Δqを与えることによって
R(p＋Δp, q＋Δq)＝0
S(p＋Δp, q＋Δq)＝0
になって欲しいので、上式を移項して代入して
∂R／∂p Δp＋∂R／∂q Δq＝－R(p, q)
∂S／∂p Δp＋∂S／∂q Δq＝－S(p, q)
というのが成り立てば良いということがわかります。
これはよく見るとΔpとΔqに関する連立方程式になっています。
R(p, q), S(p, q), ∂R／∂p, ∂R／∂q, ∂S／∂p, ∂S／∂q
の6つが分かれば、Δp, Δqが求まることになります。
求まったΔp, Δqをそれぞれp, qに足したものを新しいp, qにすればいいのですね。なにかてきとうな初期値p, qを決めて、それを更新していけばいいのです。では次にこれら6つのパラメータを求める手続きを考えましょう。

まず、R(p, q), S(p, q)を求めます。面倒なのでR, Sと書きます。
これは割りと単純に、
(x^2＋px＋q)(x＋b)＋Rx＋S
を展開して係数比較して求めればOKです。展開すると
x^3＋(b＋p)x^2＋(bp＋q＋R)x＋(qb＋S)
になりますから、x^3－5x^2＋11x－15と係数比較することで、
b＋p＝－5
bp＋q＋R＝11
qb＋S＝－15
という3本の式が出てきます。

p, qは既知ですから、それを定数として求めると
b＝－p－5
R＝p^2＋5p－q＋11
S＝pq＋5q－15
というのが求まると思います。

これでR, Sはp, qから計算できることがわかりました。

次に∂R／∂p, ∂R／∂q, ∂S／∂p, ∂S／∂qを求めましょう。
何のことはない、R, Sがわかってるので偏微分すればいいだけです。
∂R／∂p＝2p＋5
∂R／∂q＝－1
∂S／∂p＝q
∂S／∂q＝p＋5
となるでしょう。

以上より、初めの方で掲げた連立方程式の各パラメータが求まりまして
(2p＋5)Δp－Δq＝－p^2－5p＋q－11
qΔp＋(p＋5)Δq＝－pq－5q＋15
となります。後はこの式からΔp, Δqを求めるだけですね。

p, qがわかればΔp, Δqが得られるので、得られたΔp, Δqをp, qに足して、新しくなったp, qからまたΔp, Δqを得て……というのを延々くりかえせば、いずれRとSは0に収束し、2次式x^2＋px＋qで割り切れるときがくることになります。

実験

初期値はp＝q＝0としてみます。

実際に、
(2p＋5)Δp－Δq＝－p^2－5p＋q－11
qΔp＋(p＋5)Δq＝－pq－5q＋15
の連立方程式を用いてΔp, Δqを随時求めていきましょう。

1回目
5Δp－Δq＝－11
5Δq＝15
(Δp, Δq)＝(－1.6, 3)が得られますので、
(p, q)＝(－1.6, 3)に更新されます。
このとき、(R, S)＝(2.56, －4.8)なので、まだ0には程遠いです。

2回目
－1.8Δp－Δq＝－2.56
3Δp＋3.4Δq＝4.8
(Δp, Δq)≒(－0.42807, 1.7895)が得られますので、
(p, q)＝(－2.0281, 4.7895)に更新されます。
このとき、(R, S)≒(0.18319, －0.76608)なので、0に近づきました。

3回目
0.9438Δp－Δq＝－0.18319
4.7895Δp＋2.9719Δq＝0.76608
(Δp, Δq)≒(0.029188, 0.21074)が得られますので、
(p, q)＝(－1.9989, 5.0002)に更新されます。
このとき、(R, S)≒(0.000901, 0.0061)です。もういいですよね。

p≒－2, q≒5ということがわかりました。
つまり、x^2－2x＋5で割り切れるだろうということになります。
このときb＝－p－5＝－3ですので、

x^3－5x^2＋11x－15
＝(x^2－2x＋5)(x－3)
と因数分解できるということになります。
x^2－2x＋5＝0と置くとこれはx＝1±2iと求まります。
次に残った(x－3)の部分が3次以上だったら同様にして2次を取り出していくのですが、もうこれは1次式なのでこれでおしまいです。x＝3が求まります。

従って、x＝1±2i, 3であると求まります。

手計算だと長いですけれど、初期値さえうまく決めれば
手順どおりやれば全部求めることが出来るのでとっても便利だと思います。

参考URL

http://ysserve.wakasato.jp/Lecture/NumericalCalculation1/chap6.pdf