|
楼主 |
发表于 2020-6-2 15:01
|
显示全部楼层
[广告] Excel易用宝 - 提升Excel的操作效率 · Excel / WPS表格插件 ★ 免费下载 ★ ★ 使用帮助★
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 模块名:MatrixModule.bas
' 函数名:MHbergEigenv
' 功能: 用QR方法计算赫申伯格(Hessen Berg)矩阵的全部特征值
' 参数: n - Integer型变量,赫申伯格矩阵的阶数。
' dblA - Double型二维数组,体积为n x n。存放赫申伯格矩阵。
' dblUR - Double型一维数组,长度为n,存放n个特征值的实部。
' dblUI - Double型一维数组,长度为n,存放n个特征值的虚部。。
' eps - Double型变量。迭代过程中的控制精度参数。
' nMaxItNum - Integer。为求得一个特征值所允许的最大迭代次数。
' 返回值: Boolean型。False,失败无解;True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MHbergEigenv(n As Integer, dblA() As Double, dblUR() As Double, dblUI() As Double, eps As Double, nMaxItNum As Integer) As Boolean
' 局部变量
Dim i As Integer, j As Integer, k As Integer, l As Integer, m As Integer, it As Integer
Dim b As Double, c As Double, w As Double, g As Double, xy As Double, p As Double, q As Double
Dim r As Double, x As Double, s As Double, e As Double, f As Double, z As Double, y As Double
it = 0
m = n + 1
While (m <> 1)
l = m - 1
While ((l > 1) And (Abs(dblA(l, l - 1)) > eps * (Abs(dblA(l - 1, l - 1)) + Abs(dblA(l, l)))))
l = l - 1
Wend
If (l = m - 1) Then
dblUR(m - 1) = dblA(m - 1, m - 1)
dblUI(m - 1) = 0#
m = m - 1
it = 0
Else
If (l = m - 2) Then
b = -(dblA(m - 1, m - 1) + dblA(m - 2, m - 2))
c = dblA(m - 1, m - 1) * dblA(m - 2, m - 2) - dblA(m - 1, m - 2) * dblA(m - 2, m - 1)
w = b * b - 4# * c
y = Sqr(Abs(w))
If (w > 0#) Then
xy = 1#
If (b < 0#) Then xy = -1#
dblUR(m - 1) = (-b - xy * y) / 2#
dblUR(m - 2) = c / dblUR(m - 1)
dblUI(m - 1) = 0#
dblUI(m - 2) = 0#
Else
dblUR(m - 1) = -b / 2#
dblUR(m - 2) = dblUR(m - 1)
dblUI(m - 1) = y / 2#
dblUI(m - 2) = -dblUI(m - 1)
End If
m = m - 2
it = 0
Else
If (it >= nMaxItNum) Then
MHbergEigenv = False
Exit Function
End If
it = it + 1
For j = l + 2 To m - 1
dblA(j, j - 2) = 0#
Next j
For j = l + 3 To m - 1
dblA(j, j - 3) = 0#
Next j
For k = l To m - 2
If (k <> l) Then
p = dblA(k, k - 1)
q = dblA(k + 1, k - 1)
r = 0#
If (k <> m - 2) Then r = dblA(k + 2, k - 1)
Else
x = dblA(m - 1, m - 1) + dblA(m - 2, m - 2)
y = dblA(m - 2, m - 2) * dblA(m - 1, m - 1) - dblA(m - 2, m - 1) * dblA(m - 1, m - 2)
p = dblA(l, l) * (dblA(l, l) - x) + dblA(l, l + 1) * dblA(l + 1, l) + y
q = dblA(l + 1, l) * (dblA(l, l) + dblA(l + 1, l + 1) - x)
r = dblA(l + 1, l) * dblA(l + 2, l + 1)
End If
If ((Abs(p) + Abs(q) + Abs(r)) <> 0#) Then
xy = 1#
If (p < 0#) Then xy = -1#
s = xy * Sqr(p * p + q * q + r * r)
If (k <> l) Then dblA(k, k - 1) = -s
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
p = x * dblA(k, j) + e * dblA(k + 1, j)
q = e * dblA(k, j) + y * dblA(k + 1, j)
r = f * dblA(k, j) + g * dblA(k + 1, j)
If (k <> m - 2) Then
p = p + f * dblA(k + 2, j)
q = q + g * dblA(k + 2, j)
r = r + z * dblA(k + 2, j)
dblA(k + 2, j) = r
End If
dblA(k + 1, j) = q
dblA(k, j) = p
Next j
j = k + 3
If (j >= m - 1) Then j = m - 1
For i = l To j
p = x * dblA(i, k) + e * dblA(i, k + 1)
q = e * dblA(i, k) + y * dblA(i, k + 1)
r = f * dblA(i, k) + g * dblA(i, k + 1)
If (k <> m - 2) Then
p = p + f * dblA(i, k + 2)
q = q + g * dblA(i, k + 2)
r = r + z * dblA(i, k + 2)
dblA(i, k + 2) = r
End If
dblA(i, k + 1) = q
dblA(i, k) = p
Next i
End If
Next k
End If
End If
Wend
MHbergEigenv = True
End Function |
|