§6

Loading.... (view fulltext now)

Loading....

Loading....

Loading....

Loading....

全文

(1)

§6.代数方程式

[第11回]

6.1 ベアストウ法

3

n ≥

の代数方程式の数値解を求める方法の一つにベアストウ法がある. (6.1) 1 1 1

( )

n n

0

n n

f z

z

a z

a z

a

=

+

+

!

+

+

=

この式を 2 次式:

z

2

+

pz

+ q

で割ると一般に, 2 2 3 1 3 2

( )

(

)(

n n

)

n n

f z

z

pz

q z

b z

b z

b

Rz

S

− −

=

+

+

"################$

+

+

!

# ################%#

+

+

+

"##$##%

+

余り 商 (6.2) もし余り

R

=

S

=

0

ならば

f z =

( )

0

の2根が求まる.すなわち 2

4

2

p

p

z

=

− ±

q

i i

q

(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

,

,

i

p 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 i

R 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 q

R

R

R

R

p

q

S

S

S

S

p

q

=

=

=

=

ここで

(

,

)

(

,

)

i i i i

R p

p q

q

S p

p q

q

+ ∆

+ ∆ =





+ ∆

+ ∆ =



0

0

(2)

となる

p ∆q

,

を決定する.

q

1

p

q

1 (6.6)

( , )

( , )

( , )

( , )

( , )

( , )

p i i q i i i i p i i q i i i i

pR 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 i

R 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 i

p

p

p

q

q

q

+ +

∆ =



∆ =



とおけば試行値の修正値

p

i+1

,

q

i+ を得る.すなわち (6.8) 1 1 i i i i

p

p

q

q

+ +

= ∆ +



 = ∆ +



1

i

← +

i

として反復計算を行えば

p

i+1

,

q

i+ は真の値

p q

,

に収束する.

( , ), ( , )

i i i i

R 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 2

1

i i i k k i k i k n i n n i n

a

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

に関する式に書き換えると

(3)

(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 n

b

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 i

b

=

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 n

b

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 n

R

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 i

R 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 i

b

R

p

b

R

q

b

b

S

p

p

p

b

b

S

p

q

q

− − − − −

 =





=





=

+

+





 =

+





n

b

(6.14) これらの偏微分係数に関する漸化式を導く.

(4)

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 1

1

i k k i k i k n n i n i n n n i n i n

c

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 1

0

1

k k i k i k n n i n i n n n i n i n

d

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 n

R

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 n

pc

qd

b

p c

p c

b

q d

p d

b

p b

− − − − − −

∆

+ ∆

= −



∆

+

+

+ ∆

+

= − −



(6.19)式の第 1 式に

p

iをかけ,第 2 式から引くと

(5)

(6.20) 1 1 1

(

)

n n n n n n

pc

qd

b

p c

b

qd

b

− − − −

∆

+ ∆

= −



∆

+

+ ∆

= −



1 n

q

これより

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 n

b

d

c

b

b

d

c

b

b

p

q

c

d

c

d

c

b

d

c

b

d

− − − − − − − − − − −

+

∆ =

∆ =

+

+

1

q

i 1

1)

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 k1

q b

i k2

(

b

1

=

a

1

p b

i

,

0

=

を用いて

k

まで繰り返し, 2

,

,

n 1

,

n

b

!

b

b

2

= ∼

を求める. (3) を用いて

k

まで繰 り返し,

c

1 1 2

(

2 1 1

,

1 k k i k i k i

c

= −

b

p c

q c

c

= − −

b

p c c

= −

3

,

!

,

c

n−1

,

c

n

1)

= ∼

3

n

を求める. (4)

d

k

= −

b

k2

p d

i k1

q d

i k2

(

d

2

= −

1,

d

1

=

を用いて

k

まで繰り返し, 3

,

,

d !

d

n1

,

d

n

0)

= ∼

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 n

b

d

c

b

b

d

c

b

b

p

q

c

d

c

d

c

b

d

c

b

d

− − − − − − − − − − −

+

∆ =

∆ =

+

+

1 1 i

q

p

i+1

= ∆ +

p

p

i

,

q

i+1

= ∆ +

q

を計算する. (6) 判定条件:

(6)



1 1 1 n n i n

R

b

S

b

p b

− + −

 =



 =

+



R

>

ε

,

S

>

ε

( は例えば

10

)ならば

i

i

としてステップ(2)へ もどる.

ε

5

10

−10

← +

1

R

ε

,

S

ε

2 n

b

ならば

f z

( )

は 2 2 3 1 3

( )

(

)(

n n

)

0

i i n n

f z

z

p z

q z

b z

b z

b

− −

+

+

+

+

!

+

+

=

とみなして

z

2

p

0

より2 根を求める. i

z

q

i

+

+

=

(7)

n

← −

n

2

とし,

n

=

1

ならば

z

+

b

1

=

0

より1 根を求め,終了.

2

n =

ならば 2 より2 根をもとめ.終了. 1 2

0

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;

(7)

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;}

(8)

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を起動する方法および コンパイルの方法等は下記の書物に出ている.

Updating...

参照

Updating...

関連した話題 :