共振を防ぐように設計を行ったり， 振動を早く減衰させる設計を行う際， 固有値と固有ベクトルを求めることが重要

Presentation on theme: "共振を防ぐように設計を行ったり， 振動を早く減衰させる設計を行う際， 固有値と固有ベクトルを求めることが重要"— Presentation transcript:

1 共振を防ぐように設計を行ったり， 振動を早く減衰させる設計を行う際， 固有値と固有ベクトルを求めることが重要
３．固有値と固有ベクトル 共振を防ぐように設計を行ったり， 振動を早く減衰させる設計を行う際， 固有値と固有ベクトルを求めることが重要

2 （１）固有値と固有ベクトルの定義 以下の を固有値， を固有ベクトルという。 N行N列の正方行列 について （１） となる解のうち
なる解を考える。（１）式を変形すると， （２） ： N行N列の単位行列 係数行列 が逆行列を持つと，解は になる。 となる解を持つためには が逆行列を持たない， すなわち，行列式が零でならなければならない。したがって， （３）

3 式（１）の解釈 （１） ベクトルに行列 の変換操作を行っても，向きは変わらず，大きさだけが 倍となるベクトル が存在する。
固有ベクトルと呼ぶが，次のように無限のベクトルが考えられる。 そこで，単位ベクトルを の代表として選ぶ。

4 固有ベクトルを使った対角化 ベクトルに固有値ベクトルを並べた を使うと として対角化することができる。

5 簡単な行列を使った具体例（１） とすると， となる。 すなわち， 従って，固有値は または

6 簡単な行列を使った具体例（２） のとき 一方， から （符号が逆転していることに注意）

7 簡単な行列を使った具体例（２） のとき 同様に， から （同一符号であることに注意）

8 ２次元平面で図示するとλ＝１，λ＝３のいずれの組をとっても
直交することが分かる。 固有ベクトル同士の乗算を行うと直交性の確認ができる。

9 簡単な行列を使った具体例（３） 固有ベクトルを並べた行列 を使って を計算すると

10 外力がない場合の振動解析，すなわち自由振動問題となる 振動系そのものが固有に有する特性を求めること
（２）固有値問題の意味づけ 外力がない場合の運動方程式 を考える。 ここで ：質量行列， ：剛性行列， ：変位ベクトル， ：加速度ベクトル とすると この問題を解くのが固有値問題そのものである。 言い換えると 外力がない場合の振動解析，すなわち自由振動問題となる 振動系そのものが固有に有する特性を求めること

11 （３）固有値および固有ベクトルの求め方 ①最低次の振動モード 周期が長く，振幅が大きく，減衰しにくい。 ②高周波の振動モード
　　　周期が長く，振幅が大きく，減衰しにくい。 ②高周波の振動モード 　　　　騒音等の検討で重要。 したがって，工学的には， 最小固有値と最大固有値が重要

12 手法の分類 大きくは，以下の２つに分類できる。 ①反復法 　　（繰返し計算により求める） ②行列変換法 　　（行列の変換により求める）

13 反復法 ［一部の固有値のみ。固有ベクトルも求まる］ ①べき乗法（最大固有値を求める方法。累乗法ともいう）
　　　逆行列を計算し，その最大固有値を求め逆数をとることで， 　　最小固有値を計算できる。 ②逆反復法（最小固有値を求める方法） ③行列サブスペース法（逆反復法の拡張。大次元問題に対して， 　　　ごく一部の固有値を求める。なお，手法の複雑さに対して 　　　ごく一部の固有値しか求められないので，省略する） ［固有ベクトルは求まらないが，すべての固有値が求まる］ ④ＬＲ法（下三角行列と上三角行列に変換し，三角行列に収斂させる方法） ⑤ＱＲ法（直交行列と上三角行列に分けて，三角行列に収斂させる方法）

14 行列変換法 ①ヤコビ法 　　　　回転行列による直交変換を繰り返す。 ②ハウスホルダー法 　　　　三重対角行列に変換し，二分法を適用して 　　　　固有値を計算する。

