ｽﾗｲﾄﾞﾀｲﾄﾙなし

４．３複素数解
（１）ニュートン・ラプソン法
複素数においても，ニュートン・ラプソンの方法を使うことができる。
f  zk 
zk 1  zk 
f zk 
プログラミング言語によっては，複素数が使えない。
自分で定義
複素数演算のＶＢでの定義
①複素数型のユーザ定義と複素数型定数の設定
Public Type Complex
R As Double
I As Double
End Type
Public Function setCom(R, I) As Complex
setCom.R = R
setCom.I = I
End Function
②複素数の絶対値
Public Function absCom(A As Complex) As Double
absCom = Sqr(A.R * A.R + A.I * A.I)
End Function
③加算
複素数演算のＶＢでの定義
Public Function addCom(A As Complex, B As Complex) As Complex
addCom.R = A.R + B.R
addCom.I = A.I + B.I
End Function
Public Function addComReal(A As Complex, B As Double) As Complex
addComReal.R = A.R + B
addComReal.I = A.I
End Function
④減算と符号反転
複素数演算のＶＢでの定義
Public Function subCom(A As Complex, B As Complex) As Complex
subCom.R = A.R - B.R
subCom.I = A.I - B.I
End Function
Public Function subComReal(A As Complex, B As Double) As Complex
subComReal.R = A.R - B
subComReal.I = A.I
End Function
Public Function minusCom(A As Complex) As Complex
minusCom.R = -A.R
minusCom.I = -A.I
End Function
⑤乗算とべき乗
複素数演算のＶＢでの定義
Public Function multCom(A As Complex, B As Complex) As Complex
multCom.R = A.R * B.R - A.I * B.I
multCom.I = A.R * B.I + A.I * B.R
End Function
Public Function multComReal(A As Complex, B As Double) As Complex
multComReal.R = A.R * B
multComReal.I = A.I * B
End Function
Public Function powerCom(A As Complex, N) As Complex
Dim S As Complex
S = A
For I = 2 To N
S = multCom(S, A)
Next
powerCom = S
End Function
複素数演算のＶＢでの定義
⑥複素数の逆数
Public Function invCom(B As Complex) As Complex
ImageP = B.I: RealP = B.R
If Abs(ImageP) < 0.000001 Then
If Abs(RealP) < 0.000001 Then
MsgBox "逆数で0割りが起こりました"
invCom.R = 0: invCom.I = 0
Else
invCom.R = 1 / RealP: invCom.I = 0
End If
Else
S = RealP * RealP + ImageP * ImageP
If Abs(S) < 0.000001 Then
MsgBox "逆数で0割りが起こりました"
invCom.R = 0: invCom.I = 0
Else
invCom.R = RealP / S: invCom.I = -ImageP / S
End If
End If
End Function
複素数演算のＶＢでの定義
⑦複素数の除算
Public Function divCom(A As Complex, B As Complex) As Complex
divCom = multCom(A, invCom(B))
End Function
⑧複素数の平方根
Public Function sqrCom(A As Complex) As Complex
SS = Sqr(A.R * A.R + A.I * A.I)
sqrCom.R = Sqr((SS + A.R) / 2)
sqrCom.I = Sqr((SS - A.R) / 2)
End Function
⑨複素数のExponential
Public Function expCom(A As Complex) As Complex
EP = Exp(A.R): TH = A.I
expCom.R = Cos(TH) * EP: expCom.I = Sin(TH) * EP
End Function
複素数演算のＶＢでの定義
⑩複素数を文字列に変換
Public Function strCom(A As Complex) As String
S = A.I: SS = A.R
If Abs(S) < 0.000005 And Abs(SS) < 0.000005 Then
strCom = "0.0000"
ElseIf Abs(S) < 0.000005 Then
strCom = Format(A.R, "#0.0000")
ElseIf Abs(SS) < 0.000005 Then
strCom = Format(S, "#0.0000") & "*j"
Else
K = " - "
If S > 0 Then K = " + "
S = Abs(S)
strCom = Format(A.R, "#0.0000") & K & Format(S, "#0.0000") & "*j"
End If
End Function
［簡単な例題］
• 次の式の解は複素数解となる。
z  z 1  0
2
例題を Excel でやってみよう
Ｅｘｃｅｌでの定義
グラフを描いて確かめよう
収束の様子
ＶＢＡでのプログラム
①処理本体
Function 関数値(X As Complex) As Complex
関数値 = addCom(addCom(powerCom(X, 2), X), setCom(1, 0))
End Function
Function 微分値(X As Complex) As Complex
微分値 = addCom(multCom(X, setCom(2, 0)), setCom(1, 0))
End Function
Function ComplexNewton(iter, E, EPS, iterMax) As Complex
Dim X As Complex: Dim F As Complex: Dim DF As Complex
iter = 0: X = setCom(0, 1)
Do While iter < iterMax
iter = iter + 1
F = 関数値(X): E = absCom(F)
If E < EPS Then Exit Do
DF = 微分値(X)
X = subCom(X, divCom(F, DF)) ' X(k) = X(k-1) - F(x)/F'(x)
Loop
ComplexNewton = X
End Function
ＶＢＡでのプログラム
②ボタンのクリックイベントハンドラ
Sub ボタン1_Click()
Dim X As Complex
X = ComplexNewton(iter, E, 0.00001, 500)
MsgBox "繰返し回数 = " & iter & " 結果 = " & strCom(X) & _
" 誤差=" & E
End Sub
課題１
• 次の方程式の複素数解をニュートンラプソ
ン法を用いて求めよ。
z  2z  3  0
2
（２）セカント法
複素数においても，セカント法を使うことができる。
zk 1  zk
zk 1  zk  f zk 
f zk 1   f zk 1 
例題を Excel でやってみよう
Ｅｘｃｅｌでの定義
グラフを描いて確かめよう
収束の様子
ＶＢＡでのプログラム
①処理本体
Private Function 関数値(X As Complex) As Complex
関数値 = addCom(addCom(powerCom(X, 2), X), setCom(1, 0))
End Function
Function ComplexSecant(iter, E, EPS, iterMax) As Complex
Dim Z1 As Complex: Dim Z2 As Complex
Dim F1 As Complex: Dim F2 As Complex
Dim DZ As Complex: Dim DF As Complex
iter = 0: Z1 = setCom(0, 1): Z2 = setCom(0.5, 0.5)
F1 = 関数値(Z1)
Do While iter < iterMax
iter = iter + 1
F2 = 関数値(Z2): E = absCom(F2)
If E < EPS Then Exit Do
DZ = subCom(Z1, Z2)
DF = subCom(F1, F2)
Z1 = Z2
F1 = F2
Z2 = subCom(Z2, multCom(F2, divCom(DZ, DF)))
Loop
ComplexSecant = Z2
End Function
ＶＢＡでのプログラム
②ボタンのクリックイベントハンドラ
Sub ボタン6_Click()
Dim X As Complex
X = ComplexSecant(iter, E, 0.00001, 500)
MsgBox "繰返し回数 = " & iter & " 結果 = " & strCom(X) & _
" 誤差=" & E
End Sub
課題２
• 次の方程式の複素数解をセカント法を用
いて求めよ。
z  z20
2
セカント法は，
z1    f z1  


