1
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$ – Kirill May 4 '19 at 22:37
  • $\begingroup$ 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. $\endgroup$ – decarat May 5 '19 at 2:18
  • $\begingroup$ 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. $\endgroup$ – decarat May 5 '19 at 2:48
  • $\begingroup$ @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 $\endgroup$ – Chip Jul 1 '19 at 6:02
  • $\begingroup$ Not an answer but there is absolutely no way you need 100,000 Fourier modes for initial data as smooth as a Gaussian. You should be able to get spatial resolution to close to machine epsilon for NLS with periodic BCs with a a few hundred to a thousand max $\endgroup$ – whpowell96 Jun 25 '20 at 0:09
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.