§６．代数方程式
［第１１回］
6.1 ベアストウ法
n ≥ 3 の代数方程式の数値解を求める方法の一つにベアストウ法がある．
f (z ) = z n + a1z n −1 + ! + an −1z + an = 0
(6.1)
2
この式を 2 次式: z + pz + q で割ると一般に，
f (z ) = (z 2 + pz + q )("###############
z n −2 + b1z n −3 + #!
+ bn −3z + bn −2 ) +
+S
"##$##%
$
# ###############
##% Rz
商
(6.2)
余り
もし余り R = S = 0 ならば f (z ) = 0 の２根が求まる．すなわち
−p ± p 2 − 4q
2
z=
(6.3)
一般に p, q の値は R, S の値に影響を与えるので R(p, q ), S (p, q ) と考える．そこで
R(p, q ) = 0
S (p, q ) = 0
(6.4)
を満たすような p, q の値を求めるアルゴリズムを作る．上式を満たす p, q は未知であるから始めに
試行値 pi , qi を与える． pi , qi の修正値を ∆p, ∆q とすると，
R(pi + ∆p, qi + ∆q ) ≅ R(pi , qi ) + Rp (pi , qi )∆p + Rq (pi , qi )∆q
S (pi + ∆p, qi + ∆q ) ≅ S (pi , qi ) + S p (pi , qi )∆p + Sq (pi , qi )∆q
ただし
∂R
,
∂p
∂S
Sp =
,
∂p
∂R
∂q
∂S
Sq =
∂q
Rp =
Rq =
ここで
R(pi + ∆p, qi + ∆q ) = 0
S (pi + ∆p, qi + ∆q ) = 0
1
(6.5)
となる ∆p, ∆q を決定する．
∆pRp (pi , qi ) + ∆qRq (pi , qi ) = −R(pi , qi )
∆pS p (pi , qi ) + ∆qSq (pi , qi ) = −S (pi , qi )
(6.6)
上式において
Rp (pi , qi ),
Rq (pi , qi ),
R(pi , qi )
S p (pi , qi ),
Sq (pi , qi ),
S (pi , qi )
を求めることができれば ∆p, ∆q は(6.6)式より求められる．
∆p = pi +1 − pi
∆q = qi +1 − qi
(6.7)
とおけば試行値の修正値 pi +1, qi +1 を得る．すなわち
pi +1 = ∆p + pi
qi +1 = ∆q + qi
(6.8)
i ← i + 1 として反復計算を行えば pi +1, qi +1 は真の値 p, q に収束する．
R(pi , qi ), S (pi , qi ) の計算：
(6.2)式を展開し(6.1)式の係数と比較すると（ p → pi , q → qi とおく）．
a 0 = 1 = b0
a1 = b1 + pi
a2 = b2 + pib1 + qi
&
ak = bk + pibk −1 + qibk −2
&
an −1 = R + pibn −2 + qibn−3
an = S + qibn −2
上式を b に関する式に書き換えると
2
(6.9)
b0 = 1 = a 0
b = a − p
i
1
1
b2 = a2 − pib1 − qi
&
bk = ak − pibk −1 − qibk −2
&
R = an −1 − pibn −2 − qibn −3
S = an − qibn −2
(6.10)
ここで漸化式：
bk = ak − pibk −1 − qibk −2
(b1 = a1 − pi , b0 = 1)
(6.11)
を使って k = 2 ∼ n まで繰り返すと，最後の 2 つの式は
bn −1 = an −1 − pibn −2 − qibn −3
bn = an − pibn −1 − qibn −2
(6.12)
R = an −1 − pibn −2 − qibn −3 = bn −1
S = an − qibn −2 = bn + pibn −1
(6.13)
となる．これを使うと
となり， R, S は bn と bn−1 によって計算できる．
Rp (pi , qi ), Rq (pi , qi ), S p (pi , qi ), Sq (pi , qi ) の計算
(6.13)式より
∂bn −1
Rp = ∂p
i
∂b
Rq = n −1
∂qi
∂bn
+ pi
S p =
∂pi
Sq = ∂bn + pi
∂q i
これらの偏微分係数に関する漸化式を導く．
3
∂bn −1
+ bn −1
∂pi
∂bn −1
∂qi
(6.14)
∂b
ck = k
∂pi
∂bk
dk =
∂q i
(6.15)
を定義し，(6.10)式を pi , qi でそれぞれ偏微分すれば（ k = n まで），
c1 = −1
c2 = −b1 − pic1
&
c = −b − p c − q c
k −1
i k −1
i k −2
k
&
cn −1 = −bn −2 − picn −2 − qicn −3
c = −b − p c − q c
n −1
i n −1
i n −2
n
(6.16)
d1 = 0
d = −1
2
&
dk = −bk −2 − pidk −1 − qidk −2
&
dn −1 = −bn −3 − pidn −2 − qidn −3
dn = −bn −2 − pidn −1 − qidn −2
(6.17)
Rp = cn −1
R = dn −1
q
S p = cn + picn −1 + bn −1
S = dn + pidn −1
q
(6.18)
∆pcn −1 + ∆qdn −1 = −bn −1
∆p(cn + picn −1 + bn −1 ) + ∆q(dn + pidn −1 ) = −bn − pibn −1
(6.19)
よって
(6.6)式は
(6.19)式の第 1 式に pi をかけ，第 2 式から引くと
4
∆pcn −1 + ∆qdn −1 = −bn −1
∆p(cn + bn −1 ) + ∆qdn = −bn
(6.20)
これより ∆p, ∆q について解くと
−bn −1 dn −1
∆p =
−bn
dn
cn −1
dn −1
cn + bn −1
dn
∆q =
,
cn −1
−bn −1
cn + bn −1
−bn
cn −1
dn −1
cn + bn −1
dn
(6.21)
この ∆p, ∆q を(6.8)式に代入すると pi , qi の修正値 pi +1, qi +1 が求められる．
以上計算アルゴリズムを整理すると，計算のステップは以下のようになる．
（１） i = 0 とし，試行値として例えば pi = 0, qi = 0 を与える．
(2)
bk = ak − pibk −1 − qibk −2 (b1 = a1 − pi , b0 = 1) を用いて k = 2 ∼ n まで繰り返し，
b2 , ! , bn −1, bn を求める．
(3)
ck = −bk −1 − pick −1 − qick −2 (c2 = −b1 − pic1, c1 = −1) を用いて k = 3 ∼ n まで繰
り返し， c3 , ! , cn −1, cn を求める．
(4)
dk = −bk −2 − pidk −1 − qidk −2 (d2 = −1, d1 = 0) を用いて k = 3 ∼ n まで繰り返し，
d3 , ! , dn −1, dn を求める．
−bn −1 dn −1
(5)
∆p =
−bn
dn
cn −1
dn −1
cn + bn −1
dn
pi +1 = ∆p + pi ,
,
∆q =
−bn −1
cn + bn −1
−bn
cn −1
dn −1
cn + bn −1
dn
qi +1 = ∆q + qi
を計算する．
(6)
cn −1
判定条件：
5
R = bn −1
S = bn + pi +1bn −1
R > ε, S > ε
−5
（ ε は例えば 10
∼ 10−10 ）ならば i ← i + 1 としてステップ(2)へ
もどる．
R ≤ ε, S ≤ ε ならば f (z ) は
f (z ) ≅ (z 2 + pi z + qi )(z n −2 + b1z n −3 + ! + bn −3z + bn −2 ) = 0
2
とみなして z + pi z + qi = 0 より 2 根を求める．
(7) n ← n − 2 とし， n = 1 ならば z + b1 = 0 より 1 根を求め，終了．
n = 2 ならば z 2 + b1z + b2 = 0 より 2 根をもとめ．終了．
n ≥ 3 ならば a1 ← b1, a2 ← b2 , ! , an ← bn としてステップ(1)へ
もどる．
(例 1)
f (z ) = (x − 2 + i )(x − 2 − i )(x − 3)(x − 4)(x + 5)
= x 5 − 6x 4 − 10x 3 + 142x 2 − 355x + 300 = 0
この 5 次方程式の根を求めよ．
（解）
上述のアルゴリズムに基づくＣ言語のプログラム例は以下のようである．
#include<math.h>
#include<stdio.h>
int main(void)
{
int i,n,k;
long double p,q,a[6],b[6],c[6],d[6],dp,dq,R,S,e,eps,
re,im,x1,x2,delta;
n=5; a[1]=-6.0; a[2]=-10.0; a[3]=142.0;
a[4]=-355.0; a[5]=300.0; eps=1.0e-10;
6
step1:
i=0; p=0; q=0;
step2:
b[0]=1.0; b[1]=a[1]-p;
for(k=2; k<=n; k++)
b[k]=a[k]-p*b[k-1]-q*b[k-2];
/*step3:*/ c[1]=-1.0; c[2]=-b[1]-p*c[1];
for(k=3; k<=n; k++)
c[k]=-b[k-1]-p*c[k-1]-q*c[k-2];
/*step4:*/ d[1]=0; d[2]=-1.0;
for(k=3; k<=n; k++)
d[k]=-b[k-2]-p*d[k-1]-q*d[k-2];
/*step5:*/ delta=c[n-1]*d[n]-d[n-1]*(c[n]+b[n-1]);
dp=(-b[n-1]*d[n]+b[n]*d[n-1])/delta;
dq=(-c[n-1]*b[n]+b[n-1]*(c[n]+b[n-1]))/delta;
p=dp+p; q=dq+q;
/*step6:*/ R=b[n-1]; S=b[n]+p*b[n-1]; e=p*p-4.0*q;
if(R*R>eps && S*S>eps){i=i+1; goto step2;}
else if(e<0){re=-0.5*p; im=0.5*sqrt(-e);
printf("%f
i %f¥n",re,im);
printf("%f
-i %f¥n",re,im);
printf("R=%e S= %e i=%d¥n", R,S,i);}
else{x1=0.5*(-p+sqrt(e)); x2=0.5*(-p-sqrt(e));
printf("%f¥n",x1); printf("%f¥n",x2);
printf("R=%e S= %e i=%d¥n", R,S,i);};
/*step7:*/ n=n-2;
if(n==1){printf("%f¥n",-b[1]); goto owari;};
if(n==2){e=b[1]*b[1]-4.0*b[2];
if(e<0){re=-0.5*p; im=0.5*sqrt(-e);
printf("%f
i %f¥n",re,im);
printf("%f -i %f¥n",re,im);
printf("R=%e S= %e i=%d¥n",R,S,i);
goto owari;}
7
else{x1=0.5*(-b[1]+sqrt(e));
x2=0.5*(-b[1]-sqrt(e));
printf("%f¥n",x1); printf("%f¥n",x2);
printf("R=%e S= %e i=%d¥n",R,S,i);
goto owari;}
};
if(n>=3){for(k=1; k<=n; k++) a[k]=b[k]; goto step1;};
owari:
return 0;
}
実行結果：
2.000000
i 1.000000
2.000000 -i 1.000000
R=-5.860450e-006 S= 8.312085e-006 i=6
4.000000
-5.000000
R=6.405685e-010 S= 3.205671e-009 i=7
3.000001
Press any key to continue
（参考）Ｃ言語の処理系は Microsoft Visual C++ 6.0 を用いた，Ｃを起動する方法および
コンパイルの方法等は下記の書物に出ている．
“やさしいＣ”
高橋麻奈著
ソフトバンク（株）
8
© Copyright 2024 ExpyDoc