15 ①べき乗法（Power Iteration Method)の考え方（１）
n行n列の行列 の固有値と固有ベクトルを とし， の順に並んでいるものとする。 固有ベクトルの任意の組は，正規直交基底となるから任意のベクトル は次のように表すことができる。 両辺に を乗じて この処理を 回繰り返すと

16 ①べき乗法（Power Iteration Method)の考え方（２）
を大きくしていくと， であるから， すなわち， となるので， と の比は に近づく。

17 ①べき乗法（Power Iteration Method)の考え方（３）
［処理手順］ （１）初期値として，最大成分が１であるような列ベクトルを設定する。例えば， （２） を計算する。 （３） とし， として単位ベクトル化する。 （４） により収束を判定し，収束していなければ （２）～（４）を繰り返す。

18 ＶＢでプログラムを作る ①シートの定義（λx, Axの計算の式定義を入れておくこと）

19 ②データの設定部分 Private A() As Double Private D() As Double Private X() As Double Function べき乗データ設定() N = 4 ReDim A(N, N), X(N) With Worksheets("Sheet1") 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

20 ③結果の設定とイベントハンドラ Sub べき乗法結果設定(X, λ, N) With Worksheets("Sheet1") .Cells(2, 14) = λ For i = 1 To N .Cells(i + 1, 15) = X(i) Next End With End Sub Sub ボタン5_Click() N = べき乗データ設定() EPS = べき乗法 A, EPS, X, λ, N べき乗法結果設定 X, λ, N

21 ④処理本体（１） （続きあり） Public Sub べき乗法(A, EPS, X, λ, N)
Dim i, K1, K2, Iter, IterMax As Integer ' Dim T As Double Dim XTemp() As Double IterMax = 100 ReDim XTemp(2, N) K1 = 2: K2 = 1 XTemp(K2, 1) = 1 For i = 2 To N XTemp(K2, i) = 0 Next E = EPS * 100 Iter = 0 Do While Iter < IterMax And E > EPS Iter = Iter + 1 KK = K2: K2 = K1: K1 = KK ' K1：前回計算，K2:今回 （続きあり）

22 ④処理本体（２） （続きあり） ' 行列式の乗算 For i = 1 To N T = 0 For j = 1 To N
T = T + A(i, j) * XTemp(K1, j) Next XTemp(K2, i) = T 'ノルム1=1による方法 T = T + XTemp(K2, i) * XTemp(K2, i) λ = Sqr(T) If Abs(λ) < EPS Then '固有値が0の場合は収束しないものとする E = EPS * 100 Exit Do End If XTemp(K2, i) = XTemp(K2, i) / λ （続きあり）

23 ④処理本体（３） T1 = 0: T2 = 0 For i = 1 To N
DX = XTemp(K2, i) - XTemp(K1, i) T1 = T1 + DX * DX T2 = T2 + XTemp(K2, i) * XTemp(K2, i) Next If Abs(T2) < EPS Then '固有ベクトルが0のとき収束しないものとする E = EPS * 100 Exit Do End If E = T1 / T2 Loop X(i) = XTemp(K2, i) If E > EPS Then MsgBox "収束しません " & λ & " 繰返し回数" & Iter & " E=" & E Else MsgBox λ & " 繰返し回数" & Iter & " E=" & E End Sub

24 ②逆反復法（Inverse Iteration Method)の考え方
正則対称行列 の場合，逆行列を求め， とし， べき乗法により最大固有値を求め，その逆数をとれば最小固有値となる。 したがって，適当な を初期値として， を計算し， の絶対値の逆数をとればよい。 逆行列の計算のかわりに，以下の連立方程式を解いてもよい。 連立方程式の解法として三角分解を用いれば，計算量を減らすことができる。

25 （逆行列を求めた後，べき乗法を用いる方法）
ＶＢでプログラムを作る （逆行列を求めた後，べき乗法を用いる方法） ①シートの定義（べき乗法の定義に追加する）

26 ②データの設定部分 Private A() As Double Private D() As Double Private X() As Double Function データ設定() N = 4 ReDim A(N, N + N), X(N) With Worksheets("Sheet1") 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

