1. No category

情報処理演習（１０）
ｎ次方程式の複素数解
今日のテーマ
・複素数を使った計算。
・３次方程式の解の公式（カルダノの公式）
・n 次方程式の解法（DKA 法）
３次方程式の解の公式（カルダノの公式）
３次方程式 x3 + a1x 2 + a2 x + a3 = 0 の解は次の方法で求めることができる。
１． x = y −
a1
とおくことにより、２次の項を消去すると、
3
y 3 + ay + b = 0 となる。ここで、a と b は
a=−
a12
+ a2
3
b=
2a13 a1a2
−
+ a3
27
3
である。
２．次に、u と v を未知数とする方程式
u 3 + v 3 + b = 0 と 3uv + a = 0
を考える。もし、a≠０ならば、u と v は次の式で求められる。
⎛ b
b
a ⎞⎟
u = ⎜− +
+
⎜ 2
4 27 ⎟⎠
⎝
2
3
1
3
v=−
a
3u
√の中が負のときは、虚数になるため、u を計算するには、複素数の３乗根を計算しなく
てはいけない。Fortran では、変数のデータ型を複素数(complex 型）にしておけば、複素
数の３乗根（＝1/3 乗）が計算できる。
３． １の３乗根のうち、虚数のものをωとする。
⎛ 3
⎞
−1+ i 3
⎜
ω=
= 0.86602540...⎟
（i は虚数単位）
⎜
⎟
2
⎝ 2
⎠
４． このとき、 y 3 + ay + b = 0 の解は、次の３つである。
y2 = ω u + ω 2v
y1 = u + v
５．
y3 = ω 2 u + ω v
a1
の関係を使って、ｘの３つの値を出す。
3
a
a
x2 = y 2 − 1
x3 = y3 − 1
3
3
上の３つの解から、 x = y −
x1 = y1 −
a1
3
これが元の３次方程式の解である。以上の手順をカルダノの公式という。この方法で解を求
めるときには、途中の計算に必ず複素数が現れる。（このことが、数学の歴史上、複素数とい
うものを取り入れるきっかけになった。）
問題１０−１
ラムを作れ。
解答
C
カルダノの公式を用いて、３次方程式 x3 + a1x 2 + a2 x + a3 = 0 の解を求めるプログ
上の手順をプログラムにする。
real a1,a2,a3
complex a,b,u,v,y1,y2,y3,omega,x1,x2,x3
↑変数はすべて複素数型にしてある。
write(*,*) '3 次方程式 x**3+a1*x**2+a2*x+a3=0 の解を求めます'
write(*,*) 'a1 を入力してください'
read(*,*) a1
write(*,*) 'a2 を入力してください'
read(*,*) a2
write(*,*) 'a3 を入力してください'
read(*,*) a3
omega=(-0.5 , 0.86602540)
↑複素数型の定数はこのように(実数部,虚数部)の形で表す。
この omega は 1 の３乗根（虚数）である。
a= a の計算式
<--□の部分は、どんな式が入るか？
b= b の計算式
u=(-b/2+((b**2)/4+(a**3)/27)**0.5)**(1.0/3.0)
v= v の計算式
y1= y1 の計算式
y2= y2 の計算式
y3= y3 の計算式
x1= x1 の計算式
x2= x2 の計算式
x3= x3 の計算式
write(*,*) '答えは次の通りです'
write(*,*) 'x1=',x1
write(*,*) 'x2=',x2
write(*,*) 'x3=',x3
write(*,*) '検算します。'
write(*,*)'x1 を方程式に代入すると', x1 を代入した式
write(*,*)'x2 を方程式に代入すると', x2 を代入した式
write(*,*)'x3 を方程式に代入すると', x3 を代入した式
end
C
C
プログラムを実行した結果の複素数は、
(実数部,虚数部)
の形で表示される。たとえば、3+2i という複素数は、(3.0, 2.0)のように表示される。
問
上のプログラムを使って、次の３次方程式の解を求めよ。
（２）
x 3 − 6 x 2 + 11x − 6 = 0
x3 − 2 x 2 + x = 0
（３）
x3 + x + 1 = 0
（１）
次に、一般の n 次方程式の解を、反復計算によって求めよう。代表的な方法として、DKA
法を紹介する。
方程式 x n + an −1 x n −1 + L + a2 x 2 + a1 x + a0 = 0 は、複素数の範囲で n 個の解をもつ。いま、解に
近い値を適当にきめて、x1, x2, ..., xn とする。次のような式を計算する。
xk n −1 + an −1 x n −1 + L + a2 xk 2 + a1 xk + a0
（k=1,2,..,n)
xk′ = xk −
( xk − x1 )( xk − x2 ) L ( xk − xk −1 )( xk − xk +1 ) L ( xk − xn )
上の式の分母は、( xk − x j ) の積を j＝1,2,...,n（ただし j≠k）にわたってとるということである。
左辺の x'1, x'2, ..., x'n を新しい x1, x2, ..., xn の値として、再び上の式を計算することを繰り返し
ていくと、x1, x2, ..., xn はしだいに方程式の解に近づいていく。x1, x2, ..., xn の初期値は任意に
とってよい。ここでは、つぎのようにとることにする。
1⎞π
⎛
θ k = ⎜ 2(k − 1) + ⎟
とし、r を定数として
2⎠ n
⎝
a1
+ r exp(iθ k ) とする。
（k=1,2,..,n)
n
ここで、i は虚数単位であり、 exp(iθ ) = cos θ + i sin θ である。r の決め方が難しいが、ここでは
初期値を
xk = −
深く考えずに、
「適当に」とることにしよう。DKA 法の計算はすべて複素数で行う必要がある
が、FORTRAN には複素数のデータ型（COMPLEX）があるのでプログラムが簡単である。
問題
る。
ｎ次代数方程式 x n + an −1 x n −1 + L + a2 x 2 + a1 x + a0 = 0 の解を複素数の範囲ですべて求め
解答
c 代数方程式の解法（DKA 法）
c n 次方程式 a(0)+a(1)*x+a(2)*x**2+...+a(n)*x**n=0
c の n 個の複素数解 x(1),x(2),...,x(n)を求める
parameter(n=3)
<--この文の意味は、下に説明がある。
real a(0:n)
complex c,f,p,x(1:n),xb(1:n)
write(*,*) n,'次方程式 a(0)+a(1)*x+a(2)*x**2+...+a(n)*x**n=0'
write(*,*)'の解を複素数の範囲で計算します'
do k=n,0,-1
write(*,*) k,'次項の係数 a(',k,')を入力してください'
read(*,*) a(k)
end do
c 最高次の係数が１になるように，各項の係数を設定しなおす
do k=0,n
a(k)=a(k)/a(n)
end do
c x(1),x(2),...x(n)の初期値を設定する
r=10
pi=3.14159265358979
do k=1,n
c= (0.0,1.0)*pi*(2.0*(k-1)+0.5)/n
x(k)=-a(1)/n+r*exp(c)
end do
write(*,10) (' (x(',k,')の実数部, 虚数部 ) ',k=1,n)
10
format(3(A4,I1,A19))
write(*,20) (x(k),k=1,n)
c DKA 法を実行する
do m=1,50
c x(k)の計算
do k=1,n
c 分母の計算
p=1
do j=1,n
if (k .ne. j) then
p=p*(x(k)-x(j))
end if
end do
c 分子の計算
f=1
do j=n-1,0,-1
f=f*x(k)+a(j)
end do
c 新しい x(k)の値を保存しておく
xb(k)=x(k)-f/p
end do
c x(k)の値を更新する
e=0
do k=1,n
e=e+abs(xb(k)-x(k))
x(k)=xb(k)
end do
write(*,20) (x(k),k=1,3)
20
format (3('(',f8.4,',',f8.4,') '))
if (e .LE. 1.0e-6) stop
end do
end
練習
（１）
上のプログラムを用いて、次の方程式の解を求めよ
x3 + x 2 + x + 1 = 0
（２）
（３）
x3 + 1 = 0
− 4 − 2 x + x3 = 0
・parameter(n=3)の意味
この文は、特定の数値（この場合３）を表す文字（この場合ｎ）を作る働きをする。このプ
ログラムの中で、３は方程式の次数を表しているため、配列の大きさをきめるときや、DO 文
の繰り返し数など、いろいろなところに表れる。そういうときには、３という数をそのままか
くよりもｎという文字で表しておくほうが便利なのである。たとえば、上のプログラムを４次
方程式を解くように、作りなおすことを考えてみるとよい。もしｎを３と書いてしまうと、プ
ログラムのいろんな部分を変更しなければならなくなる。そのようなとき、parameter 文で
n=3 としておくと、ここだけを変更すればすむ。
問題１０−３
上のプログラムを４次方程式を解くように、書き直してみよ。そして、方程式
x + x + x + x + 1 = 0 の解を求めよ。
4
3
ヒント：
2
parameter 文、２つの format 文の計３箇所を変更。