实对称矩阵特征值和特征向量求解_C#

/// <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; } }