Download presentation
Presentation is loading. Please wait.
1
プログラミング言語によっては，複素数が使えない。
４．３ 複素数解 （１）ニュートン・ラプソン法 複素数においても，ニュートン・ラプソンの方法を使うことができる。 プログラミング言語によっては，複素数が使えない。 自分で定義
2
複素数演算のＶＢでの定義 ①複素数型のユーザ定義と複素数型定数の設定 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
3
複素数演算のＶＢでの定義 ③加算 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
4
複素数演算のＶＢでの定義 ④減算と符号反転
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
5
複素数演算のＶＢでの定義 ⑤乗算とべき乗 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
6
複素数演算のＶＢでの定義 ⑥複素数の逆数 Public Function invCom(B As Complex) As Complex
ImageP = B.I: RealP = B.R If Abs(ImageP) < Then If Abs(RealP) < 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) < Then invCom.R = RealP / S: invCom.I = -ImageP / S End Function
7
複素数演算のＶＢでの定義 ⑦複素数の除算 ⑧複素数の平方根
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
8
複素数演算のＶＢでの定義 ⑩複素数を文字列に変換
Public Function strCom(A As Complex) As String S = A.I: SS = A.R If Abs(S) < And Abs(SS) < Then strCom = "0.0000" ElseIf Abs(S) < Then strCom = Format(A.R, "#0.0000") ElseIf Abs(SS) < 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
9
［簡単な例題］ 次の式の解は複素数解となる。
10
Ｅｘｃｅｌでの定義 例題を Excel でやってみよう
11
収束の様子 グラフを描いて確かめよう
12
ＶＢＡでのプログラム ①処理本体 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
13
ＶＢＡでのプログラム ②ボタンのクリックイベントハンドラ
Sub ボタン1_Click() Dim X As Complex X = ComplexNewton(iter, E, , 500) MsgBox "繰返し回数 = " & iter & " 結果 = " & strCom(X) & _ " 誤差=" & E End Sub
14
課題１ 次の方程式の複素数解をニュートンラプソン法を用いて求めよ。
15
（２）セカント法 複素数においても，セカント法を使うことができる。
16
Ｅｘｃｅｌでの定義 例題を Excel でやってみよう
17
収束の様子 グラフを描いて確かめよう
18
ＶＢＡでのプログラム ①処理本体 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
19
ＶＢＡでのプログラム ②ボタンのクリックイベントハンドラ
ＶＢＡでのプログラム ②ボタンのクリックイベントハンドラ Sub ボタン6_Click() Dim X As Complex X = ComplexSecant(iter, E, , 500) MsgBox "繰返し回数 = " & iter & " 結果 = " & strCom(X) & _ " 誤差=" & E End Sub
20
課題２ 次の方程式の複素数解をセカント法を用いて求めよ。
21
（３）ミュラー法 （セカント法の拡張） セカント法は， ① という連立方程式を解き， としたときの を新しい近似値 とすることと同等である。
①式を解くと
22
（３）ミュラー（Muller）法 （セカント法の拡張）
・・ 連立方程式 ① の根を新しい近似値 を解いて， を求め， とする。２根が求まるが， に近いほうを採用する。 ［注］ セカント法やニュートンラプソン法では， 近似値として複素数または虚数を近似値として与えないと， 複素根を求めることができないが， ミュラー法では２次方程式を解く過程で複素数が求まるので， 初期値は，実数でも構わない。
23
ＶＢＡでのプログラム（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)
24
ＶＢＡでのプログラム ②処理本体（１） 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 = 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
25
ＶＢＡでのプログラム ③処理本体（２） ‘ 連立方程式を解く
ＶＢＡでのプログラム ③処理本体（２） ‘ 連立方程式を解く ‘ 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
26
ＶＢＡでのプログラム ④ ボタンのClickイベントハンドラ
Sub ボタン2_Click() Dim R As Complex R= Muller() MsgBox strCom(R) End Sub
27
課題３ 次の方程式の複素数解をミュラー法を用いて求めよ。
28
（４）複数の解（１） ・・ を根として という関数を考える。
29
（４）複数の解（２） ・・ 次に を根として という関数を考える。
30
（４）複数の解（３） 一般に 個の根が既知であるとき とすると， 根を として，ニュートンラプソン法を適用すると，複数解を得ることができる。
・・ 一般に 個の根が既知であるとき 根を とすると， として，ニュートンラプソン法を適用すると，複数解を得ることができる。
31
ＶＢＡでのプログラム①宣言と方程式の値， 微分値の変更
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
32
ＶＢＡでのプログラム ②処理本体 Sub multiComplexNewton(N)
ＶＢＡでのプログラム ②処理本体 Sub multiComplexNewton(N) Dim X As Complex: Dim F As Complex: Dim DF As Complex EPS = : 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
33
ＶＢＡでのプログラム ③ボタンのクリックイベントハンドラ
ＶＢＡでのプログラム ③ボタンのクリックイベントハンドラ 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
34
（５）ＤＫＡ法 ニュートン法を，デュラン（Durand），ケルナー（ Kerner ）が修正し，
・・ ニュートン法を，デュラン（Durand），ケルナー（ Kerner ）が修正し， 同方法に対してアバース（Aberth）の初期値を用いる方法。 それぞれの頭文字からＤＫＡ法と呼ばれる。
35
ニュートンラプソン法との違い ニュートン・ラプソン法では， １個ずつ解を求めるが， ＤＫＡ法では，n個の解をまとめて求める。
・・ ニュートン・ラプソン法では， １個ずつ解を求めるが， ＤＫＡ法では，n個の解をまとめて求める。 初期値は，以下のようにとる。 ： 虚数単位, ： 複素平面上で を中心としてすべての根を囲む ような円の半径
36
ＶＢＡでのプログラム①宣言，データ設定，結果表示
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 を示す。
37
ＶＢＡでのプログラム②処理用サブルーチン
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 = For j = 0 To N - 1 Th = 2# * Pi * j / N + Pi / (2# * N) X(j) = setCom(R0 * Cos(Th), R0 * Sin(Th))
38
ＶＢＡでのプログラム ③処理本体 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
39
ＶＢＡでのプログラム ④ ボタンのClickイベントハンドラ
Sub ボタン7_Click() データの設定 R = DKA(A, N, , 500) If R Then MsgBox "収斂しません" Else 結果表示 End If End Sub
40
課題４ 次の方程式のすべての解をＤＫＡ法を用いて求めよ。
Similar presentations