27 ③結果の設定 Sub 最小固有値結果設定(A, X, λ, N) With Worksheets("Sheet1") .Cells(2, 16) = λ For i = 1 To N .Cells(i + 1, 17) = X(i) Next For j = 1 To N .Cells(i + 11, j + 1) = A(i, j) End With End Sub

28 ④イベントハンドラ（べき乗法を呼び出す）
Sub ボタン6_Click() N = データ設定 EPS = ‘ 掃出法については連立方程式の項参照 If 掃出法(A, N, EPS) Then Exit Sub For i = 1 To N ‘ N+1列目移行に逆行列が入っているので移動 For j = 1 To N A(i, j) = A(i, j + N) Next EPS = 0.001 べき乗法 A, EPS, X, λ, N ‘ べき乗法の実行 最小固有値結果設定 A, X, 1 / λ, N End Sub

29 収束しない？ …だけど結果はよさそう このような場合， 収束条件について検討する必要がある。 すなわち，逆行列をとって計算する場合は，
収束しない？　…だけど結果はよさそう このような場合， 収束条件について検討する必要がある。 すなわち，逆行列をとって計算する場合は， 収束条件を別の方法にする必要がある。 たとえば，固有値が収束しているかどうか等に切り替える。

30 を解く方法 同じ行列で，何回も連立方程式を解く場合はＬＵ分解を用いる。 すなわち，
として，Lを用いた前進代入，Uを用いた後退代入を行って連立方程式を解くと，計算量が減る。

31 ＶＢでプログラムを作る による方法 ①シートの定義（べき乗法の定義に追加）

32 ②データの設定部分 Private A() As Double Private C() As Double Private X() As Double Private LU() As Double Private Ipivot() As Integer Function 逆反復法データ設定() N = 4 ReDim A(N, N), X(N) With Worksheets("Sheet1") 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

33 ③結果の設定とイベントハンドラ Sub 逆反復法結果表示(X, λ, N) With Worksheets("Sheet1") .Cells(2, 18) = λ For j = 1 To N .Cells(j + 1, 19) = X(j) Next End With End Sub Sub ボタン9_Click() N = 逆反復法データ設定() 逆反復法 A, X, λ, N 逆反復法結果表示 X, 1 / λ, N

34 ④処理本体（１） Sub 逆反復法(A, X, λ, N) Dim i, j, k, Iter, IterMax As Integer Dim T As Double Dim LU() As Double Dim B() As Double Dim Ipivot() As Integer ReDim B(N) ReDim LU(N, N) ReDim Ipivot(N) EPS = LU分解 A, LU, Ipivot, N, EPS X(1) = 1 For i = 2 To N X(i) = 0 Next IterMax = 20 Iter = 0: λ = 1 E = EPS * 100 （続きあり）

35 ④処理本体（２） Do While E > EPS And Iter < IterMax 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 B(i) = T / LU(Ipivot(i), i) （続きあり）

36 ④処理本体（３） 　　T = ‘ Bの絶対値を固有値λとする For i = 1 To N T = T + B(i) * B(i) Next λ = Sqr(T) If λ < EPS Then ‘ 固有値0のとき収束しない E = EPS * 100 Exit Do End If For i = 1 To N　　　　　　‘ 固有ベクトルの計算 X(i) = B(i) / λ T1 = 0: T2 = ‘ 収束条件の計算 DX = X(i) - B(i) T1 = T1 + DX * DX T2 = T2 + X(i) * X(i) （続きあり）

37 ④処理本体（３） If Abs(T2) < EPS Then '固有ベクトルが0のとき収束しない E = EPS * 100 Exit Do End If E = T1 / T2 Loop MsgBox λ & " " & E & " " & Iter End Sub

38 ■発案者のルーティスハウザー（Lutishauser）が， 下三角行列をＬ行列（Left triangular matrix）,
③ＬＲ分解による固有値演算 ■発案者のルーティスハウザー（Lutishauser）が， 下三角行列をＬ行列（Left triangular matrix）, 上三角行列をＲ行列（ Right triangular matrix ） と呼んでいたため，こう呼ばれる。 ■固有値問題では， ルーティスハウザーにならって， 　　ＬＵ分解のことをＬＲ分解と呼ぶ。

