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

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

一．实验目的

1．掌握定态薛定谔方程的矩阵解法。

2．掌握几种矩阵特征值问题数值解法的原理，会调用相应的子程序求解具体问题。

二．实验内容

1．问题描述

以/2ω/()m ω为长度单位，一维谐振子的哈密顿量为

2

202d H x dx

=-+, 其本征值为21n E n =+，本证波函数为

2

/2)()n n x H x ϕ=-， 其中()n H x 为厄米多项式，满足递推关系

11()2()2()n n n H x xH x nH x +-=-。

用矩阵方法求

2

22d H x x dx

=-++ 的本证能量和相应的波函数。

2．问题分析

H E ψψ=

0()|j j j t c ψϕ∞

==>∑

0||i i j i j i j c E c x Ec ϕϕ∞

=+<>=∑

11|j j j x ϕϕϕ-+>=>>

11||/2,||(1)/2j j j j x j x j ϕϕϕϕ-+<>=

<>=+ 0010010

112111,211,11,1

n n n n n n n n n n n n E x c c x E x c c E x E x c c x E c c -------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⋅=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣

⎦

3．程序编写

子程序及调用方法见《FORTRAN 常用算法程序集（第二版）》第三章 徐士良，P97

4．实验要求

◆用恰当的算法求解以上实对称三对角矩阵的特征值问题。

◆取n=8，给出H 的全部特征值和相应的特征向量。

5．实验步骤

● 启动软件开发环境Microsoft Developer Studio 。

● 创建新工作区shiyan03。

● 创建新项目xm3。

● 创建源程序文件xm3.f90，编辑输入源程序文本。

● 编译、构建、运行、调试程序。

6．实验结果

程序设计：

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)

10FORMAT(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))

DO 30 K=1,N

H=Q(K,I+1)

Q(K,I+1)=S*Q(K,I)+E*H

Q(K,I)=E*Q(K,I)-S*H

30 CONTINUE

40 CONTINUE

C(J)=S*P

B(J)=E*P

IF (ABS(C(J)).GT.D) GOTO 15

END IF

B(J)=B(J)+F

50CONTINUE

DO 80 I=1,N

K=I

P=B(I)

IF (I+1.LE.N) THEN

J=I

60 J=J+1

IF (J.LE.N) THEN