z2    f z2 
という連立方程式を解き， z  
①
（３）ミュラー法
（セカント法の拡張）
 0としたときの z
を新しい近似値
とすることと同等である。
①式を解くと
f z1   f z2 
f z2 z1  f z1 z2

, 
z1  z2
z1  z2

f z2 z1  f z1 z2
f z2 z1  f z1 z2
z    
 z2  z2 

f z1   f z2 
f z1   f z2 
f z1 z2  f z2 z2  f z2 z1  f z1 z2
 z2 
f z1   f z2 
f z2 z1  f z2 z2
f z2 z1  z2 
 z2 
 z2 
f z1   f z2 
f z1   f z2 
z3
（３）ミュラー（Muller）法
（セカント法の拡張）
z12  z1    f z1  
 2

z2  z2    f z2  ①
z 2  z    f z 
3
3 
 3
・・
連立方程式
を解いて， ,  ,  を求め， z
2
 z    0 の根を新しい近似値 z4
とする。２根が求まるが， z3 に近いほうを採用する。
［注］
セカント法やニュートンラプソン法では，
近似値として複素数または虚数を近似値として与えないと，
複素根を求めることができないが，
ミュラー法では２次方程式を解く過程で複素数が求まるので，
初期値は，実数でも構わない。
ＶＢＡでのプログラム（Excelでは煩雑になるのでＶＢＡのみ示す）
① 宣言と方程式の値，複素数型の２次方程式の根，
Public Type 根データ
R1 As Complex
R2 As Complex
End Type
Function 方程式(X As Complex) As Complex
方程式 = addCom(powerCom(X, 2), addComReal(X, 1))
End Function
Function 二次方程式の根(A As Complex, B As Complex, C As Complex) As 根データ
Dim BAC As Complex: Dim A2 As Complex: Dim B2 As Complex
A2 = multComReal(A, 2): B2 = minusCom(B)
BAC = sqrCom(subCom(powerCom(B, 2), multCom(multComReal(A, 4), C)))
二次方程式の根.R1 = divCom(addCom(B2, BAC), A2)
二次方程式の根.R2 = divCom(subCom(B2, BAC), A2)
End Function
ＶＢＡでのプログラム ②処理本体（１）
Function Muller() As Complex
Dim X1 As Complex: Dim X2 As Complex: Dim X3 As Complex
Dim F1 As Complex: Dim F2 As Complex: Dim F3 As Complex
Dim DX1 As Complex: Dim DX2 As Complex: Dim DX3 As Complex
Dim DF1 As Complex: Dim DF2 As Complex
Dim A As Complex: Dim B As Complex: Dim C As Complex
Dim XX As Complex: Dim YY As Complex
Dim R As 根データ
EPS = 0.00001
iterMax = 500: iter = 0
' 初期値
X1 = setCom(0, 0): X2 = setCom(1, 0): X3 = setCom(2, 0)
F1 = 方程式(X1): F2 = 方程式(X2): F3 = 方程式(X3)
DX1 = subCom(X1, X2): DX2 = subCom(X2, X3): DX3 = subCom(X3, X1)
DF1 = subCom(F1, F2): DF2 = subCom(F2, F3)
Do While iter < iterMax
iter = iter + 1: E = absCom(F3)
If E < EPS Then Exit Do
ＶＢＡでのプログラム ③処理本体（２）
‘ 連立方程式を解く
‘ A=((F2-F3)*(X1-X2)-(F1-F2)*(X2-X3))/((X1-X2)*(X2-X3)*(X3-X1))
'
B=(F1-F2)/(X1-X2)-A*(X1+X2)
C=F1-B*X1-A*X1*X1
XX = subCom(multCom(DF2, DX1), multCom(DF1, DX2))
YY = multCom(multCom(DX1, DX2), DX3)
A = divCom(XX, YY)
B = subCom(divCom(DF1, DX1), multCom(A, addCom(X1, X2)))
C = subCom(subCom(F1, multCom(B, X1)), multCom(A, powerCom(X1, 2)))
‘ 二次方程式の根
R = 二次方程式の根(A, B, C)
X1 = X2: X2 = X3: F1 = F2: F2 = F3
DX1 = DX2: DX2 = DX3: DF1 = DF2
E1 = absCom(subCom(R.R1, X3)): E2 = absCom(subCom(R.R2, X3))
If E1 <= E2 Then X3 = R.R1 Else X3 = R.R2
F3 = 方程式(X3)
DX3 = subCom(X3, X1): DF2 = subCom(F2, F3)
Loop
Muller = X3
MsgBox strCom(X3)
End Function
ＶＢＡでのプログラム ④ ボタンのClickイベントハンドラ
Sub ボタン2_Click()
Dim R As Complex
R= Muller()
MsgBox strCom(R)
End Sub
課題３
• 次の方程式の複素数解をミュラー法を用
いて求めよ。
z  2z  2  0
2
・・
（４）複数の解（１）
1 を根として
f z 
f1 z  
という関数を考える。
z  1 
f z z  1   f z 
f1z  
2
z  1 
f z 
f1 z 

