薛定谔方程的分步傅里叶法

计算科学 傅里叶变换 运算符拆分
2021-12-07 18:27:40

我正在尝试求解形式为薛定谔方程itψ=2x2ψ+(V(x)+α|ψ|2)ψ使用通过 MatLab 代码实现的分步傅里叶方法。为了确保它有效,我正在测试我的代码V(x)=x2,α=0ψ=Cex2/2在哪里C是归一化系数(在此代码中设置为 1,因为它对测试目的无关紧要)。我的 MatLab 代码如下:

N = 100000;                         % Number of Fourier mode
dt = .001;                          % Time step
tfinal = 5;                         % Final time
M = round(tfinal/dt);               % Total number of time steps
L = 5;                              % Total space length
h = L/N;                            % Space step
n =( -N/2:N/2-1);                   % Indices
x = 2*n*h;                          % Grid points
u_i = exp(-((x-0).^2)/2);           % Intial pulse
u = u_i;                            % Make a duplicate of the initial pulse in order to compare at the end.
k = 2*pi/N * n;                     % wavenumber values.
epsilon = 0;                        % nonlinear coefficient
optical_potential=.5*(x-0).^2;

figure(1);                          % Plotting the probability function of the initial wavefunction (not normalised)
plot(x,abs(u_i).^2);
hold on;
plot(x,optical_potential);
hold off;


for m = 1:1:M/2 % Start time loop (M/2 since each loop is 2 time steps)
    c = (fft(u));                   % Half time-step linear propagation
    c = exp(-dt/2*1i*k.^2).*c;
    u = ifft(c);                    % Full time-step nonlinear propagation
    u = exp(dt*1i*(optical_potential + epsilon *(abs(u).^2))).*u; 
    c = (fft(u));                   % Half time-step linear propagation
    c = exp(-dt/2*1i*k.^2).*c;
    u = ifft(c); 
end

u_out=u; 
figure(2);                          % Plotting the probability function of the final wavefunction (not normalised)
plot(x,abs(u_out).^2);
hold on;
plot(x,optical_potential);
hold off;

error_percent=sum(abs(((abs(u_i) - abs(u_out)))))/sum(abs((abs(u_i))))

最终图和初始图应该相同,因为ψ是势的本征态(谐波振荡器)。我的问题是: What is the correct value for k? 我相当确定我目前的 k 值不正确,尽管它在 1e-13 大小的误差范围内给出了正确答案。我已经尝试了一些,并且一些在合理的错误范围内给出了正确的答案,所以我不确定哪个是实际的正确值。

1个回答

您可以在 MATLAB 的 Trefethen 的 Spectral 方法中找到波数的详细描述和值请注意,如果您从0N,那么你应该应用fftshift/ifftshift函数。