プログラミング言語によっては，複素数が使えない。 ４．３　複素数解 （１）ニュートン・ラプソン法 複素数においても，ニュートン・ラプソンの方法を使うことができる。 プログラミング言語によっては，複素数が使えない。 自分で定義

複素数演算のＶＢでの定義 ①複素数型のユーザ定義と複素数型定数の設定 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

複素数演算のＶＢでの定義 ④減算と符号反転 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 Public Function minusCom(A As Complex) As Complex minusCom.R = -A.R minusCom.I = -A.I

複素数演算のＶＢでの定義 ⑤乗算とべき乗 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 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

複素数演算のＶＢでの定義 ⑥複素数の逆数 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 S = RealP * RealP + ImageP * ImageP If Abs(S) < 0.000001 Then invCom.R = RealP / S: invCom.I = -ImageP / S 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

［簡単な例題］ 次の式の解は複素数解となる。

Ｅｘｃｅｌでの定義 例題を 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)) 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

ＶＢＡでのプログラム　　　②ボタンのクリックイベントハンドラ Sub ボタン1_Click() Dim X As Complex X = ComplexNewton(iter, E, 0.00001, 500) MsgBox "繰返し回数 = " & iter & " 結果 = " & strCom(X) & _ " 誤差=" & E End Sub

課題１ 次の方程式の複素数解をニュートンラプソン法を用いて求めよ。

（２）セカント法 複素数においても，セカント法を使うことができる。

Ｅｘｃｅｌでの定義 例題を Excel でやってみよう

収束の様子 グラフを描いて確かめよう

ＶＢＡでのプログラム ①処理本体 Private Function 関数値(X As Complex) As Complex ＶＢＡでのプログラム　　　①処理本体 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

ＶＢＡでのプログラム ②ボタンのクリックイベントハンドラ ＶＢＡでのプログラム　　　②ボタンのクリックイベントハンドラ Sub ボタン6_Click() Dim X As Complex X = ComplexSecant(iter, E, 0.00001, 500) MsgBox "繰返し回数 = " & iter & " 結果 = " & strCom(X) & _ " 誤差=" & E End Sub

課題２ 次の方程式の複素数解をセカント法を用いて求めよ。

（３）ミュラー法 （セカント法の拡張） セカント法は， ① という連立方程式を解き， としたときの を新しい近似値 とすることと同等である。 ①式を解くと

（３）ミュラー（Muller）法 （セカント法の拡張） ・・ 連立方程式 ① の根を新しい近似値 を解いて， を求め， とする。２根が求まるが， に近いほうを採用する。 ［注］ セカント法やニュートンラプソン法では， 近似値として複素数または虚数を近似値として与えないと， 複素根を求めることができないが， ミュラー法では２次方程式を解く過程で複素数が求まるので， 初期値は，実数でも構わない。

ＶＢＡでのプログラム（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)

ＶＢＡでのプログラム ②処理本体（１） Function Muller() As Complex ＶＢＡでのプログラム　②処理本体（１） 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

課題３ 次の方程式の複素数解をミュラー法を用いて求めよ。

（４）複数の解（１） ・・ を根として という関数を考える。

（４）複数の解（２） ・・ 次に を根として という関数を考える。

（４）複数の解（３） 一般に 個の根が既知であるとき とすると， 根を として，ニュートンラプソン法を適用すると，複数解を得ることができる。 ・・ 一般に 個の根が既知であるとき 根を とすると， として，ニュートンラプソン法を適用すると，複数解を得ることができる。

ＶＢＡでのプログラム①宣言と方程式の値， 微分値の変更 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))) 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

ＶＢＡでのプログラム ②処理本体 Sub multiComplexNewton(N) ＶＢＡでのプログラム　②処理本体 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 ）が修正し， ・・ ニュートン法を，デュラン（Durand），ケルナー（ Kerner ）が修正し， 同方法に対してアバース（Aberth）の初期値を用いる方法。 それぞれの頭文字からＤＫＡ法と呼ばれる。

ニュートンラプソン法との違い ニュートン・ラプソン法では， １個ずつ解を求めるが， ＤＫＡ法では，n個の解をまとめて求める。 ・・ ニュートン・ラプソン法では， １個ずつ解を求めるが， ＤＫＡ法では，n個の解をまとめて求める。 初期値は，以下のようにとる。 ： 虚数単位, ： 複素平面上で を中心としてすべての根を囲む ような円の半径

ＶＢＡでのプログラム①宣言，データ設定，結果表示 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 を示す。

ＶＢＡでのプログラム②処理用サブルーチン 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 初期値用半径 = 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))

ＶＢＡでのプログラム ③処理本体 Function DKA(A, N, delta, iterMax) As Boolean ＶＢＡでのプログラム　　③処理本体 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 If MaxError < delta Then DKA = False: Exit Function End If DKA = True End Function

ＶＢＡでのプログラム ④ ボタンのClickイベントハンドラ Sub ボタン7_Click() データの設定 R = DKA(A, N, 0.00000001, 500) If R Then MsgBox "収斂しません" Else 結果表示 End If End Sub

課題４ 次の方程式のすべての解をＤＫＡ法を用いて求めよ。