计算物理第三次作业new.docVIP

计算物理第三次作业 ——用不同方法求矩阵的本征值和本征向量 问题描述： 如上图所示的一个对称矩阵，我们要求取它的本值和本征向量有多种方法，此处，我们采用Jacobi Method（Givens）和Triple Diagonalization Method来求解相关的本征值和本征向量。这两种方法既有相似的地方也有不同之处：两种方法都采用一个正交的变换矩阵来不断转换原矩阵方式实现的，最终达到一个对角或三对角的矩阵；不同之处在于第一种直接将矩阵对角化，得到本征值及向量，第二种则先将矩阵三对角化，再利用转换方式对角化。因此两种方法中的变换方式有所不同，具体稍后详叙，此处不再赘述。 主要算法与公式： 1．Jacobi Method. 由于我们处理的是一个对称矩阵，这种方法的实现是通过找到一个矩阵Q，使其能将原矩阵变成一个非对角元为零的矩阵，此时对角元即为本征值，Q的每一个对应的列向量为本征向量，如下图所示为矩阵变换： 关键的问题归结到如何得到矩阵Q，如下： 首先，我们确定如下形式的矩阵（Givens）： 其中（p，q）对应变换后可以使矩阵中第p行第q个元素变成零，选取方法可对应该矩阵绝对值最大非对角元： 然后使矩阵按以下的方式作变换，最终成为只有对角元非零形式： 至此，变换后的矩阵本征值为D的对角元，Q的列向量为对应的本征向量。 2．Triple Diagonalization Method. 如上图，本方法的第一步是将矩阵三对角化，与第一种方法操作过程相似，不同之处在于矩阵Q形式稍有差别，我们采用Givens Transformation，即： （p，q）可理解变换后矩阵的（p-1，q）元素为零，取值取遍三对角以外所有值，矩阵Q中相关元素值： 经过若干次变换，就可以得到一个三对角矩阵及变换矩阵Q； 第二步是将三对角的矩阵变换为对角矩阵，此处我们用Givens方法，同第一种方法相似，稍有不同的是： 经过第二步，得到的对角矩阵对角元为本征值，将所有矩阵连乘得到的就是本征向量矩阵。 源代码： 1.Jacobi Method. #includeiostream #includeiomanip using namespace std; void main() { double A[5][5]={4,2,2,5,8,2,5,1,3,4,2,1,6,2,6,5,3,2,1,3,8,4,6,3,3}; double QS[5][5]={1,0,0,0,0,0,1,0,0,0, 0,0,1,0,0,0,0,0,1,0,0,0,0,0,1}; //运用 Givens Transformation将矩阵A对角化； int num=0; do { num+=1; int i,j; double d=0; for(int l=0;l4;l++) { for(int m=l+1;m5;m++) { if(abs(A[l][m])d) { d=abs(A[l][m]); i=l; j=m; } } }//寻找绝对值最大的非对角元； double s=(A[i][i]-A[j][j])/(2*A[i][j]),t=-pow(s*s+1,0.5)-s; double sino=t*pow(1+t*t,-0.5),coso=pow(1+t*t,-0.5); double QT[5][5]={1,0,0,0,0,0,1,0,0,0, 0,0,1,0,0,0,0,0,1,0,0,0,0,0,1}, Q[5][5]={1,0,0,0,0,0,1,0,0,0, 0,0,1,0,0,0,0,0,1,0,0,0,0,0,1}, b[5][5]; { Q[i][i]=coso,QT[i][i]=coso, Q[j][j]=coso,QT[j][j]=coso, Q[j][i]=sino,QT[j][i]=-sino, Q[i][j]=-sino,QT[i][j]=sino; }//赋值，计算矩阵Q及QT； { for(int l=0;l5;l++) { for(int m=0;m5;m++) { b[l][m]=0; for(int k=0;k5;k++) { b[l][m]+=QS[l][k]*Q[k][m]; } } } for(int l=0;l5;l++) { for(int m=0