39 ＬＲ分解の考え方 固有値の定義から 左辺の係数行列 をＬＵ分解したときの下三角行列，上三角行列を とすると 左から
を乗じて，次のように置くこととする。 上式から，次の関係をもつ反復式を考えることができる。 ここで は上三角行列に収束し，対角要素が固有値になることが， 知られている。これを利用して固有値を求める方法である。 　固有ベクトルを求めることはできないが， 　工学的には，固有値のみを求めれば良い場面も多く， 　そのような場合には，単純で分かりやすい方法である。

40 ＬＲ分解の手順 ①初期値を とする。 ②行列 に対して，ＬＵ分解し， を求め次の計算を行う。 ③行列 が上三角行列になるまで，②を繰り返す。
終了条件は，たとえば行列 の 行 列の要素を とすると とする。

41 ＶＢでプログラムを作る ①シートの定義 （単純化のためにＬＵ分解はクラウト法による方法を用いる）

42 With Worksheets("Sheet1") For i = 1 To 4 For j = 1 To 4
②データの設定部分 Private A(4, 4) As Double Sub データ設定() With Worksheets("Sheet1") For i = 1 To 4 For j = 1 To 4 A(i, j) = .Cells(i + 1, j + 1) Next End With End Sub

43 With Worksheets("Sheet1") For i = 1 To 4 For j = 1 To 4
③結果の設定とイベントハンドラ Sub 結果設定(A) With Worksheets("Sheet1") For i = 1 To 4 For j = 1 To 4 .Cells(i + 1, j + 7) = A(i, j) Next End With End Sub Sub ボタン1_Click() データ設定 EPS = R = LR分解固有値計算(A, EPS, 500, 4) 結果設定 A

44 C(i, j) = C(i, j) + A(i, k) * B(k, j) Next
④共通ルーチン Sub 行列乗算(A, B, C, N1, N2) For i = 1 To N1 For j = 1 To N2 C(i, j) = 0 For k = 1 To N2 C(i, j) = C(i, j) + A(i, k) * B(k, j) Next End Sub

45 Function LR分解固有値計算(A, EPS, iterMax, N)
⑤処理本体 Function LR分解固有値計算(A, EPS, iterMax, N) Dim L() As Double: ReDim L(N, N) Dim U() As Double: ReDim U(N, N) R = 0: E = EPS * 100: E2 = EPS * EPS Do While R < iterMax And E > E2 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

46 ■ＬＲ法と処理手順は同じであるが，Ｌ行列ではなく，直交行列を用いる。
④ＱＲ分解による固有値演算 ■ＬＲ法と処理手順は同じであるが，Ｌ行列ではなく，直交行列を用いる。 ■直交化法としては，先に示したグラム・シュミットの方法を使うこととする。 ■最近では，ＬＲ法よりＱＲ法がよく用いられている。

47 ＱＲ分解の考え方 固有値の定義から 左辺の係数行列 をＱＲ分解したときの直交行列と上三角行列を とすると 左から
を乗じて，次のように置くと， 上式から，次の関係をもつ反復式を考えることができる。 ＬＲ法と同様， は上三角行列に収束し，対角要素が固有値になることが， 知られている。これを利用して固有値を求める方法である。 手順は，ＬＲ法とほぼ同様である。

48 ＶＢでプログラムを作る ①シートの定義 （ＱＲ分解はグラム・シュミットの直交化法を用いる）
（データ設定，結果設定，行列乗算は，ＬＲ法と同じ）

49 Sub ボタン1_Click() データ設定 EPS = 0.00001 R = QR分解固有値計算(A, EPS, 500, 4)
②イベントハンドラ Sub ボタン1_Click() データ設定 EPS = R = QR分解固有値計算(A, EPS, 500, 4) MsgBox R 結果設定 A End Sub

50 Function QR分解固有値計算(A, EPS, iterMax, N)
③固有値計算 Function QR分解固有値計算(A, EPS, iterMax, N) Dim Q() As Double: ReDim Q(N, N) Dim R() As Double: ReDim R(N, N) k = 0: E = EPS * 100: E2 = EPS * EPS Do While k < iterMax And E > E2 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

51 （４）行列変換法 （すべての固有値，固有ベクトルが求まる） 　　①ヤコビ法 　　②ハウスホルダ法

