Fortran Coder

查看: 12918|回复: 13

[其他行业算法] 差分法解一维含时薛定谔方程,急...

[复制链接]

8

帖子

1

主题

0

精华

入门

F 币
50 元
贡献
24 点
发表于 2017-1-1 17:15:04 | 显示全部楼层 |阅读模式
简单解释:偏微分,长得像倒e的那个符号,看作导数符号。基本就是差分法直接代替两个偏微分(一个一阶偏微分,一个二阶偏微分),化简一下,即得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))

先用Fortran写的程序,之后用origin画的图(图见附件),可是psy波函数解出来后应该是随时间而变化的,但我的基本没变化。(非物理方向专业可略去不看——这个方程类似一维热传导方程,只是相对多了一个虚数单位和一个psy函数,不知道会不会出现问题,方法二(未在图中展示)是我把psy(念为pusai)函数分解为实部和虚部,代回到薛定谔方程所得到的两个对应项的方程,然后仍旧是使用差分解法。)理论上两种方法基本没有区别,而实际上我的图也确实无异,但这更加使我疑惑,不知道问题出现在了哪里。
老师说我图像不对的原因是——M(空间的步数),N(时间的步数),Tao(即dx,空间步长)和h(即dt,时间步长)的值设置的不合理,而我试了一个晚上,终究还是没出来结果。其他考试临近,实在没太多时间执着于这个问题了,万望高人指点
[Fortran] 纯文本查看 复制代码
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

基本上就是初始值随t轴垂直于x轴平移,而理想中的应该不是垂直平移,而应该有一定的夹角 ... ...

基本上就是初始值随t轴垂直于x轴平移,而理想中的应该不是垂直平移,而应该有一定的夹角 ... ...

正态分布式的初始psy(归一化之前的波函数)-x图像

正态分布式的初始psy(归一化之前的波函数)-x图像

一维薛定谔方程

一维薛定谔方程

8

帖子

1

主题

0

精华

入门

F 币
50 元
贡献
24 点
 楼主| 发表于 2017-1-1 21:57:26 | 显示全部楼层
这就是一维薛定谔方程,简单解释:偏微分,长得像倒e的那个符号,看作导数符号。基本就是差分法直接代替两个偏微分(一个一阶偏微分,一个二阶偏微分),化简一下,即得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))
Screenshot_2017-01-01-21-46-56_com.tencent.mtt_14.jpg

790

帖子

2

主题

0

精华

大宗师

F 币
3765 元
贡献
2255 点
发表于 2017-1-2 20:04:43 | 显示全部楼层
需要U和miu的表达式

8

帖子

1

主题

0

精华

入门

F 币
50 元
贡献
24 点
 楼主| 发表于 2017-1-3 11:15:16 | 显示全部楼层
li913 发表于 2017-1-2 20:04
需要U和miu的表达式

我的u是以V的形式给的,miu是mu

790

帖子

2

主题

0

精华

大宗师

F 币
3765 元
贡献
2255 点
发表于 2017-1-3 15:24:53 | 显示全部楼层
本帖最后由 li913 于 2017-1-3 16:13 编辑

1、我改了初始条件;2、tao参数不合适。把t1改成20后,第30个点随时间变化如图;
3、你的差分公式貌似不对。
[Fortran] 纯文本查看 复制代码
program time_dependent_schrodinger_equation
    implicit none
    integer,parameter::M=150,N=500
    real(8)::L=150d0, T1=20, hba=1d0,mu=0.5d0
    integer::i,j
    real,external::V !势能函数
    real(8) ::psy2(m,n)
    complex(8):: cj=cmplx(0.0,1.0), psy(M,N)
    real(8)::h, tao, pi=3.1415926d0
    h=L/M
    tao=T1/N
!给出初始边界及初始值
    psy = 0
    psy2 = 0
    do i=2,M-1
        psy(i,1)=sin(2*pi/m*i)   !!阱内给出初始值
    end do                       

!中心差法解薛定谔方程
    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))    !!核心内容,将中心差法代入薛定谔方程并整理
            psy(i,j+1) = v(i)*psy(i,j) - (psy(i+1,j)-2*psy(i,j)+psy(i-1,j))/2.0/mu + cj*h/tao*psy(i,j)
            psy(i,j+1) = psy(i,j+1) *tao/h/cj
            psy2(i,j+1)=(abs(psy(i,j+1)))**2
            
        end do
    end do
    do i=1,n
      write(10,*) i,psy2(30,i)
    end do
pause
    end
!给出势函数
real    function V(i)
        implicit none
        integer::i
        v=0
        if (abs(i-75)-5<1.e-6) then
          V=10.
        end if
end
QQ截图20170103160402.jpg

评分

参与人数 1贡献 +9 收起 理由
vvt + 9

查看全部评分

59

帖子

2

主题

0

精华

大师

F 币
810 元
贡献
476 点
发表于 2017-1-3 20:11:06 | 显示全部楼层
时间步长选着不合适
你可以找本时域有限差分法(FDTD)看看,
附件是一个解释 未标题-1.jpg


评分

参与人数 1贡献 +9 收起 理由
vvt + 9 赞一个!

查看全部评分

8

帖子

1

主题

0

精华

入门

F 币
50 元
贡献
24 点
 楼主| 发表于 2017-1-5 12:07:16 | 显示全部楼层
li913 发表于 2017-1-3 15:24
1、我改了初始条件;2、tao参数不合适。把t1改成20后,第30个点随时间变化如图;
3、你的差分公式貌似不对 ...

麻烦你再推导一下你的公式,我推导过程如图。另外,我想要的是三维图像(x-t-psy的平方),当将你的数据输出修改后,图像差别太大。你为何把归一化条件给删掉了呢?
Screenshot_2017-01-05-12-02-43_com.miui.gallery.p.JPG

8

帖子

1

主题

0

精华

入门

F 币
50 元
贡献
24 点
 楼主| 发表于 2017-1-5 12:10:37 | 显示全部楼层
kerb 发表于 2017-1-3 20:11
时间步长选着不合适
你可以找本时域有限差分法(FDTD)看看,
附件是一个解释

我根据c(光速)=dx/dt作了图。M=15000,N=50,l=1500,T1=0.05,可是图像也并未移动。另外,在我的程序里,我应该是给光速暗含了重定义,并不等与传统的10**8量级,所以应该并不作影响吧?

8

帖子

1

主题

0

精华

入门

F 币
50 元
贡献
24 点
 楼主| 发表于 2017-1-5 12:15:38 | 显示全部楼层
会不会边界条件已经暗示了我的波函数只能是驻波?
不知道为何,就连我给的热传导程序,也都是个静止的图像。
[Fortran] 纯文本查看 复制代码
program main
    implicit none
    integer,parameter::M=150,N=500
    real,dimension(M)::x,x1
    real,dimension(N)::t
    real,dimension(M,N)::y
    integer::i,j
    real::l=150,t1=5,pi,h,tao
    pi=asin(1.)*2
    h=l/m
    tao=t1/n
    do j=1,N
        y(1,j)=0
        y(m,j)=0
    ENd do
    do i=2,M-1
        x1(i)=i*h
        y(i,1)=exp(((-1)*(x1(i)*0.2-5)**2./2.)/(2.*pi))
    end do
    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
    end program

59

帖子

2

主题

0

精华

大师

F 币
810 元
贡献
476 点
发表于 2017-1-5 19:56:19 | 显示全部楼层
你搜一下courant condition
双曲类偏微分方程差分法求解,对dx,dy,dt等有一定要求,你自己搜一下
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-3-28 23:56

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表