我想使用分裂算子方法对谐振子中的位移高斯进行数值量子动力学（有关算法，请参见Hal Evans 的这些注释的底部）。

我对步骤 (3) 有疑问。我不确定如何表示动量项。它只是循环的索引吗？$p^2$

创建初始状态：

  unsigned int n=1400; ///> rank of the state vector psi(x). (x_0, x_1, ..., x_1399)^T
  const double dx=0.01;///> step of the state vector psi(x). dx = x_1-x_0;

  array1<Complex> psix(n,align); ///> array1 is the array class from fftw package that is wrapped into fftwpp. psi(x) is the state vector in the position space

  double x0 = 4; ///> the mean (center) of the gaussian of the initial state.

  for(unsigned int i=0; i < n; i++) { ///> go over all points in the position space
    double x = (i-0.5*n)*dx; ///> for each point, calculate x
    x+=x0; ///> shift the gaussian back by x0
    psix[i]=exp(-pow(x,2)); ///> create the gaussian for the initial state
  }

这是执行上述算法的代码部分。

  /// step 1
  for (size_t i = 0; i < n; ++i) {
    double x = (i-0.5*n)*dx; ///> get value of x for each element of position space vector
    ax[i]=exp(-ii*u(x)*dt*0.5*invhb)*psix[i];///> propagate half of the step using potential energy operator U:  a(x)=e^(-i*u*dt/(2*hb))*s(x,t)
  }
  /// step 2
  Forward.fft(ax,fs);///> do a fft into momentum space.
  /// step 3
  for (size_t i = 0; i < n; ++i) {
    double p=i;///> THE PROBLEM IS HERE!!! , assign the value of the momentum to be used in kinetic energy operator. I just used the omega from http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Definition
    bp[i]=exp(-ii*p*p*dt*0.5*invm)*fs[i];///> B(p) ≡ e^( −i p2 Δt / 2m) A(p) ///> act with a kinetic energy operator onto the momentum-space state vector.
  }
  /// step 4
  Backward.fftNormalized(bp,gfs);///> inverse fourier transform back into the position space
  /// step 5
  for (size_t i = 0; i < n; ++i) {
    double x = (i-0.5*n)*dx;///> convert indices of the position-space state vector into values of position
    psix[i]=exp(-ii*u(x)*dt*0.5*invhb)*gfs[i];///> propagate with potential energy operator through the other half-step: psi(x, t+Δt)  =  e^(−i U Δt / 2ℏ) B(x)
  }

这是我完整的 C++ 代码（基于 fftwpp 包装器示例 到 fftw 的教程）：

#include <limits>
#include <vector>
#include <complex>
//#include "gnuplot_i.hpp"
#include "Array.h"
#include "fftw++.h"

// instal fftw and 
// Compile with:
// g++ -std=c++11 -I fftw++-1.13 -fopenmp main.cc fftw++-1.13/fftw++.cc -lfftw3 -lfftw3_omp; ./a.out

using namespace std;
using namespace Array;
using namespace fftwpp;

double u(double x) {
  return x*x;
}
int main()
{
  fftw::maxthreads=get_max_threads();
  const Complex ii(0.0,1.0);
  const double invhb = 1.0;
  const double invm = 1.0;

  unsigned int n=1400;
  const double dx=0.01;
  const double dt=0.000001;
  unsigned int np=n/2+1;
  size_t align=sizeof(Complex);

  array1<Complex> psix(n,align);
  array1<Complex> fs(n,align);
  array1<Complex> gfs(n,align);

  fft1d Forward(n,-1,psix,fs);
  fft1d Backward(n,1,fs,gfs);

  double x0 = 4;

  for(unsigned int i=0; i < n; i++) {
    double x = (i-0.5*n)*dx;
    x+=x0;
    psix[i]=exp(-pow(x,2));
  }

  ofstream ofs("psix.dat");
  for (size_t i = 0; i < n; ++i) {
    double x = (i-0.5*n)*dx;
    ofs << x << "  " << real(psix[i]) << "  " << imag(psix[i]) << "\n";
  }

  int tn=10;
  //// following http://hep.physics.indiana.edu/~hgevans/p411-p610/material/08_pde_iv/schrodinger.html
  array1<Complex> ax(n,align); ///> a(x)=e^(-i*u*dt/(2*hb))*s(x,t)
  array1<Complex> bp(n,align); ///> B(p) ≡ e^( −i p2 Δt / 2m) A(p)
  for (size_t l = 0; l < 100; ++l) {
    for (size_t j = 0; j < tn; ++j) {
      /// step 1
      for (size_t i = 0; i < n; ++i) {
        double x = (i-0.5*n)*dx;
        ax[i]=exp(-ii*u(x)*dt*0.5*invhb)*psix[i];///> a(x)=e^(-i*u*dt/(2*hb))*s(x,t)
      }
      /// step 2
      Forward.fft(ax,fs);
      /// step 3
      for (size_t i = 0; i < n; ++i) {
        double p=i;
        //double p=i*dx;
        bp[i]=exp(-ii*p*p*dt*0.5*invm)*fs[i];///> B(p) ≡ e^( −i p2 Δt / 2m) A(p)
      }
      /// step 4
      Backward.fftNormalized(bp,gfs);
      /// step 5
      for (size_t i = 0; i < n; ++i) {
        double x = (i-0.5*n)*dx;
        psix[i]=exp(-ii*u(x)*dt*0.5*invhb)*gfs[i];///>psi(x, t+Δt)  =  e^(−i U Δt / 2ℏ) B(x)
      }
    }
    //// to plot gnuplots on the fly using gnuplot++, uncomment the lines below.
    //vector<double> xv,psir;
    //for (size_t i = 0; i < n; ++i) {
    //  xv.push_back((i-0.5*n)*dx);
    //  psir.push_back(real(psix[i]));
    //}
    //Gnuplot g("lines");
    //g.plot_xy(xv,psir);
    //std::cin.ignore( std::numeric_limits <std::streamsize> ::max(), '\n' );
  }

  ofstream offs("fs.dat");
  for (size_t i = 0; i < n; ++i) {
    double p=i*dx;
    offs << p << "  " << real(fs[i]) << "  " << imag(fs[i]) << "\n";
  }

  ofstream ofgfs("psixt.dat");
  for (size_t i = 0; i < n; ++i) {
    double x = (i-0.5*n)*dx;
    ofgfs << x << "  " << real(psix[i]) << "  " << imag(psix[i]) << "\n";
  }
}