§６

6. 代数方程式 [ 第 回 ] 6. ベアストウ法 3 の代数方程式の数値解を求める方法の一つにベアストウ法がある. fz () z + az +! + a z+ a 0 この式を 次式 : z + pz +q で割ると一般に, 3 fz () ( z + pz+ q)( "###############$# z + bz +! ############## + b 3z+ b ) + #%# Rz "##$##% + S もし余りR S 0 ならば fz () 0の 根が求まる. すなわち z 商 p ± p 一般に p, q の値はR, S の値に影響を与えるのでRp (,), q Spq (,) と考える. そこで 4q 余り (6.) (6.) (6.3) Rpq (,) 0 Spq (,) 0 (6.4) を満たすような pq, の値を求めるアルゴリズムを作る. 上式を満たす pq, は未知であるから始めに 試行値 p q を与える., の修正値を p とすると,, p q, q R( p + p, q + q) R( p, q) + Rp( p, q) p + Rq( p, q) q S( p + p, q + q) S( p, q) + Sp( p, q) p + Sq( p, q) q (6.5) ただし ここで R S p p R R, Rq p q S S, Sq p q Rp ( + pq, + q) 0 Sp ( + pq, + q) 0

となる p, q を決定する. prp( p, q) + qrq( p, q) R( p, q) psp( p, q) + qsq( p, q) S( p, q) (6.6) 上式において R ( p, q ), R ( p, q ), R( p, q ) p q S ( p, q ), S ( p, q ), S( p, q ) p q を求めることができれば p, q は (6.6) 式より求められる. p p+ p q q+ q (6.7) とおけば試行値の修正値 p, + q+ を得る. すなわち p+ p + p q+ q + q (6.8) + として反復計算を行えば p は真の値 pq に収束する., + q+, Rp (, q), Sp (, q) の計算 : (6.) 式を展開し (6.) 式の係数と比較すると ( p p q q とおく )., a0 b0 a b + p a b + pb + q a b + pb + qb a R + pb + qb a S + qb k k k k 3 (6.9) 上式を b に関する式に書き換えると

b0 a0 b a p b a pb q b a pb qb R a pb qb S a qb k k k k 3 (6.0) ここで漸化式 : b a pb qb ( b a p, b ) k k k k 0 (6.) を使って k まで繰り返すと, 最後の つの式は b a pb qb b a pb qb 3 (6.) となる. これを使うと R a pb qb b S a qb b + pb 3 (6.3) となり,R, S はb とb によって計算できる. Rp( p, q), Rq( p, q), Sp( p, q), Sq( p, q) の計算 (6.3) 式より b Rp p b Rq q b b Sp + p + b p p b b これらの偏微分係数に関する漸化式を導く. Sq + p q q (6.4) 3

bk ck p bk dk q を定義し,(6.0) 式を p, でそれぞれ偏微分すれば (k まで ), q (6.5) c c b pc ck bk pc k qc k c b pc qc c b pc qc 3 d 0 d dk bk pd k qd k d b pd qd d b pd qd 3 3 (6.6) (6.7) よって (6.6) 式は Rp c Rq d S c + pc + b Sq d + pd p (6.8) pc + qd b pc ( + pc + b ) + qd ( + pd ) b pb (6.9) 式の第 式に p をかけ, 第 式から引くと (6.9) 4

pc + qd b pc ( + b ) + qd b (6.0) これより p, q について解くと b d c b b d c + b b p, q c d c d c + b d c + b d (6.) この p, を (6.8) 式に代入すると, の修正値が求められる. q p q p, + q+ 以上計算アルゴリズムを整理すると, 計算のステップは以下のようになる. () 0 とし, 試行値として例えば p 0, q 0 を与える. () b a pb qb ( b a p, b ) を用いてk まで繰り返し, k k k k 0 b,!, b, b を求める. (3) c b pc qc ( c b pc, c ) を用いてk 3 まで繰 k k k k り返し,c,!, 3 c, c を求める. (4) d b pd qd ( d, d 0) を用いてk 3 まで繰り返し, k k k k d 3,!, d, dを求める. (5) b d c b b d c + b b p, q c d c d c + b d c + b d p p + p, q q + q + + を計算する. (6) 判定条件 : 5

R b S b + p b + R > ε, S > ε ( ε は例えば0 5 0 0 ) ならば + としてステップ () へ もどる. R ε, S εならば fz () は fz () ( z + pz+ q)( z + bz +! + b z+ b ) 0 3 3 + z + q 0 とみなしてz p より 根を求める. (7) とし, ならばz + b 0 より 根を求め, 終了. ならばz + bz + b より 根をもとめ. 終了. 0 3 ならばa b, a b,!, a b としてステップ () へ もどる. ( 例 ) fz ( ) ( x + )( x )( x 3)( x 4)( x+ 5) 5 4 3 x x x + x x + 6 0 4 355 300 0 この 5 次方程式の根を求めよ. ( 解 ) 上述のアルゴリズムに基づく C 言語のプログラム例は以下のようである. #clude<math.h> #clude<stdo.h> t ma(vod) { t,,k; log double p,q,a[6],b[6],c[6],d[6],dp,dq,r,s,e,eps, re,m,x,x,delta; 5; a[]-6.0; a[]-0.0; a[3]4.0; a[4]-355.0; a[5]300.0; eps.0e-0; 6

step: 0; p0; q0; step: b[0].0; b[]a[]-p; for(k; k<; k++) b[k]a[k]-p*b[k-]-q*b[k-]; /*step3:*/ c[]-.0; c[]-b[]-p*c[]; for(k3; k<; k++) c[k]-b[k-]-p*c[k-]-q*c[k-]; /*step4:*/ d[]0; d[]-.0; for(k3; k<; k++) d[k]-b[k-]-p*d[k-]-q*d[k-]; /*step5:*/ deltac[-]*d[]-d[-]*(c[]+b[-]); dp(-b[-]*d[]+b[]*d[-])/delta; dq(-c[-]*b[]+b[-]*(c[]+b[-]))/delta; pdp+p; qdq+q; /*step6:*/ Rb[-]; Sb[]+p*b[-]; ep*p-4.0*q; f(r*r>eps S*S>eps){+; goto step;} else f(e<0){re-0.5*p; m0.5*sqrt(-e); prtf("%f %f ",re,m); prtf("%f - %f ",re,m); prtf("r%e S %e %d ", R,S,);} else{x0.5*(-p+sqrt(e)); x0.5*(-p-sqrt(e)); prtf("%f ",x); prtf("%f ",x); prtf("r%e S %e %d ", R,S,);}; /*step7:*/ -; f(){prtf("%f ",-b[]); goto owar;}; f(){eb[]*b[]-4.0*b[]; f(e<0){re-0.5*p; m0.5*sqrt(-e); prtf("%f %f ",re,m); prtf("%f - %f ",re,m); prtf("r%e S %e %d ",R,S,); goto owar;} 7

} owar: else{x0.5*(-b[]+sqrt(e)); x0.5*(-b[]-sqrt(e)); prtf("%f ",x); prtf("%f ",x); prtf("r%e S %e %d ",R,S,); goto owar;} }; f(>3){for(k; k<; k++) a[k]b[k]; goto step;}; retur 0; 実行結果 :.000000.000000.000000 -.000000 R-5.860450e-006 S 8.3085e-006 6 4.000000-5.000000 R6.405685e-00 S 3.0567e-009 7 3.00000 Press ay key to cotue ( 参考 )C 言語の処理系は Mcrosoft Vsual C++ 6.0 を用いた,C を起動する方法および コンパイルの方法等は下記の書物に出ている. やさしい C 高橋麻奈著ソフトバンク ( 株 ) 8