一、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+1 sgn(β)(2)
我们可以得到： cos ⁡ ( θ ) = 1 t 2 + 1 ( 3 ) \cos(\theta)=\frac{1}{\sqrt{t^2+1}}\qquad (3) cos(θ)=t2+1 1(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)

Created with Raphaël 2.2.0 开始矩阵A相似对角化主对角元素平方和不变结束 yes no

  该矩阵对角元素即为该是实对称阵的特征值

二、VB程序

代码如下（示例）：

1.以下为预定义部分：

Public Function Jacobi(ByVal A(,) As Double) As Double()Dim n As Double = A.GetUpperBound(0)Dim max As DoubleDim P As IntegerDim Q As IntegerDim U(n, n) As DoubleDim fail As DoubleConst PI As Double = 3.1415926Dim t As DoubleDim c As DoubleDim result(n) As Double

2.选取非对角线元素的主元素（偏离主对角最大的元素）

Domax = Math.Abs(A(0, 1))For I = 0 To nFor J = I To nIf I <> J ThenIf Math.Abs(A(I, J)) >= max Thenmax = Math.Abs(A(I, J))P = IQ = JEnd IfEnd IfNextNext

3.构造旋转矩阵（通过上述分析两种情况）

For I = 0 To nU(I, I) = 1NextIf A(P, P) = A(Q, Q) Thenfail = 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）Elsec = (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.5U(P, Q) = -t / (1 + t ^ 2) ^ 0.5U(Q, P) = -U(P, Q)U(Q, Q） = U(P, P)End If

4.迭代直至收敛

 A = Similaritydiagon_alization(A, U)For I = 0 To nFor J = 0 To nU（I, J) = 0NextNextmax = Math.Abs(A(0, 1))For I = 0 To nFor J = I To nIf I <> J ThenIf Math.Abs(A(I, J)) >= max Thenmax = Math.Abs(A(I, J))P = IQ = JEnd IfEnd IfNextNextFor I = 0 To nresult(I) = A(I, I)NextLoop Until Math.Abs(A(P, Q)) < 0.00001Return resultEnd 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 DoubleDim temp(,) As DoubleUT = Transposition(U)temp = Multiple(UT, A)Return Multiple(temp, U)End Function

Function sgn(ByVal a As Double) As DoubleIf a > 0 ThenReturn 1ElseIf a < 0 ThenReturn -1ElseReturn 0End IfEnd Function

jacobi旋转法的VB实现相关推荐

利用特征值与特征向量求解弹性力学中的主应力与主平面问题
利用特征值与特征向量求解弹性力学中的主应力与主平面问题前言一.二向应力状态 1. 莫尔圆图解法 2. 特征值与特征向量解法二.三向应力状态前言已知物体在任意一点的六个应力分量(σx,σy,σ ...
ktt算法约化_矩阵特征与特征向量的计算
矩阵特征与特征向量的计算第三章第三章矩阵特征与特征向量的计算矩阵特征与特征向量的计算3.1 引言引言在科学技术的应用领域中,许多问题都归为求解一个特征系统.如动力学系统和结构系统中的振动问题,求 ...
matlab幂法与反幂法,matlab位移反幂法
本次实验所用的幂法和反幂法分别是求解最大特征值和最小特征值,并根据它们的结果求解二条件数.幂法和反幂法的 Matlab 程序很好的解决了手算时所会遇到的...... 类似幂法和反幂法可以写出按模最小 ...
vb中可视对象的操作
问题 : 在调试机房结账的部分,这两部分总是出问题,实时错误424. 错误解释: 未找到窗体(错误 424) 后来通过大量的查阅,找到了答案. MSHFlexGrid1是一个"控件" ...
机房收费系统【VB版】——前期准备
前言: 没有源码和参考的机房收费系统,很犯怵的开始,完全不懂如何下手,经过后来小伙伴的交流和巨人的博客. 准备: 1.安装机房收费系统程序 1.1添加ODBC数据源--添加文件DSN--附加数据库-- ...
【VB】学生信息管理系统6——错误调试
因为站在了巨人的肩膀上,在理解代码意思后的调试中,用到之前的别人的CSDN.所以原理查的不是很透彻.这里总结一下我的问题! 1.VB(如下代码)中mrc.EOF = False应该怎么理解呢? Set ...
【VB】学生信息管理系统5——数据库代码
这次学生信息管理系统在代码的理解过程中遇到了一些问题.总结如下: 1. sql server的安装过程各个步骤的意思.在安装SQL Server的时候按照网上的步骤,我觉得这个需要学完整个数据库再返回 ...
【VB】学生信息管理系统4——数据库的发展
由于连接数据的时候出现了很多不懂得问题,为什么要连接,它是怎么连接的,查着查着,就越看越多.又不舍得就这么放过这些问题,所以就耐心看看究竟是怎么回事! 1.自从出现数据库,人们渴望用数据和应用程序做交 ...
【VB】学生信息管理系统2——窗体设计
这次学生系统是照着书敲的,先敲完然后开始调试!中途遇到了很多问题,查了很多,这里不容易系统的总结!所以就针对各个问题,各个击破! 问题一:VB 6.0中,状态栏控件(sbstatusbar):右击选项 ...

jacobi旋转法的VB实现