| 并不是太理解您的意思,但我还是让电脑算了一下,结果和以前一样(昨天可能运行其他内容,所以速度较慢吧,今天挺快的)。库朗条件貌似没有应用到复数范围内的吧?fdtd倒是有篇论文,或许可以借鉴一下,正在研究 |
|
你是用迭代法求解,你只迭代了一次,对于AX=B这样的方程至少要迭代A的阶数次才能得到近似解,你的那两层循环只是把边值在整个区域“匀和”了一次,你需要“匀和”N次结果才能比较接近解 或说回来,如果你先形成细数矩阵A,然后用直接法求解,X=A^{-1}B就不会有问题 |
kerb 发表于 2017-1-5 22:39 加个k的目的是什么?感觉除了一直算下去、然后输出同样的结果,并没有什么区别呀?另外,电脑已经算10分钟了 ![]() 明天起床看看库朗条件吧,谢谢啦 |
|
[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode
program time_dependent_schrodinger_equation
implicit none
integer,parameter::M=150,N=500
integer::i,j
real*16::hba=1.,mu=1./2
!hba即(普朗克常数/2/pi),mu为质量
real*16::l=150,T1=5
real,external::V !势能函数
real*16,dimension(M)::x,x1
real*16,dimension(N)::t,normalize,psy4
real*16,dimension(M,N)::psy2,psy3 !psy2为对复数psy取模并求其平方的大小,psy3,即为归一化后的波恩诠释下的概率。
complex*16,dimension(M,N)::psy
complex::cj
real*16::h,tao,pi
h=l/M
tao=T1/N
pi=asin(1.)*2
cj=(0.,1.) !虚数单位
!给出初始边界及初始值
do j=1,N
psy(1,j)=0
psy(M,j)=0 !边界设定为无限深势阱
end do
do i=2,M-1
x1(i)=i*h
psy(i,1)=exp(((-1)*(x1(i)*0.2-5)**2./2.)/(2.*pi)) !阱内给出初始值
end do !x(i)比着正常多乘上个0。2,平移了波的位置
!中心差法解薛定谔方程
do j=1,N-1
do i=2,m-1
psy(i,j+1)=psy(i,j)+cj*tao*((psy(i+1,j)-2*psy(i,j)+psy(i-1,j))/(h**2)-V(i)*psy(i,j)) !核心内容,将中心差法代入薛定谔方程并整理
psy2(i,j+1)=(abs(psy(i,j+1)))**2
! write(13,*) psy2(i,j+1)
end do
end do
!归一化
do j=2,n
normalize(j)=(sum(psy2(1:m,j)))!**0.5
end do
do j=2,N
do i=2,m-1
x(i)=i*h
t(j)=j*tao
psy3(i,j)=psy2(i,j)/(normalize(j))!**2
write(13,*) x(i),t(j),psy3(i,j)
end do
psy4(j)=sum(psy3(1:M,j))
! write(14,*) psy4(j) !检验归一化的正确性
end do
stop
end
!给出势函数
real function V(i)
implicit none
integer::i
v=0
if (abs(i-75)-5<1.e-6) then
V=10.
end if
! if abs(x-40)<3&abs(x-60)<3
! V=10
return
end这样求解不行 试一下这样 [Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode do j=1,n-1
do i=2,m-1
x(i)=i*h
t(j+1)=(j+1)*tao
y(i,j+1)=y(i,j)+tao*((y(i+1,j)-2*y(i,j)+y(i-1,j))/(h**2))
write(22,*) x(i),t(j+1),y(i,j+1)
end do
end do |
|
你搜一下courant condition 双曲类偏微分方程差分法求解,对dx,dy,dt等有一定要求,你自己搜一下 |
|
会不会边界条件已经暗示了我的波函数只能是驻波? 不知道为何,就连我给的热传导程序,也都是个静止的图像。 [Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode do k=1,m*n
do j=1,n-1
do i=2,m-1
x(i)=i*h
t(j+1)=(j+1)*tao
y(i,j+1)=y(i,j)+tao*((y(i+1,j)-2*y(i,j)+y(i-1,j))/(h**2))
if(k==m*n)write(22,*) x(i),t(j+1),y(i,j+1)
end do
end do
enddo |
kerb 发表于 2017-1-3 20:11 我根据c(光速)=dx/dt作了图。M=15000,N=50,l=1500,T1=0.05,可是图像也并未移动。另外,在我的程序里,我应该是给光速暗含了重定义,并不等与传统的10**8量级,所以应该并不作影响吧? |
捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )
GMT+8, 2025-11-1 19:21