Fortran Coder

标题: 差分法解一维含时薛定谔方程，急... [打印本页]

作者: 薛定谔方程 时间: 2017-1-1 17:15
标题: 差分法解一维含时薛定谔方程，急...

简单解释：偏微分，长得像倒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

5秒内psy的运动.png (18.17 KB, 下载次数: 475)

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

初始值.png (6.05 KB, 下载次数: 480)

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

Screenshot_2017-01-01-21-46-56_com.tencent.mtt_14.jpg (116.93 KB, 下载次数: 491)

一维薛定谔方程

作者: 薛定谔方程 时间: 2017-1-1 21:57
这就是一维薛定谔方程，简单解释：偏微分，长得像倒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 (116.93 KB, 下载次数: 416)

作者: li913 时间: 2017-1-2 20:04
需要U和miu的表达式

作者: 薛定谔方程 时间: 2017-1-3 11:15

li913 发表于 2017-1-2 20:04
需要U和miu的表达式

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

作者: li913 时间: 2017-1-3 15:24
本帖最后由 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 (10.9 KB, 下载次数: 402)

作者: kerb 时间: 2017-1-3 20:11
时间步长选着不合适
你可以找本时域有限差分法（FDTD）看看，
附件是一个解释

作者: 薛定谔方程 时间: 2017-1-5 12:07

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 (6.35 KB, 下载次数: 478)

作者: 薛定谔方程 时间: 2017-1-5 12:10

kerb 发表于 2017-1-3 20:11
时间步长选着不合适
你可以找本时域有限差分法（FDTD）看看，
附件是一个解释

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

作者: 薛定谔方程 时间: 2017-1-5 12:15
会不会边界条件已经暗示了我的波函数只能是驻波？
不知道为何，就连我给的热传导程序，也都是个静止的图像。

[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

作者: kerb 时间: 2017-1-5 19:56
你搜一下courant condition
双曲类偏微分方程差分法求解，对dx,dy,dt等有一定要求，你自己搜一下

作者: kerb 时间: 2017-1-5 22:39

[Fortran] 纯文本查看 复制代码

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

这样求解不行
试一下这样

[Fortran] 纯文本查看 复制代码

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

作者: 薛定谔方程 时间: 2017-1-5 23:22

kerb 发表于 2017-1-5 22:39
do j=1,n-1
do i=2,m-1
x(i)=i*h

加个k的目的是什么？感觉除了一直算下去、然后输出同样的结果，并没有什么区别呀？另外，电脑已经算10分钟了

明天起床看看库朗条件吧，谢谢啦

作者: kerb 时间: 2017-1-6 11:48
你是用迭代法求解，你只迭代了一次，对于AX=B这样的方程至少要迭代A的阶数次才能得到近似解，你的那两层循环只是把边值在整个区域“匀和”了一次，你需要“匀和”N次结果才能比较接近解
或说回来，如果你先形成细数矩阵A，然后用直接法求解，X=A^{-1}B就不会有问题

作者: 薛定谔方程 时间: 2017-1-6 16:26
并不是太理解您的意思，但我还是让电脑算了一下，结果和以前一样（昨天可能运行其他内容，所以速度较慢吧，今天挺快的）。库朗条件貌似没有应用到复数范围内的吧？fdtd倒是有篇论文，或许可以借鉴一下，正在研究

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)