jacobi旋转法的VB实现

一、jacobi方法概要

1.jacobi方法的基本思想:

jacobi旋转法是通过一组平面旋转变换（构建一个平面旋转矩阵）对实对称方阵A进行相似对角化，化其为对角阵，进而求出特征值与特征向量的方法。

2.jacobi方法的基本步骤:

<1>.选取非对角线元素的主元素(这里只画出了上三角阵）$a_i(p,q)$
$$\left(\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \cdots & a_{1n}\\ 0& a_{22}& a_{23}& \cdots & a_{2n}\\ ⋮& ⋮& ⋮& a_{pq}& ⋮\\ 0& 0& 0& \cdots & a_{nn}\end{array}\right)$$

<2> 构造如下形式的旋转矩阵

其中：(1) $a_{pp}=a_{qq}$时，$$\tan (2\theta )=\infty$$于是：
$$\theta =\frac{\pi }{4}sgn(a_{pq})(1)$$
(2)$a_{pp}\ne a_{qq}$时, $$t=\frac{sgn(\beta )}{\Vert \beta \Vert +\sqrt{{\beta }^2+1}}(2)$$
我们可以得到：$$\cos (\theta )=\frac{1}{\sqrt{t^2+1}}(3)$$
$$\sin (\theta )=c\ast t(4)$$

<3> 构造如下形式的迭代$$A_k=R^T{A}_{k-1}R(5)$$

<4> $$\dots \dots$$

判断收敛的条件为：

$$\sum _{i=1}^n{a}_i^2=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

另外补充了一个子函数，$Similaritydiagonalization$和$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