# Split-step Fourier method applied on Schrodinger equation

I'm trying to solve a Schrodinger equation of the form $$i\frac{\partial}{\partial t}\psi=-\frac{\partial^2}{\partial x^2}\psi + (V(x)+\alpha|\psi|^2)\psi$$ using the split-step Fourier method implemented via MatLab code. In order to make sure it works, I'm testing my code with $$V(x)=x^2,\alpha=0$$ and $$\psi=Ce^{-x^2/2}$$ where $$C$$ is the normalization coefficient (it's set to 1 in this code because it shouldn't matter for testing purposes). My MatLab code is as follows:

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))))


The final and initial plots should be the same since $$\psi$$ is an eigenstate of the potential (harmonic oscillator). My question is: $$\text{What is the correct value for k?}$$ I'm fairly certain what I currently have as my k value is incorrect, despite it giving the correct answer within an error on the magnitude of 1e-13. I've tried a few, and a few give the correct answer within reasonable error, so I'm not sure which is the actual correct value.

• Isn't exp(-dt/2*1i*k.^2) just a constant in your code? Have you tried integrating $i\psi_t = -\psi_{xx}$ for a few steps and checking whether it's correctly integrating that component first, before integrating the potential term also? – Kirill May 4 at 22:37
• The issue isn't if it's integrating correctly, it's what the correct value(s) of k is so that it does integrate correctly. – decarat May 5 at 2:18
• I'm now using k = (-pi:2*pi/N:pi-2*pi/N); for k, and it seems to be working. Still not 100% certain if this is correct tho. – decarat May 5 at 2:48
• @decarat The above seems to be just the same as you had initially in your code. Looks valid. You have to shift your frequencies by $-\pi$, in order to have BOTH negative and positive ones. Take one array random and do a forward, backward FFT to see if they are identical. Sometime, one needs to divide by $N$. Also, to my knowledge, standard split-step does not apply to nonlinear PDEs, you may want to read more on that, eg hplgit.github.io/fdm-book/doc/pub/book/sphinx/._book018.html – Chip Jul 1 at 6:02

You can find a detailed description and the values for the wavenumbers in Trefethen's Spectral methods in MATLAB. Notice that if you generate the wavenumbers from $$0$$ to $$N$$, then you should apply fftshift/ifftshift functions.