Private A As Double Private D As Double
②データの設定部分 Private A() As Double Private D() As Double Private X() As Double Function べき乗データ設定() N = 4 Re. Dim A(N, N), X(N) With Worksheets("Sheet 1") For i = 1 To N For j = 1 To N A(i, j) = Val(. Cells(i + 1, j + 1)) Next End With べき乗データ設定 = N End Function
④処理本体(1) Public Sub べき乗法(A, EPS, X, λ, N) Dim i, K 1, K 2, Iter. Max As Integer ' Dim T As Double Dim XTemp() As Double Iter. Max = 100 Re. Dim XTemp(2, N) K 1 = 2: K 2 = 1 XTemp(K 2, 1) = 1 For i = 2 To N XTemp(K 2, i) = 0 Next E = EPS * 100 Iter = 0 Do While Iter < Iter. Max And E > EPS Iter = Iter + 1 KK = K 2: K 2 = K 1: K 1 = KK ' K 1:前回計算,K 2: 今回 (続きあり)
④処理本体(2) ' 行列式の乗算 For i = 1 To N T = 0 For j = 1 To N T = T + A(i, j) * XTemp(K 1, j) Next XTemp(K 2, i) = T Next 'ノルム 1=1による方法 T = 0 For i = 1 To N T = T + XTemp(K 2, i) * XTemp(K 2, i) Next λ = Sqr(T) If Abs(λ) < EPS Then '固有値が0の場合は収束しないものとする E = EPS * 100 Exit Do End If For i = 1 To N XTemp(K 2, i) = XTemp(K 2, i) / λ Next (続きあり)
④処理本体(3) T 1 = 0: T 2 = 0 For i = 1 To N DX = XTemp(K 2, i) - XTemp(K 1, i) T 1 = T 1 + DX * DX T 2 = T 2 + XTemp(K 2, i) * XTemp(K 2, i) Next If Abs(T 2) < EPS Then '固有ベクトルが0のとき収束しないものとする E = EPS * 100 Exit Do End If E = T 1 / T 2 Loop For i = 1 To N X(i) = XTemp(K 2, i) Next If E > EPS Then Msg. Box "収束しません " & λ & " 繰返し回数" & Iter & " E=" & E Else Msg. Box λ & " 繰返し回数" & Iter & " E=" & E End If End Sub
②データの設定部分 Private A() As Double Private D() As Double Private X() As Double Function データ設定() N = 4 Re. Dim A(N, N + N), X(N) With Worksheets("Sheet 1") For i = 1 To N For j = 1 To N A(i, j) = Val(. Cells(i + 1, j + 1)) Next End With データ設定 = N End Function
③結果の設定 Sub 最小固有値結果設定(A, X, λ, N) With Worksheets("Sheet 1"). Cells(2, 16) = λ For i = 1 To N. Cells(i + 1, 17) = X(i) Next For i = 1 To N For j = 1 To N. Cells(i + 11, j + 1) = A(i, j) Next End With End Sub
②データの設定部分 Private A() As Double Private C() As Double Private X() As Double Private LU() As Double Private Ipivot() As Integer Function 逆反復法データ設定() N = 4 Re. Dim A(N, N), X(N) With Worksheets("Sheet 1") For i = 1 To N For j = 1 To N A(i, j) = Val(. Cells(i + 1, j + 1)) Next End With 逆反復法データ設定 = N End Function
④処理本体(1) Sub 逆反復法(A, X, λ, N) Dim i, j, k, Iter. Max As Integer Dim T As Double Dim LU() As Double Dim B() As Double Dim Ipivot() As Integer Re. Dim B(N) Re. Dim LU(N, N) Re. Dim Ipivot(N) EPS = 0. 00000001 LU分解 A, LU, Ipivot, N, EPS X(1) = 1 For i = 2 To N X(i) = 0 Next Iter. Max = 20 Iter = 0: λ = 1 E = EPS * 100 (続きあり)
④処理本体(2) Do While E > EPS And Iter < Iter. Max Iter = Iter + 1 For i = 1 To N ‘前進代入 T = X(Ipivot(i)) For j = 1 To i - 1 T = T - LU(Ipivot(i), j) * B(j) Next B(i) = T Next ‘後退代入 For i = N To 1 Step -1 T = B(i) For j = i + 1 To N T = T - LU(Ipivot(i), j) * B(j) Next B(i) = T / LU(Ipivot(i), i) Next (続きあり)
④処理本体(3) If Abs(T 2) < EPS Then '固有ベクトルが0のとき収束しない E = EPS * 100 Exit Do End If E = T 1 / T 2 Loop Msg. Box λ & " " & E & " " & Iter End Sub
②データの設定部分 Private A(4, 4) As Double Sub データ設定() With Worksheets("Sheet 1") For i = 1 To 4 For j = 1 To 4 A(i, j) =. Cells(i + 1, j + 1) Next End With End Sub
③結果の設定とイベントハンドラ Sub 結果設定(A) With Worksheets("Sheet 1") For i = 1 To 4 For j = 1 To 4. Cells(i + 1, j + 7) = A(i, j) Next End With End Sub ボタン 1_Click() データ設定 EPS = 0. 00001 R = LR分解固有値計算(A, EPS, 500, 4) 結果設定 A End Sub
④共通ルーチン Sub 行列乗算(A, B, C, N 1, N 2) For i = 1 To N 1 For j = 1 To N 2 C(i, j) = 0 For k = 1 To N 2 C(i, j) = C(i, j) + A(i, k) * B(k, j) Next End Sub
⑤処理本体 Function LR分解固有値計算(A, EPS, iter. Max, N) Dim L() As Double: Re. Dim L(N, N) Dim U() As Double: Re. Dim U(N, N) R = 0: E = EPS * 100: E 2 = EPS * EPS Do While R < iter. Max And E > E 2 R = R + 1 LU分解 A, L, U, N 行列乗算 U, L, A, N, N E = 0 For i = 2 To N For j = 1 To i - 1 E = E + A(i, j) * A(i, j) Next Loop LR分解固有値計算 = R End Function
③固有値計算 Function QR分解固有値計算(A, EPS, iter. Max, N) Dim Q() As Double: Re. Dim Q(N, N) Dim R() As Double: Re. Dim R(N, N) k = 0: E = EPS * 100: E 2 = EPS * EPS Do While k < iter. Max And E > E 2 k = k + 1 QR分解 A, Q, R, N 行列乗算 R, Q, A, N, N E = 0 For i = 2 To N For j = 1 To i - 1 E = E + A(i, j) * A(i, j) Next Loop QR分解固有値計算 = k End Function
②データの設定部分 Private A() As Double Private D() As Double Function データ設定() N = 4 Re. Dim A(N, N), D(N) With Worksheets("Sheet 1") For i = 1 To N For j = 1 To N A(i, j) = Val(. Cells(i + 1, j + 1)) Next End With データ設定 = N End Function
③結果の設定とイベントハンドラ Sub 結果表示(A, D, N) With Worksheets("Sheet 1") For i = 1 To N For j = 1 To N. Cells(j + 1, i + 7) = A(i, j) Next. Cells(i + 1, 13) = D(i) Next End With End Sub ボタン 1_Click() N = データ設定() Jacobi A, D, N 結果表示 A, D, N End Sub
④処理本体(1) Sub Jacobi(A, D, N) Dim i, j, k As Integer Dim T, C, S, X, Y As Double Dim Iter, Max. Iter As Integer Dim OK As Boolean Dim V() As Double: Re. Dim V(N) Dim W() As Double: Re. Dim W(N, N) For k = 1 To N For j = 1 To N W(k, j) = 0 Next W(k, k) = 1: D(k) = A(k, k) Next Iter = 0: Max. Iter = 50 Do Iter = Iter + 1: OK = True For j = 1 To N - 1 For k = j + 1 To N T = Abs(D(j)) + Abs(D(k)) If Abs(A(j, k)) + T <> T Then OK = False (続きあり)
④処理本体(2) T = (D(k) - D(j)) / (2 * A(j, k)) If T >= 0 Then T = 1 / (T + Sqr(T * T + 1)) _ Else T = 1 / (T - Sqr(T * T + 1)) C = 1 / Sqr(T * T + 1): S = T * C: T = T * A(j, k) D(j) = D(j) - T: D(k) = D(k) + T: A(j, k) = 0 For i = 1 To j - 1 X = A(i, j): Y = A(i, k) A(i, j) = X * C - Y * S: A(i, k) = X * S + Y * C Next For i = j + 1 To k - 1 X = A(j, i): Y = A(i, k) A(j, i) = X * C - Y * S: A(i, k) = X * S + Y * C Next For i = k + 1 To N X = A(j, i): Y = A(k, i) A(j, i) = X * C - Y * S: A(k, i) = X * S + Y * C Next For i = 1 To N X = W(j, i): Y = W(k, i) W(j, i) = X * C - Y * S: W(k, i) = X * S + Y * C Next End If
④処理本体(3) Next Loop Until OK Or (Iter > Max. Iter) If OK Then For i = 1 To N - 1 k = i For j = i + 1 To N If D(j) > D(k) Then k = j Next T = D(i): D(i) = D(k): D(k) = T For j = 1 To N T = W(i, j): W(i, j) = W(k, j): W(k, j) = T Next A = W Else Msg. Box "収斂しません" End If End Sub
②データの設定部分 Public A() As Double Public D() As Double Public E() As Double Function データ設定() N = 4 Re. Dim A(N, N), D(N), E(N) With Worksheets("Sheet 1") For I = 1 To N For J = 1 To N A(I, J) = Val(. Cells(I + 1, J + 1)) Next End With データ設定 = N End Function
③三重対角化結果の設定とイベントハンドラ Sub 結果表示(A, D, E, N) With Worksheets("Sheet 1") For I = 1 To N. Cells(I + 1, 14) = D(I). Cells(I + 1, 15) = E(I) Next End With End Sub ボタン 1_Click() N = データ設定() EPS = 0. 00001 三重対角化 A, D, E, EPS, N 結果表示 A, D, E, N End Sub
④三重対角処理本体(1) Sub 三重対角化(A, D, E, EPS, N) Dim I, J, K As Integer Dim S As Double Dim V() As Double Dim W() As Double Re. Dim V(N), W(N) For K = 1 To N - 2 D(K) = A(K, K) For I = K + 1 To N W(I) = A(I, K) Next S = 0 For I = K + 1 To N S = S + W(I) * W(I) Next (続きあり)
④三重対角処理本体(2) If Abs(S) > EPS Then If W(K + 1) >= 0 Then S = Sqr(S) Else S = -Sqr(S) E(K + 1) = -S: W(K + 1) = W(K + 1) + S S = 1 / Sqr(W(K + 1) * S) For I = K + 1 To N W(I) = W(I) * S Next For I = K + 1 To N S = 0 For J = K + 1 To I S = S + A(I, J) * W(J) Next For J = I + 1 To N S = S + A(J, I) * W(J) Next V(I) = S Next (続きあり)
④三重対角処理本体(4) D(N - 1) = A(N - 1, N - 1): D(N) = A(N, N) E(1) = 0: E(N) = A(N, N - 1) '以下,固有ベクトル計算用 For K = N To 1 Step -1 If K < N - 1 Then For I = K + 1 To N W(I) = A(I, K) Next For I = K + 1 To N S = 0 For J = K + 1 To N S = S + A(I, J) * W(J) Next V(I) = S Next For I = K + 1 To N For J = K + 1 To N A(I, J) = A(I, J) - V(I) * W(J) Next End If (続きあり)
④三重対角処理本体(5) For I = 1 To N A(I, K) = 0 Next A(K, K) = 1 Next End Sub
⑤固有値と固有値ベクトルのためのデータ設定 Function 固有値用データ設定() N = 4 Re. Dim A(N, N), D(N) With Worksheets("Sheet 1") For I = 1 To N For J = 1 To N A(I, J) = Val(. Cells(I + 1, J + 1)) Next End With 固有値用データ設定 = N End Function
⑦固有値と固有値ベクトルの処理本体(1) Function Zero(D, E, I) As Boolean If I > 1 Then T = Abs(D(I)) + Abs(D(I - 1)) Else T = 0 Zero = ((T + Abs(E(I))) = T) End Function Sub 三重対角化による固有値(A, D, N) Dim Max. Iter, I, K, H, L As Integer Dim S, C, T, W, X, Y, Z As Double Dim E() As Double Max. Iter = 20 Re. Dim E(N) EPS = 0. 00001 三重対角化 A, D, E, EPS, N E(1) = 1: H = N Do While Zero(D, E, H) H = H - 1 Loop (続きあり)
⑦固有値と固有値ベクトルの処理本体(2) Do While H > 1 E(1) = 0: L = H - 1 Do While Not Zero(D, E, L) L = L - 1 Loop Iter = 0 Do Iter = Iter + 1 If Iter > Max. Iter Then Msg. Box "収斂しません" Exit Sub End If W = (D(H - 1) - D(H)) / 2 S = Sqr(W * W + E(H) * E(H)) If W < 0 Then S = -S X = D(L) - D(H) + E(H) * E(H) / (W + S) Y = E(L + 1) Z = 0 (続きあり)
⑦固有値と固有値ベクトルの処理本体(3) For K = L To H - 1 If Abs(X) >= Abs(Y) Then T = -Y / X: C = 1 / Sqr(T * T + 1): S = T * C Else T = -X / Y: S = 1 / Sqr(T * T + 1) If T < 0 Then S = -S C = T * S End If W = D(K) - D(K + 1): T = (W * S + 2 * C * E(K + 1)) * S D(K) = D(K) - T: D(K + 1) = D(K + 1) + T E(K) = C * E(K) - S * Z E(K + 1) = E(K + 1) * (C * C - S * S) + W * S * C For I = 1 To N X = A(K, I): Y = A(K + 1, I) A(K, I) = C * X - S * Y: A(K + 1, I) = S * X + C * Y Next If K < N - 1 Then X = E(K + 1): Y = -S * E(K + 2): Z = Y E(K + 2) = C * E(K + 2) End If Next Loop Until Zero(D, E, H) (続きあり)
⑦固有値と固有値ベクトルの処理本体(4) E(1) = 1: H = H - 1 Do While Zero(D, E, H) H = H - 1 Loop '以下,固有ベクトル For K = 1 To N - 1 L = K For I = K + 1 To N If D(I) > D(L) Then L = I Next T = D(K): D(K) = D(L): D(L) = T For I = 1 To N T = A(K, I): A(K, I) = A(L, I): A(L, I) = T Next End Sub
- Slides: 88