函數格式:Math_Matrix_EigenValue(K1,n,LoopNumber,Errro,Ret,IsHess)
K1:要求特徵值的方陣
n:方陣K1的階數
LoopNumber:循環次數
Errro:誤差控制變量
Ret:返回的特徵值,Ret是是n*2的數組,第一列是實數部分,第2列為虛數部分
IsHess:K1是否已經是上Hessenberg矩陣
函數說明:函數成功返回True,失敗返回False
特別說明:
:將方陣轉化為上Hessenberg矩陣
源代碼:
<code>Public Function Math_Matrix_EigenValue(ByVal K11(,) As Double, ByVal n As Int16, ByVal LoopNumber As Int16, ByVal Errro As Int16, ByRef Ret(,) As Double, ByVal IsHess As Boolean) As Boolean 'ret裡是n*2的數組，第一列是實數部分，第2列為虛數部分
Dim i As Int16 = K11.Length / n
If n * n <> K11.Length Then '只有方陣才有特徵值
Return False
End If
Dim j As Int16
Dim k As Int16
Dim t As Int16
Dim m As Int16
Dim A(0, 0) As Double
ReDim Ret(n - 1, 1) 'u v
Dim erro As Double = Math.Pow(0.1, Errro)
Dim b As Double
Dim c As Double
Dim d As Double
Dim g As Double
Dim xy As Double
Dim p As Double
Dim q As Double
Dim r As Double
Dim x As Double
Dim s As Double
Dim e As Double
Dim f As Double
Dim z As Double
Dim y As Double
Dim loop1 As Int16 = LoopNumber
Dim K1(n - 1, n - 1) As Double
For m = 0 To n - 1
For t = 0 To n - 1
K1(m, t) = K11(m, t)
Next
Next
If IsHess Then
i -= 1
ReDim A(i, n - 1)
For j = 0 To i
For k = 0 To n - 1
A(j, k) = K1(j, k)
Next
Next
Else
Math_Matrix_Hessenberg(K1, n, A) '將方陣K1轉化成上Hessenberg矩陣A
End If
m = n
While m <> 0
t = m - 1
While t > 0
If Math.Abs(A(t, t - 1)) > erro * (Math.Abs(A(t - 1, t - 1)) + Math.Abs(A(t, t))) Then
t -= 1
Else
Exit While
End If
End While
If t = m - 1 Then
Ret(m - 1, 0) = A(m - 1, m - 1)
Ret(m - 1, 1) = 0
m -= 1
loop1 = LoopNumber
ElseIf t = m - 2 Then
b = -(A(m - 1, m - 1) + A(m - 2, m - 2))
c = A(m - 1, m - 1) * A(m - 2, m - 2) - A(m - 1, m - 2) * A(m - 2, m - 1)
d = b * b - 4 * c
y = Math.Pow(Math.Abs(d), 0.5)
If d > 0 Then
xy = 1
If b xy = -1
End If
Ret(m - 1, 0) = -(b + xy * y) / 2
Ret(m - 1, 1) = 0
Ret(m - 2, 0) = c / Ret(m - 1, 0)
Ret(m - 2, 1) = 0
Else
Ret(m - 1, 0) = -b / 2
Ret(m - 2, 0) = Ret(m - 1, 0)
Ret(m - 1, 1) = y / 2
Ret(m - 2, 1) = -Ret(m - 1, 1)
End If
m -= 2
loop1 = LoopNumber
Else
If loop1 Return False
End If
loop1 -= 1
j = t + 2
While j A(j, j - 2) = 0
j += 1
End While
j = t + 3
While j A(j, j - 3) = 0
j += 1
End While
k = t
While k If k <> t Then
p = A(k, k - 1)
q = A(k + 1, k - 1)
If k <> m - 2 Then
r = A(k + 2, k - 1)
Else
r = 0
End If
Else
b = A(m - 1, m - 1)
c = A(m - 2, m - 2)
x = b + c
y = c * b - A(m - 2, m - 1) * A(m - 1, m - 2)
p = A(t, t) * (A(t, t) - x) + A(t, t + 1) * A(t + 1, t) + y
q = A(t + 1, t) * (A(t, t) + A(t + 1, t + 1) - x)
r = A(t + 1, t) * A(t + 2, t + 1)
End If
If p <> 0 Or q <> 0 Or r <> 0 Then
If p xy = -1
Else
xy = 1
End If
s = xy * Math.Pow(p * p + q * q + r * r, 0.5)
If k <> t Then
A(k, k - 1) = -s
End If
e = -q / s
f = -r / s
x = -p / s
y = -x - f * r / (p + s)
g = e * r / (p + s)
z = -x - e * q / (p + s)
For j = k To m - 1
b = A(k, j)
c = A(k + 1, j)
p = x * b + e * c
q = e * b + y * c
r = f * b + g * c
If k <> m - 2 Then
b = A(k + 2, j)
p += f * b
q += g * b
r += z * b
A(k + 2, j) = r
End If
A(k + 1, j) = q
A(k, j) = p
Next
j = k + 3
If j >= m - 1 Then
j = m - 1
End If
For i = t To j
b = A(i, k)
c = A(i, k + 1)
p = x * b + e * c
q = e * b + y * c
r = f * b + g * c
If k <> m - 2 Then
b = A(i, k + 2)
p += f * b
q += g * b
r += z * b
A(i, k + 2) = r
End If
A(i, k + 1) = q
A(i, k) = p
Next
End If
k += 1
End While
End If
End While
Return True
End Function/<code>
分享到:
閱讀更多 愛碼廝小妖 的文章