f z 
z  1 


f1z  f z z  1   f z 
 1 
f z   f z 
2

z  1 
z  1 
・・
（４）複数の解（２）
次に
1, 2 を根として
f z 
f 2 z  
z  1 z  2 
という関数を考える。
f z z  1 z   2   f z z  1   z   2 
f 2z  
z  1 2 z  2 2
f z 
f 2 z 

z  1 z   2 

f 2z  f z z  1 z   2   f z z  1   z   2 
z  1 2 z  2 2
f z 

 1
1 
f z   f z 






z


z


1
2 

・・
（４）複数の解（３）
n 個の根が既知であるとき
i (i  1, 2,, n) とすると，
一般に
根を
f z 
f n z  

z  1 z  2 z  2 
f n z 
f z 

n
f nz 
f z   f z 
f z 
n
 z   
i 1
i
1
j 1 z  1 
として，ニュートンラプソン法を適用すると，複数解を得ることができる。
ＶＢＡでのプログラム①宣言と方程式の値，微分値の変更
Public Z(3) As Complex
Function 関数値(X As Complex) As Complex
関数値 = addCom(addCom(powerCom(X, 3), _
addCom(multCom(powerCom(X, 2), setCom(2, 0)), X)), setCom(1, 0))
End Function
Function 微分値(X As Complex) As Complex
微分値 = addCom(multCom(powerCom(X, 2), setCom(3, 0)), _
addCom(multCom(X, setCom(4, 0)), setCom(1, 0)))
End Function
Function 変更微分値(DF As Complex, F As Complex, X As Complex, K) As Complex
Dim DFD As Complex
If K >= 2 Then
DFD = setCom(0, 0) ' F'-F*Σ(1/(X-Zi))
(i=1～N-1)の計算
For I = 1 To K - 1
DFD = addCom(DFD, invCom(subCom(X, Z(I))))
Next
変更微分値 = subCom(DF, multCom(F, DFD))
Else
変更微分値 = DF
End If
End Function
ＶＢＡでのプログラム ②処理本体
Sub multiComplexNewton(N)
Dim X As Complex: Dim F As Complex: Dim DF As Complex
EPS = 0.00001: iterMax = 500
For K = 1 To N
iter = 0: X = setCom(1, 1)
Do While iter < iterMax
iter = iter + 1: F = 関数値(X): E = absCom(F)
If E < EPS Then Exit Do
DF = 変更微分値(微分値(X), F, X, K)
X = subCom(X, divCom(F, DF))
Loop
Z(K) = X
Next
End Sub
ＶＢＡでのプログラム ③ボタンのクリックイベントハンドラ
Sub ボタン3_Click()
multiComplexNewton 3
S = ""
For I = 1 To 3
S = S & "結果(" & I & ") = " & strCom(Z(I)) & Chr(13) & Chr(10)
Next
MsgBox S
End Sub
・・
（５）ＤＫＡ法
ニュートン法を，デュラン（Durand），ケルナー（ Kerner ）が修正し，
同方法に対してアバース（Aberth）の初期値を用いる方法。
それぞれの頭文字からＤＫＡ法と呼ばれる。
z
( m1)
j
z
( m)
j

 
f z (jm)
z
n
k 1
k j
( m)
k
 z (jm)

