§6.代数方程式
[第11回]6.1 ベアストウ法
3
n ≥
の代数方程式の数値解を求める方法の一つにベアストウ法がある. (6.1) 1 1 1( )
n n0
n nf z
z
a z
−a z
a
−=
+
+
!
+
+
=
この式を 2 次式:z
2+
pz
+ q
で割ると一般に, 2 2 3 1 3 2( )
(
)(
n n)
n nf z
z
pz
q z
−b z
−b z
b
Rz
S
− −=
+
+
"################$
+
+
!
# ################%#
+
+
+
"##$##%
+
余り 商 (6.2) もし余りR
=
S
=
0
ならばf z =
( )
0
の2根が求まる.すなわち 24
2
p
p
z
=
− ±
−
q
i iq
(6.3) 一般にp
, q
の値はR
, S
の値に影響を与えるのでR p
( , ), ( , )
q S p q
と考える.そこで (6.4)( , )
0
( , )
0
R p q
S p q
=
=
を満たすような の値を求めるアルゴリズムを作る.上式を満たす は未知であるから始めに 試行値 を与える. の修正値を∆
とすると,,
p q
p q
,
,
ip q
p q
i,
p
,
∆
(6.5)(
,
)
( , )
( , )
( , )
(
,
)
( , )
( , )
( , )
i i i i p i i q i i i i i i p i i q i iR p
p q
q
R p q
R p q
p
R p q
q
S p
p q
q
S p q
S p q
p
S p q
q
+ ∆
+ ∆ ≅
+
∆ +
∆
+ ∆
+ ∆ ≅
+
∆ +
∆
ただし,
,
p q p qR
R
R
R
p
q
S
S
S
S
p
q
∂
∂
=
=
∂
∂
∂
∂
=
=
∂
∂
ここで(
,
)
(
,
)
i i i iR p
p q
q
S p
p q
q
+ ∆
+ ∆ =
+ ∆
+ ∆ =
0
0
となる
∆
p ∆q
,
を決定する.q
1p
q
1 (6.6)( , )
( , )
( , )
( , )
( , )
( , )
p i i q i i i i p i i q i i i ipR p q
qR p q
R p q
pS p q
qS p q
S p q
∆
+ ∆
= −
∆
+ ∆
= −
上式において( , ),
( , ),
( , )
( , ),
( , ),
( , )
p i i q i i i i p i i q i i i iR p q
R p q
R p q
S p q
S p q
S p q
を求めることができれば∆
p ∆
,
は(6.6)式より求められる. (6.7) 1 1 i i i ip
p
p
q
q
q
+ +∆ =
−
∆ =
−
とおけば試行値の修正値p
i+1,
q
i+ を得る.すなわち (6.8) 1 1 i i i ip
p
q
q
+ +
= ∆ +
= ∆ +
1
i
← +
i
として反復計算を行えばp
i+1,
q
i+ は真の値p q
,
に収束する.( , ), ( , )
i i i iR p q
S p q
の計算: (6.2)式を展開し(6.1)式の係数と比較すると(p
→
p q
i,
→
q
iとおく). 3 i n− (6.9) 0 0 1 1 2 2 1 1 2 1 2 21
i i i k k i k i k n i n n i na
b
a
b
p
a
b
p b
q
a
b
p b
q b
a
R
p b
q b
a
S
q b
− − − − − = =
= +
= +
+
=
+
+
=
+
+
= +
&
&
上式をb
に関する式に書き換えると(6.10) 0 0 1 1 2 2 1 1 2 1 2 2
1
i i i k k i k i k n i n i n n i nb
a
b
a
p
b
a
p b
q
b
a
p b
q b
R
a
p b
q b
S
a
q b
− − − − − = =
= −
= − −
=
−
−
=
−
−
= −
&
&
3 −1)
n
3 − 1 − ここで漸化式: (6.11) 1 2(
1 1,
0 k k i k i k ib
=
a
−
p b
−−
q b
−b
=
a
−
p b
=
を使ってk
= ∼
2
まで繰り返すと,最後の 2 つの式は (6.12) 1 1 2 1 2 n n i n i n n n i n i nb
a
p b
q b
b
a
p b
q b
− − − − −
=
−
−
= −
−
となる.これを使うと (6.13) 1 2 3 2 1 n i n i n n n i n n i nR
a
p b
q b
b
S
a
q b
b
p b
− − − − − =
−
−
=
= −
=
+
となり,R
, S
はb
nとb
n−1によって計算できる.( , ),
( , ),
( , ),
( , )
p i i q i i p i i q i iR p q
R p q
S p q
S p q
の計算 (6.13)式より 1 1 1 1 1 n p i n q i n n p i i i n n q i i ib
R
p
b
R
q
b
b
S
p
p
p
b
b
S
p
q
q
− − − − −
∂
=
∂
∂
=
∂
∂
∂
=
+
+
∂
∂
∂
∂
=
+
∂
∂
nb
(6.14) これらの偏微分係数に関する漸化式を導く.k k i k k i
b
c
p
b
d
q
∂
=
∂
∂
=
∂
(6.15) を定義し,(6.10)式をp
i,
q
iでそれぞれ偏微分すれば(k
まで), 3 2 − 3 2 − 1 −n
=
(6.16) 1 2 1 1 1 1 2 1 2 2 1 11
i k k i k i k n n i n i n n n i n i nc
c
b
p c
c
b
p c
q c
c
b
p c
q c
c
b
p c
q c
− − − − − − − − − = −
= − −
= − −
−
= −
−
−
= −
−
−
&
&
(6.17) 1 2 2 1 2 1 3 2 2 10
1
k k i k i k n n i n i n n n i n i nd
d
d
b
p d
q d
d
b
p d
q d
d
b
p d
q d
− − − − − − − − − =
= −
= − −
−
= −
−
−
= −
−
−
&
&
よって (6.18) 1 1 1 1 1 p n q n p n i n n q n i nR
c
R
d
S
c
p c
b
S
d
p d
− − − − −
=
=
= +
+
= +
(6.6)式は (6.19) 1 1 1 1 1 1(
)
(
)
n n n n i n n n i n n i npc
qd
b
p c
p c
b
q d
p d
b
p b
− − − − − −∆
+ ∆
= −
∆
+
+
+ ∆
+
= − −
(6.19)式の第 1 式にp
iをかけ,第 2 式から引くと(6.20) 1 1 1
(
)
n n n n n npc
qd
b
p c
b
qd
b
− − − −∆
+ ∆
= −
∆
+
+ ∆
= −
1 nq
これより∆
p ∆
,
について解くと 1 1 1 1 1 1 1 1 1 1,
n n n n n n n n n n n n n n n n n n nb
d
c
b
b
d
c
b
b
p
q
c
d
c
d
c
b
d
c
b
d
− − − − − − − − − − −−
−
−
+
−
∆ =
∆ =
+
+
1q
i 11)
n
(6.21) この∆ ∆
p
,
を(6.8)式に代入するとp q
i,
の修正値p
i+1,
q
i+ が求められる. 以上計算アルゴリズムを整理すると,計算のステップは以下のようになる. (1)i =
0
とし,試行値として例えばp
i=
0,
q
i=
0
を与える. (2)b
k=
a
k−
p b
i k−1−
q b
i k−2(
b
1=
a
1−
p b
i,
0=
を用いてk
まで繰り返し, 2,
,
n 1,
nb
!
b
−b
2
= ∼
を求める. (3) を用いてk
まで繰 り返し,c
1 1 2(
2 1 1,
1 k k i k i k ic
= −
b
−−
p c
−−
q c
−c
= − −
b
p c c
= −
3,
!
,
c
n−1,
c
n1)
= ∼
3
n
を求める. (4)d
k= −
b
k−2−
p d
i k−1−
q d
i k−2(
d
2= −
1,
d
1=
を用いてk
まで繰り返し, 3,
,
d !
d
n−1,
d
n0)
= ∼
3
n
を求める. (5) 1 1 1 1 1 1 1 1 1,
n n n n n n n n n n n n n n n n n n nb
d
c
b
b
d
c
b
b
p
q
c
d
c
d
c
b
d
c
b
d
− − − − − − − − − − −−
−
−
+
−
∆ =
∆ =
+
+
1 1 iq
p
i+1= ∆ +
p
p
i,
q
i+1= ∆ +
q
を計算する. (6) 判定条件:
1 1 1 n n i nR
b
S
b
p b
− + − =
=
+
R
>
ε
,
S
>
ε
( は例えば10
)ならばi
i
としてステップ(2)へ もどる.ε
510
−∼
−10← +
1
R
≤
ε
,
S
≤
ε
2 nb
ならばf z
( )
は 2 2 3 1 3( )
(
)(
n n)
0
i i n nf z
z
p z
q z
−b z
−b z
b
− −≅
+
+
+
+
!
+
+
=
とみなしてz
2p
0
より2 根を求める. iz
q
i+
+
=
(7)n
← −
n
2
とし,n
=
1
ならばz
+
b
1=
0
より1 根を求め,終了.2
n =
ならば 2 より2 根をもとめ.終了. 1 20
z
+
b z
+
b
=
3
n ≥
ならばa
1←
b
1,
a
2←
b
2,
!
,
a
n←
としてステップ(1)へ もどる. (例 1) 5 4 3 2( )
(
2
)(
2
)(
3)(
4)(
5)
6
10
142
355
300
0
f z
x
i x
i x
x
x
x
x
x
x
x
=
− +
− −
−
−
+
=
−
−
+
−
+
=
この 5 次方程式の根を求めよ. (解) 上述のアルゴリズムに基づくC言語のプログラム例は以下のようである. #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;
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;}
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
(参考)C言語の処理系はMicrosoft Visual C++ 6.0 を用いた,Cを起動する方法および コンパイルの方法等は下記の書物に出ている.