用C语言求一n阶实对称矩阵的特征矩阵和特征值

Python011

用C语言求一n阶实对称矩阵的特征矩阵和特征值,第1张

阶数通过修改 #define N 3

#include <math.h>

#include <alloc.h>

#include "c\elmhes.c "

#include "c\hqr.c "

#define N 3

void main()

{

double b[N+1][N+1]={{0,0,0,0},{0,3,1,-1},{0,2,2,-1},{0,2,2,0}}/*这里矩阵无需对称,若对称,则特征值均为实数*/

double **a

double wr[N+1],wi[N+1]/*特征值的实部和虚部*/

int i,j

a=(double **) malloc((unsigned) (N+1)*sizeof(double*))

for(i=0i <=Ni++)

a[i]=b[i] /*赋值完毕*/

printf( "Original A:\n ")

for(i=1i <=Ni++)

{

for(j=1j <=Nj++)

printf( "%f ",a[i][j])

printf( "\n ")

}

elmhes(a,N) /*变换至上Hessenburg型,特征值不变*/

printf( "After Elmhes:\n ")

for(i=1i <=Ni++)

for(j=1j <(i-1)j++)

a[i][j]=0

for(i=1i <=Ni++)

{

for(j=1j <=Nj++)

printf( "%f ",a[i][j])

printf( "\n ")

}

hqr(a,N,wr,wi)/*用QR方法求上Hessenburg型矩阵的特征值*/

printf( "Eigenvalue:\n ")

for(i=1i <=Ni++)

printf( "%f + %f * i\n ",wr[i],wi[i])

用C++或者VB编程很烦人的，matlab中命令：[a,b]=eig(A)就是求解矩阵A的特征值和特征值对应的向量，他们分别会构成一个由特征值组成的对角矩阵b和一个由对应特征值的特征列向量组成的a矩阵。或者命令a=eig[A]就只有特征值组成的对角矩阵a，别去想用C++和VB之类的，这些软件用来求解矩阵和matlab相差太远了。我之前也想过编程解决，人家一个命令就能解决的问题何不取巧呢？

方法1：推导出det(aA-I)=0的解析式，这应该是个四次方程，因为只有4阶，不是很困难的，写出后就可以用方程求根的方法求解（如newton迭代法）

方法2：如果你是对角优势阵，也就是对角线上的值的绝对值，比同行所有其他元素的绝对值的和还大，可以通过局部旋转的方法把矩阵“能量”集中到对角线

这个是方法，你可以自己去写一下试试~