求实对称阵的 特征值 和 特征向量(转)

Posted 2020-06-25 撑起一片天

 /// <summary>
        /// 求实对称阵的 特征值 和 特征向量
        /// </summary>
        /// <param name="data">实对称阵</param>
        /// <param name="num">维数</param>
        /// <param name="eigenvalue">引用参数 特征值 回传</param>
        /// <param name="eigenvector">引用参数 特征向量 回传</param>
        /// <returns>是否成功</returns>
        private bool GetEigenValueAndEigenVector(double[,] data, int num,ref double [] eigenvalue,ref double [,] eigenvector)
        {
            try
            {
                double[,] A = data;

                //E 单位标准矩阵   存储 特征向量--------------------------------------------
                double[,] V = new double[num, num];
                for (int iv = 0; iv < num; iv++)
                {
                    for (int iv2 = 0; iv2 < num; iv2++)
                    {
                        if (iv == iv2)
                        {
                            V[iv, iv2] = 2;
                        }
                        else
                        {
                            V[iv, iv2] = 2;
                        }
                    }
                }
                //----------------------------------------------


                double[] eigsv = new double[num];//存储 特征值
                for (int ieigsv = 0; ieigsv < num; ieigsv++)
                {
                    eigsv[ieigsv] = 0;
                }


                double epsl = 0.0001;
                int maxt = 10;
                int n = num;


                double tao, t, cn, sn; // 临时变量
                double maxa;// 记录非对角线元素最大值

                //------------------------------------------------------------------------------------------------
                for (int it = 0; it < maxt; it++)
                {
                    maxa = 0;
                    for (int p = 0; p < n - 1; p++)
                    {
                        for (int q = p + 1; q < n; q++)
                        {
                            if (Math.Abs(A[p, q]) > maxa) // 记录非对角线元素最大值
                            {
                                maxa = Math.Abs(A[p, q]);
                            }
                            if (Math.Abs(A[p, q]) > epsl) // 非对角线元素非0时才执行Jacobi变换
                            {
                                // 计算Givens旋转矩阵的重要元素:cos(theta), 即cn, sin(theta), 即sn
                                tao = 0.5 * (A[q, q] - A[p, p]) / A[p, q];

                                if (tao >= 0) // t为方程 t^2 + 2*t*tao - 1 = 0的根, 取绝对值较小的根为t
                                {
                                    t = -tao + Math.Sqrt(1 + tao * tao);
                                }
                                else
                                {
                                    t = -tao - Math.Sqrt(1 + tao * tao);
                                }
                                cn = 1 / Math.Sqrt(1 + t * t);
                                sn = t * cn;

                                // Givens旋转矩阵之转置左乘A, 即更新A的p行和q行
                                for (int j = 0; j < n; j++)
                                {
                                    double apj = A[p, j];
                                    double aqj = A[q, j];
                                    A[p, j] = cn * apj - sn * aqj;
                                    A[q, j] = sn * apj + cn * aqj;
                                }

                                // Givens旋转矩阵右乘A, 即更新A的p列和q列
                                for (int i = 0; i < n; i++)
                                {
                                    double aip = A[i, p];
                                    double aiq = A[i, q];
                                    A[i, p] = cn * aip - sn * aiq;
                                    A[i, q] = sn * aip + cn * aiq;
                                }

                                // 更新特征向量存储矩阵V, V=J0×J1×J2...×Jit, 所以每次只更新V的p, q两列
                                for (int i2 = 0; i2 < n; i2++)
                                {
                                    double vip = V[i2, p];
                                    double viq = V[i2, q];
                                    V[i2, p] = cn * vip - sn * viq;
                                    V[i2, q] = sn * vip + cn * viq;
                                }
                            }

                        }

                    }




                    if (maxa < epsl) // 非对角线元素已小于收敛标准，迭代结束
                    {

                        break;
                    }
                }
                //-----------------------------------------------------------------------------------------------------

                // 特征值向量  排序
                for (int j2 = 0; j2 < n; j2++)
                {
                    eigsv[j2] = A[j2, j2];
                    //fprintf(fp2, "%f ",eigsv[j2]);
                }
                // 对特征值向量从大到小进行排序, 并调整特征向量顺序 (直接插入法)
                double[] tmp = new double[n];
                for (int j = 1; j < n; j++)
                {
                    int i = j;
                    double a = eigsv[j];
                    for (int k = 0; k < n; k++)
                    {
                        tmp[k] = V[k, j];
                    }
                    while (i > 0 && eigsv[i - 1] < a)
                    {
                        eigsv[i] = eigsv[i - 1];
                        for (int k = 0; k < n; k++)
                        {
                            V[k, i] = V[k, i - 1];
                        }
                        i--;
                    }
                    eigsv[i] = a;
                    for (int k2 = 0; k2 < n; k2++)
                    {
                        V[k2, i] = tmp[k2];
                    }
                }
                //----------------------------------------------------------------------------------------------------------
                //打印特征值 与 特征向量    数据回传
                for (int ivc = 0; ivc < num; ivc++)
                {
                    for (int ivc2 = 0; ivc2 < num; ivc2++)
                    {

                        //fprintf(fp2, "%f%s", V[ivc, ivc2], "\n");
                        MessageBox.Show(V[ivc, ivc2].ToString());
                        eigenvector[ivc, ivc2] = V[ivc, ivc2];

                    }

                }


                for (int ieigsvc = 0; ieigsvc < num; ieigsvc++)
                {

                    //fprintf(fp4, "%f%s", eigsv[ieigsvc], "\n");
                    MessageBox.Show(eigsv[ieigsvc].ToString());
                    eigenvalue[ieigsvc] = eigsv[ieigsvc];

                }
                //----------------------------------------------------------------------------------------------------------
                return true;
            }
            catch
            {
                return false;
            }
        }