[PR] コピー機

エクセルで操る！ヤコビ法による固有値計算、固有ベクトル計算

Excelを用いた科学技術計算が第2版になりました　30年10月！

amazon：Excelで操る! ここまでできる科学技術計算第2版

楽天： Excelで操る! ここまでできる科学技術計算 [ 神足　史人 ]

固有値と固有ベクトル

　n×n行列Aに対して

　　Ax = λx　　x ≠0 　　　・・・ (1)

を満たす時、λをAの固有値、xを固有ベクトルという。

たった1日で即戦力になるExcelの教科書

入門者のExcel VBA―初めての人にベストな学び方 (ブルーバックス)

行列の対角化

　行列Aは相似変換によって対角行列にすることができる。Aが対称行列ならば、各固有値に対応する固有ベクトルは直交する。
　対角行列：対角成分以外の成分が0の行列
　対称行列：p行q列とq行p列の成分が一致する行列

　ある正則な行列Sによって、下記の行列を作る。

　　D = S^-1AS　　・・・(2)

　この変換を相似変換といい、相似変換により固有値や固有ベクトルは変化しない。

　Aの固有値をλ1,λ2,・・・,λnとし、対応する固有ベクトルをx1,x2,・・・,xnとすると、次式が成立する。

　　Axi = λixi　(i=1,2,・・・,n)　・・・(3)

　ここで、固有ベクトルを使いSを作る。

　　S = [x1,x2,・・・,xn]　　・・・(4)

　この時、Dは対角行列となり、その対角成分は固有値となる。

　Aが対称行列の場合、固有ベクトルxiの大きさを1(規格化)にすると、Sは直行行列となり、その逆行列は転置行列となる。

　　S^-1 = S^t　　・・・(5)

　よって、このSを求めることのより、固有値と固有ベクトルが算出できる。

ヤコビ法を用いた固有値、固有値計算

　ヤコビ法は実対称行列Aの固有値と固有ベクトルを求める方法である。

①　マトリックスAの非対角成分で絶対値が最大の成分(apq)を探す。

②　成分(apq)をピボットにし、面内回転を行う。

　回転角をθとすると、次式となる。

　tan(2θ) = 2apq/(app - aqq)

　θ= atan2(app-aqq, 2apq)/2 　　・・・ (6)

　注) Excel関数　atan(y/x)を用いると、x=0のときエラーとなる。それを避けるため、Excel関数 atan2(x,y)を使用する。

　回転マトリックスS1

　単位マトリックスのｐ,qの行列に下記回転成分を追加する。

　spp = cos(θ)
　spq = -sin(θ)
　sqp = sin(θ)
　sqq = cos(θ)

　回転を行う。

　A2 = S1^tAS1　　　・・・・　(7)

③　Aiに対し、①、②の処理をAiが対角行列になるまで、繰り返す。

　Ar+1 = Sr^tArSr 　　・・・・ (8)

　S = S1S2・・・Sr　　・・・・(9)

　これで、固有値と固有ベクトルが求まった。

Excelを用いたヤコビ法の１回目の計算

　図1にExcelを用いたヤコビ法の１回目の計算例題を示す。

　　図１　Excelを用いたヤコビ法の１回目の計算例題

ヤコビ法のVBA

　ヤコビ法の繰返し計算を行うためにVBAを作成した、そのソースコードを示す。

' ヤコビ法による固有値解析
Sub yacobi()
Dim i As Long
Dim n As Long
Dim j As Long
Dim k As Long
Dim im As Long
Dim jm As Long
Dim amax As Double
Dim q As Double
Dim a() As Double
Dim sr() As Double
n = Cells(3, 3).Value
ReDim a(n, n)
ReDim sr(n, n)
'　初期設定
For i = 1 To n
　For j = 1 To n
　　a(i, j) = Cells(i + 5, j + 1).Value
　Next
Next
For i = 1 To n
　For j = 1 To n
　　If i = j Then
　　　sr(i, j) = 1
　　Else
　　　sr(i, j) = 0
　　End If
　Next
Next
For i = 1 To n
　For j = 1 To n
　　Cells(i + 21, j + 1).Value = sr(i, j)
　Next
Next
'　繰返し
For k = 1 To 100
' 最大値の探索
　amax = 0#
　For i = 1 To n - 1
　　For j = i + 1 To n
　　　If Abs(a(i, j)) > amax Then
　　　　amax = Abs(a(i, j))
　　　　im = i
　　　　jm = j
　　　End If
　　Next
　Next
　If amax < 0.0000001 Then Exit For
　' Arの設定
　For i = 1 To n
　　For j = 1 To n
　　　Cells(i + 11, j + 1).Value = a(i, j)
　　Next
　Next
　' 最大値の取出し
　Cells(6, 6).Value = a(im, im)
　Cells(6, 7).Value = a(im, jm)
　Cells(7, 6).Value = a(jm, im)
　Cells(7, 7).Value = a(jm, jm)
　' Srの設定
　q = Cells(9, 7).Value
　For i = 1 To n
　　For j = 1 To n
　　　If i = j Then
　　　　sr(i, j) = 1
　　　Else
　　　　sr(i, j) = 0
　　　End If
　　Next
　Next
　sr(im, im) = Cos(q)
　sr(im, jm) = -Sin(q)
　sr(jm, im) = Sin(q)
　sr(jm, jm) = Cos(q)
　For i = 1 To n
　　For j = 1 To n
　　　Cells(i + 11, j + 5).Value = sr(i, j)
　　Next
　Next

　' Sのコピー
　For i = 1 To n
　　For j = 1 To n
　　　sr(i, j) = Cells(i + 21, j + 5).Value
　　Next
　Next
　For i = 1 To n
　　For j = 1 To n
　　　Cells(i + 21, j + 1).Value = sr(i, j)
　　Next
　Next
　' Ar+1の取出し
　For i = 1 To n
　　For j = 1 To n
　　　a(i, j) = Cells(i + 16, j + 1).Value
　　Next
　Next

Next
' 繰返し終了
End Sub

エクセルで操る！ヤコビ法による固有値、固有ベクトル計算

　
　ボタンにマクロを登録し、計算実行を行った。計算結果の対角マトリックスより、固有値は6.82843, 2, 1.17157となる。
固有ベクトルはマトリックスSの列となる。

　図２　ヤコビ法による固有値計算、固有ベクトル計算結果

エクセルで操る！ベアストウ法による固有値計算、固有ベクトル計算

エクセルで操る！べき乗法による固有値計算、固有ベクトル計算

EXCELを用いた科学技術計算

 EXCELのソルバーを用いた科学技術計算

〒６７３－００３６　兵庫県明石市松江６２－１４
(有)ゴッドフット企画