52 ①ヤコビ法 ■実対称行列の固有値を求める方法 ■工学系で現れる行列のほとんどは，実対称行列である。 ［覚えておきたい定理］ 正方行列
, 正則行列 に対して， と （および ） の固有値は一致する。 （固有値を変換させない変換を相似変換という）

53 直交行列とは 正則行列 の転置行列 と逆行列 逆行列 が一致する行列。 すなわち， ［ヤコビ法］ に対して直交行列で相似変換を行い，
対角成分以外をすべて 0 にする方法。 対角成分の値が固有値となる。

54 考え方（1） 二つの行に回転の成分をもち，それ以外の対角成分が １ ，非対角成分が 0 となる直交行列を考える。 Ｚ ［最も簡単な例］ Ｘ
Ｙ Ｚ軸回りの回転を意味する。したがってベクトルの大きさは変わらない。

55 考え方（２） ［Aが対象行列の場合］ 2行1列目または1行2列目を0とおくと 一方， から のとき

56 ただし， 考え方（３） のとき とおけばよい。 その他の項は以下のようにして計算すればよい

57 考え方（４） ここで とおくと すなわち対角項は

58 ウイルキンソン（Wilkinson）の公式
二つの行に回転の成分をもち，それ以外の対角成分が １ ，非対角成分が 0 となる直交行列を考える。 第 p 行 第 q 行

59 ウイルキンソン（Wilkinson）の公式
と置いて， ① ② ③ のとき，①を適用すると

60 ＶＢでプログラムを作る ①シートの定義（Ａ・Ｗ=Ｄ・Ｗの式定義も行う）

61 ②データの設定部分 Private A() As Double Private D() As Double Function データ設定() N = 4 ReDim A(N, N), D(N) With Worksheets("Sheet1") 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

62 ③結果の設定とイベントハンドラ Sub 結果表示(A, D, N) With Worksheets("Sheet1") 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) End With End Sub Sub ボタン1_Click() N = データ設定() Jacobi A, D, N 結果表示 A, D, N

63 ④処理本体（１） （続きあり） Sub Jacobi(A, D, N) Dim i, j, k As Integer
Dim T, C, S, X, Y As Double Dim Iter, MaxIter As Integer Dim OK As Boolean Dim V() As Double: ReDim V(N) Dim W() As Double: ReDim 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) Iter = 0: MaxIter = 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 （続きあり）

64 ④処理本体（２） 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 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 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 End If

65 ④処理本体（３） Next Next Loop Until OK Or (Iter > MaxIter) 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 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 A = W Else MsgBox "収斂しません" End If End Sub

66 ②三重対角化法 ■固有値問題とは，いわば座標変換により行列を対角化すること。
■いきなり対角化することは困難であるが，三重対角化は有限ステップで可能。 ■三重対角化行列とは，対角項とその隣以外が零となるような行列である。 ■最も有名な方法にハウスホルダー（Householder）法がある。

67 鏡映変換とは ■面を鏡の面として，できた像の位置にベクトルを移す変換を鏡映変換という。 （ベクトルの大きさは変わらない）
　　（ベクトルの大きさは変わらない） ■ハウスホルダー変換は鏡映変換のひとつである。 ［簡単な例］

68 ハウスホルダ変換とは 変換法 ：単位行列 に垂直な面

69 三重対角化するハウスホルダ変換（１） 実対称行列 を列ベクトル を列とする行列とする。 ①第1段のハウスホルダ変換は ② であるから
とすると の符号は，桁落ちがおきないよう と同じにする。

70 三重対角化するハウスホルダ変換（２） ③ として計算できる。 ④ は， の第1項が0であるから， （＊は 0 と限らない部分） の形になる。
⑤　次の計算を行う。

71 三重対角化するハウスホルダ変換（３） ⑥ の形と の対称性から の第1行目も，これを転置した形になる。 ⑦ ところで， とおくと さらに
したがって をあらかじめ計算しておくと，単純化できる。

