jacobi旋转法的VB实现
一、jacobi方法概要
1.jacobi方法的基本思想:
jacobi旋转法是通过一组平面旋转变换(构建一个平面旋转矩阵)对实对称方阵A进行相似对角化,化其为对角阵,进而求出特征值与特征向量的方法。
2.jacobi方法的基本步骤:
<1>.选取非对角线元素的主元素(这里只画出了上三角阵)
a
i
(
p
,
q
)
a_i(p,q)
ai(p,q)
(
a
11
a
12
a
13
⋯
a
1
n
0
a
22
a
23
⋯
a
2
n
⋮
⋮
⋮
a
p
q
⋮
0
0
0
⋯
a
n
n
)
\begin{pmatrix} a_{11}& a_{12 }& a_{13} & \cdots & a_{1n }\\ 0 & a_{22} & a_{23} & \cdots & a_{2n} \\ \vdots & \vdots & \vdots &\color{#00F}{a_{pq}} & \vdots \\ 0 & 0 & 0 & \cdots & a_{nn}\\ \end{pmatrix}
⎝⎜⎜⎜⎛a110⋮0a12a22⋮0a13a23⋮0⋯⋯apq⋯a1na2n⋮ann⎠⎟⎟⎟⎞
<2> 构造如下形式的旋转矩阵
其中:(1)
a
p
p
=
a
q
q
a_{pp}=a_{qq}
app=aqq时,
tan
(
2
θ
)
=
∞
\tan(2\theta)=\infty
tan(2θ)=∞于是:
θ
=
π
4
s
g
n
(
a
p
q
)
(
1
)
\theta=\frac \pi 4 sgn(a_{pq})\qquad (1)
θ=4πsgn(apq)(1)
(2)
a
p
p
≠
a
q
q
a_{pp}\neq a_{qq}
app=aqq时,
t
=
s
g
n
(
β
)
∥
β
∥
+
β
2
+
1
(
2
)
t=\frac{sgn(\beta)}{\lVert \beta \lVert +\sqrt{\beta^2+1}}\qquad (2)
t=∥β∥+β2+1sgn(β)(2)
我们可以得到:
cos
(
θ
)
=
1
t
2
+
1
(
3
)
\cos(\theta)=\frac{1}{\sqrt{t^2+1}}\qquad (3)
cos(θ)=t2+11(3)
sin
(
θ
)
=
c
∗
t
(
4
)
\sin(\theta)=c*t\qquad (4)
sin(θ)=c∗t(4)
<3> 构造如下形式的迭代 A k = R T A k − 1 R ( 5 ) A_k=R^TA_{k-1}R\qquad (5) Ak=RTAk−1R(5)
<4> … … \ldots\ldots ……
判断收敛的条件为:
∑ i = 1 n a i 2 = c ( 6 ) \sum_{i=1}^n a_i^2\quad =c\qquad (6) i=1∑nai2=c(6)
该矩阵对角元素即为该是实对称阵的特征值
二、VB程序
代码如下(示例):
1.以下为预定义部分:
Public Function Jacobi(ByVal A(,) As Double) As Double()
Dim n As Double = A.GetUpperBound(0)
Dim max As Double
Dim P As Integer
Dim Q As Integer
Dim U(n, n) As Double
Dim fail As Double
Const PI As Double = 3.1415926
Dim t As Double
Dim c As Double
Dim result(n) As Double
2.选取非对角线元素的主元素(偏离主对角最大的元素)
Do
max = Math.Abs(A(0, 1))
For I = 0 To n
For J = I To n
If I <> J Then
If Math.Abs(A(I, J)) >= max Then
max = Math.Abs(A(I, J))
P = I
Q = J
End If
End If
Next
Next
3.构造旋转矩阵(通过上述分析两种情况)
For I = 0 To n
U(I, I) = 1
Next
If A(P, P) = A(Q, Q) Then
fail = 0.25 * PI * sgn(A(P, Q))
U(P, P) = Math.Cos(fail)
U(P, Q) = -Math.Sin(fail)
U(Q, P) = Math.Sin(fail)
U(Q, Q) = Math.Cos(fail)
Else
c = (A(P, P) - A(Q, Q)) / (2 * A(P, Q))
t = sgn(c) / (Math.Abs(c) + (c ^ 2 + 1) ^ 0.5)
U(P, P) = 1 / (1 + t ^ 2) ^ 0.5
U(P, Q) = -t / (1 + t ^ 2) ^ 0.5
U(Q, P) = -U(P, Q)
U(Q, Q) = U(P, P)
End If
4.迭代直至收敛
A = Similaritydiagon_alization(A, U)
For I = 0 To n
For J = 0 To n
U(I, J) = 0
Next
Next
max = Math.Abs(A(0, 1))
For I = 0 To n
For J = I To n
If I <> J Then
If Math.Abs(A(I, J)) >= max Then
max = Math.Abs(A(I, J))
P = I
Q = J
End If
End If
Next
Next
For I = 0 To n
result(I) = A(I, I)
Next
Loop Until Math.Abs(A(P, Q)) < 0.00001
Return result
End Function
另外补充了一个子函数, S i m i l a r i t y d i a g o n a l i z a t i o n Similaritydiagon alization\quad Similaritydiagonalization和 s g n sgn sgn
'Jacobi相似对角算子’
Public Function Similaritydiagon_alization(ByVal A(,) As Double, ByVal U(,) As Double) As Double(,)
Dim n As Double = A.GetUpperBound(0)
Dim UT(,) As Double
Dim temp(,) As Double
UT = Transposition(U)
temp = Multiple(UT, A)
Return Multiple(temp, U)
End Function
Function sgn(ByVal a As Double) As Double
If a > 0 Then
Return 1
ElseIf a < 0 Then
Return -1
Else
Return 0
End If
End Function