, ( j  1, 2, , n)
・・
ニュートンラプソン法との違い
ニュートン・ラプソン法では，
１個ずつ解を求めるが，
ＤＫＡ法では，n個の解をまとめて求める。
初期値は，以下のようにとる。
z
i
( 0)
j
  2( j 1)  
a1

 r  expi
 
na0
n
2n 

：虚数単位,
r：複素平面上で
a1

を中心としてすべての根を囲む
na0
ような円の半径
ＶＢＡでのプログラム①宣言，データ設定，結果表示
Private Const N = 4
Private A(N) As Double
Private X(N - 1) As Complex
Sub データの設定()
A(0) = 1: A(1) = -9: A(2) = -7: A(3) = -35:
A(4) = 50
End Sub
Sub 結果表示()
S = ""
For I = 0 To N - 1
S = S & strCom(X(I)) & Chr(13) & Chr(10)
Next
MsgBox S
End Sub
Z 4  9z3  7z 2  35x  50  0
を示す。
ＶＢＡでのプログラム②処理用サブルーチン
Sub 係数初期化(A, N)
For j = 1 To N
A(I) = A(I) / A(0)
Next
End Sub
Function 初期値用半径(A, N) As Double
R0 = 0
For j = 2 To N
Rn = N * (Abs(A(j)) ^ (1# / j))
If Rn > R0 Then R0 = Rn
Next
初期値用半径 = R0
End Function
Sub 初期値計算(R0, N)
Pi = 3.1415926
For j = 0 To N - 1
Th = 2# * Pi * j / N + Pi / (2# * N)
X(j) = setCom(R0 * Cos(Th), R0 * Sin(Th))
Next
End Sub
ＶＢＡでのプログラム
③処理本体
Function DKA(A, N, delta, iterMax) As Boolean
Dim Xkj As Complex: Dim Fx As Complex
Dim beforX As Complex: Dim Adj As Complex
係数初期化 A, N
初期値計算初期値用半径(A, N), N
For I = 0 To iterMax
MaxError = 0
For j = 0 To N - 1
Xkj = setCom(1, 0): Fx = setCom(1, 0): beforX = X(j)
For K = 0 To N - 1
Fx = addCom(multCom(Fx, beforX), setCom(A(K + 1), 0))
If K <> j Then Xkj = multCom(Xkj, subCom(beforX, X(K)))
Next
Adj = divCom(Fx, Xkj): X(j) = subCom(beforX, Adj)
E = absCom(Fx)
If E > MaxError Then MaxError = E
Next
If MaxError < delta Then
DKA = False: Exit Function
End If
Next
DKA = True
End Function
ＶＢＡでのプログラム
④ ボタンのClickイベントハンドラ
Sub ボタン7_Click()
データの設定
R = DKA(A, N, 0.00000001, 500)
If R Then
MsgBox "収斂しません"
Else
結果表示
End If
End Sub
課題４
• 次の方程式のすべての解をＤＫＡ法を用い
て求めよ。
z  8z  6z  32z  40  0
4
3
2

Download Report