四阶龙格库塔法解含时薛定谔方程 | 寒山烟雨
现在的位置: 首页 > 小试身手 > 正文

四阶龙格库塔法解含时薛定谔方程

2012年04月19日 小试身手 ⁄ 共 1987字 ⁄ 字号 评论 4 条 ⁄ 阅读 5,647 views 次

四阶龙格库塔法解含时薛定谔方程,一维线性谐振子附加含时驱动力。

 Fortran | 
 
 copy code |
?

001
external f
002
parameter(n=100001,m=11)
003
Dimension Y(m),D(m),z(m,n),b(m),c(m,n)
004
complex(8) y,d,z,b,h
005
real(8) t,a
006
t=0.0d0</p>
007
do j=1,m
008
Y(j)=(0.0,0.0)
009
enddo
010
Y(1)=(1.0d0,0.0d0)
011
 
012
h=2500.d0/(n-1)
013
call grkt1(t,y,m,h,n,z,f,d,b)
014
write(*,*)
015
call Psi(z)
016
end
017
 
018
subroutine Psi(z)
019
parameter(n=100001,m=11,pi=3.1415926d0)
020
complex(8) z(m,n),p(20,1000)
021
real(8) a(m,n)
022
integer t
023
 
024
open(200,File='test.txt')
025
open(201,File='test1.txt')
026
open(202,File='test2.txt')
027
open(203,File='test3.txt')
028
 
029
do k=1,20,1
030
do j=1,1000
031
x=-5.0d0+j*0.01d0
032
P(k,j)=(0.d0,0.d0)
033
do i=1,m
034
p(k,j)=p(k,j)+z(i,k*5000)/sqrt((sqrt(pi)*(2**i)*JC(i)))*exp(-x*x/2.0d0)*H(i,x*1.0)*exp((0.d0,-1.d0)*(2.d0*i+1.d0)*k*5000*(2500.d0/(n-1)))
035
!p(k,j)=p(k,j)+z(i,k)/sqrt((sqrt(pi)*(2**i)*JC(i)))*exp(-x*x/2.0d0)*H(i,x*1.0)*exp((0.d0,-1.d0)*(2.d0*i+1.d0)*k*(2500.d0/(n-1)))
036
enddo
037
enddo
038
enddo
039
 
040
do j=1,1000
041
write(200,*)(-5.d0+j*0.01d0),(abs(p(i,j))**2,i=1,2)
042
write(201,*)(abs(p(i,j))**2,i=3,5)
043
write(202,*)(abs(p(i,j))**2,i=6,8)
044
write(203,*)(abs(p(i,j))**2,i=9,11)
045
enddo
046
 
047
end subroutine
048
 
049
subroutine F(t,y,m,d)
050
 
051
Dimension Y(m),D(m)
052
complex(8) y,d
053
real(8) t
054
D(1)=-(0.d0,1.d0)*g(t)*(Y(2)*sqrt(1.d0/2.d0)*exp(-2.d0*(0,1)*t))
055
do j=2,m-1
056
D(j)=-(0,1)*g(t)*(Y(j-1)*sqrt((j-1)/2.d0)*exp(2.d0*(0,1)*t)+Y(j+1)*sqrt(j/2.)*exp(-2.d0*(0,1)*t))
057
enddo
058
D(m)=-(0,1)*g(t)*(Y(m-1)*sqrt((m-1)/2.d0)*exp(-2.d0*(0,1)*t))
059
return
060
end
061
 
062
function g(t)
063
real(8) g,t
064
pi=3.1415926535d0
065
if(t&lt;=2000.d0)then
066
g=0.5*(1+sin(t/2000.*pi-pi/2.))
067
!    g=t/2000.d0
068
else
069
g=1.00d0
070
endif
071
end
072
 
073
real recursive function H(i,x) result(Hi)
074
select case(i)
075
case(1)
076
Hi=1.0
077
case(2)
078
Hi=2.0*x
079
case(3:)
080
Hi=2*x*h(i-1,x)-2*(i-1)*H(i-2,x)
081
end select
082
end
083
 
084
integer Recursive function JC(i) result(J)
085
if(i==0)then
086
J=1
087
else
088
J=i*JC(i-1)
089
endif
090
end
091
 
092
SUBROUTINE GRKT1(T,Y,M,H,N,Z,F,D,B)
093
DIMENSION Y(M),D(M),Z(M,N),A(4),B(M)
094
complex(8) Y,D,Z,A,B,H,X,TT
095
real(8) t
096
A(1)=H/2.0d0
097
A(2)=A(1)
098
A(3)=H
099
A(4)=H
100
DO I=1,M
101
Z(I,1)=Y(I)
102
enddo
103
X=T
104
DO J=2,N
105
CALL F(T,Y,M,D)
106
DO I=1,M
107
B(I)=Y(I)
108
enddo
109
DO K=1,3
110
DO I=1,M
111
Y(I)=Z(I,J-1)+A(K)*D(I)
112
B(I)=B(I)+A(K+1)*D(I)/3.0
113
enddo
114
TT=T+A(K)
115
CALL F(TT,Y,M,D)
116
enddo
117
DO I=1,M
118
Y(I)=B(I)+H*D(I)/6.0
119
enddo
120
DO I=1,M
121
Z(I,J)=Y(I)
122
enddo
123
T=T+H
124
enddo
125
T=X
126
RETURN
127
END

0

目前有 4 条留言    访客:3 条, 博主:1 条

  1. tusooa 2012年04月20日 下午10:19  @回复  Δ-49楼 回复
    Firefox Firefox GNU/Linux GNU/Linux

    蛋痛。

  2. lainme 2012年04月20日 上午10:21  @回复  Δ-48楼 回复
    Firefox Firefox GNU/Linux GNU/Linux

    传说中的fortran77风格?


    • 管理员
      hsyyf 2012年04月20日 下午12:07  @回复  ∇地下1层 回复
      Firefox Firefox Windows Windows

      主程序是我的free风格,子程序是现成的77风格。。。

  3. YeLee 2012年04月19日 下午10:14  @回复  Δ-47楼 回复
    Firefox Firefox GNU/Linux GNU/Linux

    说了蛋疼你信了吧。

给我留言

留言无头像?


×