6. 代数方程式

6.1 概要

代数方程式とは多項式のゼロ点を求める問題です.
\[
p(x) = a_0x^n + a_1x^{n-1} + \dots + a_{n-1}x + a_n = 0
\] よく知られているように 2, 3, 4 次方程式については解の公式で解くことができます. 5 次以上(高次代数方程式)は一般的には数値解法(反復法)で解く必要があります. なお, 3, 4 次であっても数値解法で解いた方がよい場合もあります.

6.2 高次代数方程式の解法

高次代数方程式については多くの解法が提案されてきました. 例えば, ベアストウ法, グレッフェ法, ベルヌーイ法, ジェンキンス・トラウブ法などですが, 決定版というような解法はなかったようです.

特定の解を１つだけ求めればよい場合には 1 変数の非線形方程式としてニュートン法を使うことができます. ただし, よい初期値がわかっている必要があります.

すべての解を求めたい場合には次のような解法が使われています.

6.2.1 連立法

すべての解を同時に求める反復法です.

\(k\) 番目の反復の近似解を上付きの \((k)\) で表すことにすると, それぞれの解をニュートン法で求めるときの反復式は次のようになります.
\[
x_i^{(k+1)} = x_i^{(k)} – p(x_i^{(k)})/p'(x_i^{(k)})
\] ここで, \(p'(x_i^{(k)}) \simeq \prod_{j=1}^{n, j \neq i}(x_i^{(k)} – x_j^{(k)})\) で近似すると次のようになります.
\[
x_i^{(k+1)} = x_i^{(k)} – p(x_i^{(k)})/\prod_{j=1}^{n, j \neq i}(x_i^{(k)} – x_j^{(k)})
\] これは連立法 (2 次法) と呼ばれ, 重根がなければ 2 次収束します. デュラン・ケルナー法と呼ばれることもあります.

近似の仕方や初期値の与え方にはバリエーションがあります. 少し計算量が増えるものの, 重根がなければ 3 次収束する近似式もあり 3 次法と呼ぶことがあります. アバースが提案した初期値を採用したものはデュラン・ケルナー・アバース (DKA) 法と呼ばれます.

6.2.2 随伴行列法

固有値問題において「固有値 \(\lambda\) に関する n 次代数方程式を特性方程式と呼び, これを解けば固有値が求められる」というような解説がされることがあります. 実際には逆で, 代数方程式を直接解く方が難しいため, 解くべき代数方程式に対応する固有値問題を組み立てて (随伴行列という) それを QR 法などの固有値問題の解法で解くことにより間接的に代数方程式の解を求める方法があります.

6.3 XLPack を使った代数方程式の解き方

VBA サブルーチン Rpzero2 は実数係数の代数方程式を連立法(2 次法)で解くプログラムです. 複素解があればそれも求めることができます. ワークシート関数としては WRpzero2 が使用できます.

例題

次の 5 次方程式を解く.
\[
x^5 + 2x^3 + 2x^2 – 15x + 10 = 0
\] これは因数分解できて,
\[
(x – 1)^2(x + 2)(x^2 + 5) = 0
\] となるので, 解は \(1(重根), -2, \pm\sqrt{5}i\) となることがわかっています.

6.3.1 ワークシート関数を使用した解き方

ワークシートの適当な場所に代数方程式の係数 \(a_i\) を入力します(オレンジ色のセル). 次に, 解 \((x(r), x(i))\) および 誤差限界 \(s(i)\) を入れるためのセルを \((n + 1) \times 3\) セル選択してワークシート関数 WRpzero2 を入力します(緑色のセル).

WRpzero2 の必要なパラメータは N および A です. N は方程式の次数(この例では 5), A は方程式の係数のセル範囲です. ここで, 方程式を \(a_0 x^n + a_1x^{n-1} + \dots + a_{>n-1}x + a_n = 0\) とすると, 係数は実数 \(a_0, a_1, \dots, a_n\) で表されます.

入力が終了したら Ctrl + Shift + Enter を押します.

これで, 5 つの解が正しく求められました. なお, 一番右の列 s(i) は誤差限界, 一番下の行は左からリターンコード, 収束に要した反復回数を表します.

重根 (1) については精度が半分程度に落ちているのがわかります. 代数方程式の数値計算においては, \(m\) 重根を求める場合に解の有効桁数がおおよそ \(1/m\) になるので注意が必要です.

6.3.2 VBAプログラムを使用した解き方

上と同じ例題を VBA プログラムにより解きます. Rpzero2 を使ったプログラム例を示します.

Sub Start()
    Const NMax = 10
    Dim N As Long
    Dim A(NMax) As Double, XR(NMax) As Double, XI(NMax) As Double, S(NMax) As Double
    Dim IFlag As Long, Info As Long, Iter As Long, I As Long
    '--- Input data
    N = 5
    For I = 0 To N
        A(I) = Cells(5 + I, 1)
    Next
    '--- Compute roots of equation
    IFlag = 0
    Call Rpzero2(N, A(), XR(), XI(), IFlag, S(), Info, Iter)
    Cells(10, 2) = Info
    Cells(10, 3) = Iter
    '--- Output roots
    For I = 0 To N - 1
        Cells(5 + I, 2) = XR(I)
        Cells(5 + I, 3) = XI(I)
        Cells(5 + I, 4) = S(I)
    Next
End Sub

所定の位置(下のオレンジ色のセル)に係数データを入力し, マクロ Start を実行すると次の結果が得られます. ワークシート関数を使用したときと異なり, 係数データを入力しただけでは結果が得られず, VBA プログラムを実行してやる必要があります.