Option Explicit On Option Strict On Public Class Form1 Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click Dim i As Integer Dim j As Integer Dim k As Integer Dim nf As Integer Dim a(,) As Double Dim b(,) As Double Dim ev() As Double Dim vec(,) As Double Dim dat As String Dim ar(,) As Double Dim br(,) As Double Dim c(,) As Double Dim d(,) As Double nf = 3 ReDim a(nf, nf) ReDim b(nf, nf) ReDim ev(nf) ReDim vec(nf, nf) ReDim ar(nf, nf) ReDim br(nf, nf) ReDim c(nf, nf) ReDim d(nf, nf) 'Case-1 a(1, 1) = 4 : a(1, 2) = 2 : a(1, 3) = 3 a(2, 1) = 2 : a(2, 2) = 6 : a(2, 3) = 1 a(3, 1) = 3 : a(3, 2) = 1 : a(3, 3) = 5 b(1, 1) = 1 : b(1, 2) = 0 : b(1, 3) = 0 b(2, 1) = 0 : b(2, 2) = 1 : b(2, 3) = 0 b(3, 1) = 0 : b(3, 2) = 0 : b(3, 3) = 1 'a(1, 1) = 106.482 : a(1, 2) = -59.228 : a(1, 3) = 10.816 'a(2, 1) = -59.228 : a(2, 2) = 89.1 : a(2, 3) = -42.618 'a(3, 1) = 10.816 : a(3, 2) = -42.618 : a(3, 3) = 33.348 'b(1, 1) = 0.01 : b(1, 2) = 0.0 : b(1, 3) = 0.0 'b(2, 1) = 0.0 : b(2, 2) = 0.01 : b(2, 3) = 0.0 'b(3, 1) = 0.0 : b(3, 2) = 0.0 : b(3, 3) = 0.01 'a(1, 1) = 3 : a(1, 2) = 2 : a(1, 3) = -2 'a(2, 1) = 2 : a(2, 2) = 3 : a(2, 3) = -2 'a(3, 1) = -2 : a(3, 2) = -2 : a(3, 3) = 8 'b(1, 1) = 1 : b(1, 2) = 0 : b(1, 3) = 0 'b(2, 1) = 0 : b(2, 2) = 1 : b(2, 3) = 0 'b(3, 1) = 0 : b(3, 2) = 0 : b(3, 3) = 1 'a(1, 1) = 5 : a(1, 2) = 4 : a(1, 3) = 3 : a(1, 4) = 2 : a(1, 5) = 1 'a(2, 1) = 4 : a(2, 2) = 4 : a(2, 3) = 3 : a(2, 4) = 2 : a(2, 5) = 1 'a(3, 1) = 3 : a(3, 2) = 3 : a(3, 3) = 3 : a(3, 4) = 2 : a(3, 5) = 1 'a(4, 1) = 2 : a(4, 2) = 2 : a(4, 3) = 2 : a(4, 4) = 2 : a(4, 5) = 1 'a(5, 1) = 1 : a(5, 2) = 1 : a(5, 3) = 1 : a(5, 4) = 1 : a(5, 5) = 1 'b(1, 1) = 1 : b(1, 2) = 0 : b(1, 3) = 0 : b(1, 4) = 0 : b(1, 5) = 0 'b(2, 1) = 0 : b(2, 2) = 1 : b(2, 3) = 0 : b(2, 4) = 0 : b(2, 5) = 0 'b(3, 1) = 0 : b(3, 2) = 0 : b(3, 3) = 1 : b(3, 4) = 0 : b(3, 5) = 0 'b(4, 1) = 0 : b(4, 2) = 0 : b(4, 3) = 0 : b(4, 4) = 1 : b(4, 5) = 0 'b(5, 1) = 0 : b(5, 2) = 0 : b(5, 3) = 0 : b(5, 4) = 0 : b(5, 5) = 1 For i = 1 To nf For j = 1 To nf ar(i, j) = a(i, j) br(i, j) = b(i, j) Next j Next i Call GJACOBI(nf, a, b, ev, vec) Console.WriteLine("*** Eigen value ***") dat = String.Format("{0,16:0.000000e+000}", ev(1)) For j = 2 To nf dat = dat & String.Format("{0,16:0.000000e+000}", ev(j)) Next j Console.WriteLine(dat) Console.WriteLine("*** Eigen vector ***") For i = 1 To nf dat = String.Format("{0,16:0.000000e+000}", vec(i, 1)) For j = 2 To nf dat = dat & String.Format("{0,16:0.000000e+000}", vec(i, j)) Next j Console.WriteLine(dat) Next i 'Case-1 出力結果 '*** Eigen value *** '9.000000e+000 4.732051e+000 1.267949e+000 '*** Eigen vector *** '9.227825e-001 -2.633891e-001 1.060632e+000 '9.227825e-001 9.829816e-001 -2.841954e-001 '9.227825e-001 -7.195925e-001 -7.764364e-001 '検算 For i = 1 To nf For j = 1 To nf c(i, j) = 0.0 d(i, j) = 0.0 For k = 1 To nf c(i, j) = c(i, j) + ar(i, k) * vec(k, j) d(i, j) = d(i, j) + br(i, k) * vec(k, j) * ev(j) Next k Next j Next i Console.WriteLine("*** 検算 ***") For i = 1 To nf dat = String.Format("{0,16:0.000000e+000}", c(i, 1)) For j = 2 To nf dat = dat & String.Format("{0,16:0.000000e+000}", c(i, j)) Next j dat = dat & " " dat = dat & String.Format("{0,16:0.000000e+000}", d(i, 1)) For j = 2 To nf dat = dat & String.Format("{0,16:0.000000e+000}", d(i, j)) Next j Console.WriteLine(dat) Next i End Sub Private Sub GJACOBI(ByVal nf As Integer, ByRef a(,) As Double, ByRef b(,) As Double, ByRef ev() As Double, ByRef vec(,) As Double) '************************************ '一般化JACOBI法による固有値計算 '(パソコン構造力学，宮澤健二著) '************************************ Dim i As Integer Dim j As Integer Dim k As Integer Dim fmk As Double Dim fmg As Double Dim fact As Double = 1.0E-20 Dim epsk As Double Dim epsg As Double Dim kmax As Integer Dim ip As Integer Dim iq As Integer Dim amax As Double Dim dam As Double Dim de1 As Double Dim de2 As Double Dim de3 As Double Dim dd As Double Dim ss As Double Dim x As Double Dim gg As Double Dim ww As Double Dim im As Integer fmk = 0.0 fmg = 0.0 For i = 1 To nf For j = 1 To nf If Math.Abs(a(i, j)) > fmk Then fmk = Math.Abs(a(i, j)) If Math.Abs(b(i, j)) > fmg Then fmg = Math.Abs(b(i, j)) Next j Next i epsk = fmk * fact epsg = fmg * fact For i = 1 To nf For j = 1 To nf vec(i, j) = 0.0 Next j vec(i, i) = 1.0 Next i kmax = 5 * nf * nf For k = 1 To kmax ip = 1 iq = 2 amax = a(1, 2) ^ 2 / (a(1, 1) ^ 2 + a(2, 2) ^ 2) For i = 1 To nf - 1 For j = i + 1 To nf dam = a(i, j) ^ 2 / (a(i, i) ^ 2 + a(j, j) ^ 2) If amax < dam Then ip = i : iq = j : amax = dam End If Next j Next i For i = 1 To nf - 1 For j = i + 1 To nf dam = b(i, j) ^ 2 / (b(i, i) ^ 2 + b(j, j) ^ 2) If amax < dam Then ip = i : iq = j : amax = dam End If Next j Next i de1 = a(ip, ip) * b(ip, iq) - a(ip, iq) * b(ip, ip) de2 = a(iq, iq) * b(ip, iq) - a(ip, iq) * b(iq, iq) de3 = a(ip, ip) * b(iq, iq) - a(iq, iq) * b(ip, ip) dd = de3 * de3 + 4.0 * de1 * de2 ss = Math.Sqrt(dd) If de3 < 0 Then x = 0.5 * (de3 - ss) Else x = 0.5 * (de3 + ss) End If If x = 0.0 Then dam = 0.0 : gg = -a(ip, iq) / a(iq, iq) Else dam = de2 / x : gg = -de1 / x End If For j = 1 To nf ww = a(ip, j) a(ip, j) = a(ip, j) + gg * a(iq, j) a(iq, j) = a(iq, j) + dam * ww ww = b(ip, j) b(ip, j) = b(ip, j) + gg * b(iq, j) b(iq, j) = b(iq, j) + dam * ww Next j For i = 1 To nf ww = a(i, ip) a(i, ip) = a(i, ip) + gg * a(i, iq) a(i, iq) = a(i, iq) + dam * ww ww = b(i, ip) b(i, ip) = b(i, ip) + gg * b(i, iq) b(i, iq) = b(i, iq) + dam * ww ww = vec(i, ip) vec(i, ip) = vec(i, ip) + gg * vec(i, iq) vec(i, iq) = vec(i, iq) + dam * ww Next i If amax < epsk And amax < epsg Then Exit For Next k For i = 1 To nf ev(i) = a(i, i) / b(i, i) Next i '並び替え 'ev(k)：最大固有値k=1，最小固有値k=nf 'vec(i,k)：k番目の固有値に対応する固有ベクトル For k = 1 To nf - 1 amax = -1.0E+30 For i = k To nf If ev(i) >= amax Then im = i : amax = ev(i) End If Next i ww = ev(k) : ev(k) = ev(im) : ev(im) = ww For i = 1 To nf ww = vec(i, k) vec(i, k) = vec(i, im) vec(i, im) = ww Next i Next k End Sub End Class