C语言方形障碍数值模拟中如何使高斯包移动

Question

我正在尝试使用高斯包通过 Trotter-Suzuki 公式和快速傅立叶变换 (FFT) 来研究遇到方形障碍时的传输概率，就像所做的那样 in this Quantum Python article。但我需要用C来实现它。原则上，波函数在与方形障碍物碰撞之前将保持其形状。但我发现波函数在与方形屏障碰撞之前随时间急剧变平。有人发现以下代码有问题吗？

此处，创建了两个文件 - result 和 psi.txt - 用于存储初始波函数和演化波函数。每个的前两个数据是 x 坐标，波函数在 x 处的概率。文件 result 中每行的第三个数据是方形障碍分布。我使用的 FFT 显示为 in this C program.

#include <stdio.h>
#include <math.h>

#define h_bar 1.0
#define pi 3.1415926535897932385E0
#define m0 1.0

typedef double real;
typedef struct { real Re; real Im; }  complex;

extern void fft(complex x[], int N, int flag);

complex complex_product(complex x, real y_power, real y_scale)
{//x*exp(i*y_power)*y_scale
    real Re, Im;

    Re = (x.Re*cos(y_power)-x.Im*sin(y_power))*y_scale;
    Im = (x.Re*sin(y_power)+x.Im*cos(y_power))*y_scale;

    x.Re = Re; x.Im = Im;
    return x;
}

real potential(real x, real a)
{
    return (x<0 || x>=a) ? 0 : 1;
}


void main()
{
    int t_steps=20, i, N=pow(2,10), m, n;
    complex psi[N];
    real x0=-2, p0=1, k0=p0/h_bar, x[N], k[N], V[N];
    real sigma=0.5, a=0.1, x_lower=-5, x_upper=5;
    real dt=1, dx=(x_upper-x_lower)/N, dk=2*pi/(dx*N);
    FILE  *file;

    file = fopen("result", "w");
    //initialize
    for (n=0; n<N; n++)
    {
        x[n] = x_lower+n*dx;
        k[n] = k0+(n-N*0.5)*dk;
        V[n] = potential(x[n], a);
        psi[n].Re = exp(-pow((x[n]-x0)/sigma, 2)/2)*cos(p0*(x[n]-x0)/h_bar);
        psi[n].Im = exp(-pow((x[n]-x0)/sigma, 2)/2)*sin(p0*(x[n]-x0)/h_bar);
    }

    for (m=0; m<N; m++)
        fprintf(file, "%g %g %g\n", x[m],     psi[m].Re*psi[m].Re+psi[m].Im*psi[m].Im, V[m]);
    fclose(file);

    for (i=0; i<t_steps; i++)
    {

        printf("t_steps=%d\n", i);
        for (n=0; n<N; n++)
        {
            psi[n]=complex_product(psi[n], -V[n]*dt/h_bar, 1);
            psi[n]=complex_product(psi[n], -k[0]*x[n], dx/sqrt(2*pi));//x--->x_mod
        }


        fft(psi, N, 1);//psi: x_mod--->k_mod

        for (m=0; m<N; m++)
        {
            psi[m]=complex_product(psi[m], -m*dk*x[0], 1);//k_mod--->k
            psi[m]=complex_product(psi[m], -h_bar*k[m]*k[m]*dt/(2*m0), 1./N);
            psi[m]=complex_product(psi[m], m*dk*x[0], 1);//k--->k_mod
        }

        fft(psi, N, -1);
        for (n=0; n<N; n++)
            psi[n] = complex_product(psi[n], k[0]*x[n], sqrt(2*pi)/dx);//x_mod--->x
    }

    file = fopen("psi.txt", "w");
    for (m=0; m<N; m++)
        fprintf(file, "%g %g 0\n", x[m], pow((psi[m]).Re, 2)+pow((psi[m]).Im, 2));
    fclose(file);   

}

我使用以下 Python 代码绘制初始和最终演化波函数：

call: `>>> python plot.py result psi.txt`

import matplotlib.pyplot as plt
from sys import argv

for filename in argv[1:]:
    print filename
    f = open(filename, 'r')
    lines = [line.strip(" \n").split(" ") for line in f]
    x = [float(line[0]) for line in lines]
    y = [float(line[2]) for line in lines]
    psi = [float(line[1]) for line in lines]
    print "x=%g, max=%g" % (x[psi.index(max(psi))], max(psi))

    plt.plot(x, y, x, psi)
#plt.xlim([-1.0e-10, 1.0e-10])
plt.ylim([0, 3])
plt.show()

Answer 1

此处未初始化变量 i:

 k[n] = k0+(i-N*0.5)*dk;

Answer 2

你的代码几乎是正确的，除了你在实域中缺少 initial/final 半步和一些不必要的操作 (k_mod -> k and back), 但主要问题是你的初始条件真的选错了。高斯波包的时间演化导致不确定性随时间呈二次方分布：

根据您选择的粒子质量和初始波包宽度，大括号中的项等于 1 + 4 t²。在一个时间步之后，波包已经比最初宽得多，而在另一个时间步之后变得比整个模拟盒更宽。使用 FFT 所隐含的周期性会导致空间和频率混叠，再加上过大的时间步长，这就是最终波函数看起来那么奇怪的原因。

我建议您尝试完全复制 Python 程序的条件，包括整个系统处于深势阱中的事实 (V_border -> +oo).

C语言方形障碍数值模拟中如何使高斯包移动

How to make gaussian package move in numerical simulation of a square barrier in C

c

physics

numerical-methods