72 三重対角化するハウスホルダ変換（４） ⑧その結果， ⑨第2列目以降は，次のようになるので，1行・1列目に影響を与えない。
の第1行と第1列を除いた部分 ⑨第2列目以降は，次のようになるので，1行・1列目に影響を与えない。 したがって，1列目，第2列目を除いた小行列に対して処理することになる。

73 三重対角化するハウスホルダ変換（５） ⑩その結果，以下の行列が得られる。
この変換によっても固有値は変化しない。したがって固有値を求める前に，この変換を行っておくことで，Ｎが大きい場合，ヤコビ法より速く固有値を求めることができる。

74 ＶＢでプログラムを作る ①シートの定義（Ａ・Ｗ=Ｄ・Ｗの式定義も行う）

75 ②データの設定部分 Public A() As Double Public D() As Double Public E() As Double Function データ設定() N = 4 ReDim A(N, N), D(N), E(N) With Worksheets("Sheet1") 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

76 ③三重対角化結果の設定とイベントハンドラ
Sub 結果表示(A, D, E, N) With Worksheets("Sheet1") For I = 1 To N .Cells(I + 1, 14) = D(I) .Cells(I + 1, 15) = E(I) Next End With End Sub Sub ボタン1_Click() N = データ設定() EPS = 三重対角化 A, D, E, EPS, N 結果表示 A, D, E, N

77 ④三重対角処理本体（１） （続きあり） Sub 三重対角化(A, D, E, EPS, N) Dim I, J, K As Integer
Dim S As Double Dim V() As Double Dim W() As Double ReDim 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 S = S + W(I) * W(I) （続きあり）

78 ④三重対角処理本体（２） （続きあり） 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 S = 0 For J = K + 1 To I S = S + A(I, J) * W(J) For J = I + 1 To N S = S + A(J, I) * W(J) V(I) = S （続きあり）

79 ④三重対角処理本体（３） （続きあり） S = 0 For I = K + 1 To N S = S + W(I) * V(I) Next
V(I) = V(I) - S * W(I) For J = K + 1 To I A(I, J) = A(I, J) - W(I) * V(J) - W(J) * V(I) ' 固有ベクトルを求めるときWの内容をAに移す A(I, K) = W(I) Else E(K + 1) = 0 End If （続きあり）

80 ④三重対角処理本体（４） （続きあり） 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 S = 0 For J = K + 1 To N S = S + A(I, J) * W(J) V(I) = S A(I, J) = A(I, J) - V(I) * W(J) End If （続きあり）

81 ④三重対角処理本体（５） For I = 1 To N A(I, K) = 0 Next A(K, K) = 1 End Sub

82 ⑤固有値と固有値ベクトルのためのデータ設定
Function 固有値用データ設定() N = 4 ReDim A(N, N), D(N) With Worksheets("Sheet1") 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

83 ⑥固有値と固有値ベクトルのためのイベントハンドラ
Sub 固有値結果表示(A, D, N) With Worksheets("Sheet1") For I = 1 To N .Cells(I + 1, 13) = D(I) For J = 1 To N .Cells(I + 1, J + 7) = A(I, J) Next End With End Sub Sub ボタン5_Click() N = 固有値用データ設定() EPS = 三重対角化による固有値 A, D, N 固有値結果表示 A, D, N

84 ⑦固有値と固有値ベクトルの処理本体（１） （続きあり） 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 MaxIter, I, K, H, L As Integer Dim S, C, T, W, X, Y, Z As Double Dim E() As Double MaxIter = 20 ReDim E(N) EPS = 三重対角化 A, D, E, EPS, N E(1) = 1: H = N Do While Zero(D, E, H) H = H - 1 Loop （続きあり）

85 ⑦固有値と固有値ベクトルの処理本体（２） （続きあり） 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 > MaxIter Then MsgBox "収斂しません" 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 （続きあり）

86 ⑦固有値と固有値ベクトルの処理本体（３） （続きあり） 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) Loop Until Zero(D, E, H) （続きあり）

87 ⑦固有値と固有値ベクトルの処理本体（４） E(1) = 1: H = H - 1 Do While Zero(D, E, H)
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 End Sub

88 （Ｘ）演習 （問題１）次の行列の固有値と固有ベクトルをサンプルプログラムを利用して求めよ。