实验三定态薛定谔方程的矩阵解法.docVIP

实验三 定态薛定谔方程的矩阵解法 实验目的 掌握定态薛定谔方程的矩阵解法。 掌握几种矩阵特征值问题数值解法的原理，会调用相应的子程序求解具体问题。 实验内容 问题描述 以为能量单位，为长度单位，一维谐振子的哈密顿量为, 其本征值为，本证波函数为 ， 其中为厄米多项式，满足递推关系 。 用矩阵方法求 的本证能量和相应的波函数。 问题分析 程序编写 子程序及调用方法见《FORTRAN常用算法程序集（第二版）》第三章 徐士良，P97 实验要求 用恰当的算法求解以上实对称三对角矩阵的特征值问题。 取n=8，给出H的全部特征值和相应的特征向量。 实验步骤 启动软件开发环境Microsoft Developer Studio。 创建新工作区shiyan03。 创建新项目xm3。 创建源程序文件xm3.f90，编辑输入源程序文本。 编译、构建、运行、调试程序。 实验结果 程序设计： DIMENSION Q(9,9),B(9),C(9) DOUBLE PRECISION Q,B,C,I N=9 Q=0 DO I=1,N Q(I,I)=1 END DO DO I=1,N B(I)=2*(I-1)+1 ENDDO DO I=1,N-1 C(I)=SQRT(I/2.) ENDDO ESP=0.000001 CALL CSSTQ(N,B,C,Q,EPS,L) !WRITE(*,*) ! WRITE(*,10) 10 FORMAT(1X,MAT A IS:) !WRITE(*,50) ((A(I,J),J=1,N),I=1,N) IF (L.NE.0) THEN !WRITE(*,*) ! WRITE(*,30) 30 FORMAT(1X,MAT Q IS:) WRITE(*,50) ((Q(I,J),J=1,N),I=1,N) WRITE(*,*) ! WRITE(*,40) 40 FORMAT(1X,MAT B IS:) WRITE(*,50) (B(I),I=1,N) 50 FORMAT(1X,9F8.3) END IF WRITE(*,*) END SUBROUTINE CSSTQ(N,B,C,Q,EPS,L) DIMENSION B(N),C(N),Q(N,N) DOUBLE PRECISION B,C,Q,D,H,P,R,F,E,S,G C(N)=0.0 D=0.0 F=0.0 DO 50 J=1,N IT=0 H=EPS*(ABS(B(J))+ABS(C(J))) IF (H.GT.D) D=H M=J-1 10 M=M+1 IF (M.LE.N) THEN IF (ABS(C(M)).GT.D) GOTO 10 END IF IF (M.NE.J) THEN 15 IF (IT.EQ.60) THEN L=0 WRITE(*,18) 18 FORMAT(1X, FAIL) RETURN END IF IT=IT+1 G=B(J) P=(B(J+1)-G)/(2.0*C(J)) R=SQRT(P*P+1.0) IF (P.GE.0.0) THEN B(J)=C(J)/(P+R) ELSE B(J)=C(J)/(P-R) END IF H=G-B(J) DO 20 I=J+1,N 20 B(I)=B(I)-H F=F+H P=B(M) E=1.0 S=0.0 DO 40 I=M-1,J,-1 G=E*C(I) H=E*P IF (ABS(P).GE.ABS(C(I))) THEN E=C(I)/P R=SQRT(E*E+1.0) C(I+1)=S*P*R S=E/R E=1.0/R ELSE E=P/C(I) R=SQRT(E*E+1.0) C(I+1)=S*C(I)*R S=1.0/R E=E/R END IF P=E*B(I)-S*G B(I+1)=H+S*(E*G+S*B(I)) D