why the two split-step fourier methods using fft and ifft is equivalent ?

clear
close all
N = 512; % Number of Fourier modes
dt = .005; % Size of time step
tfinal = 2; % Final time
M = round(tfinal/dt); % Total number of time steps
J = 100; % Steps between output
L = 20; % Space period
q = -1;
mu = 2/3;
v = 1;
h = L/N; % Space step
n = (-N/2:1:N/2-1)'; % Indices
x = n*h; % Grid points
u0 = (12 - 48*x.^2 + 16*x.^4)./...
(8.*sqrt(6).*exp(x.*x/2.)*pi.^0.25);
u = u0; % Initial Condition Phi_4 of harmonic potential
k = 2*n.*pi/L; % Wavenumbers
t = 0;
for m = 1:1:M % Start time loop
u = exp(dt/2*1j*q*x.^2).*u; % Solve non-constant part of LSE
c = fftshift(fft(fftshift(u))); % Take Fourier transform
c = exp(dt/2*1j*-k.^2).*c; % Advance in Fourier Space
u = fftshift(ifft(fftshift(c))); % Return to Physical Space
end
plot(x,u0,'k',  'LineWidth', 2,'DisplayName', 'initial wavefunction')
hold on
plot(x,real(u),'r',  'LineWidth', 2,'DisplayName', 'final wavefunction real part')
plot(x,imag(u),'m',  'LineWidth', 2,'DisplayName', 'final wavefunction imag part')
legend('show');
u_an=exp(-1j*9/2*20).*u0;%*20 is not important, because it always in the eigen state and the modulus is constant
%just the real and imag of wavefunction is changed.
figure;
plot(x, abs(u), 'r.', 'MarkerSize', 12, 'DisplayName', '|u|');
hold on;
plot(x, abs(u_an), 'b-', 'LineWidth', 1.5, 'DisplayName', '|u_{an}|');
hold off;
legend('show');

second:@Paul told me that fft must need input from x=0, I can not figure out why the two code gives the same result.

clear
close all
N = 512; % Number of Fourier modes
dt = .005; % Size of time step
tfinal = 2; % Final time
M = round(tfinal/dt); % Total number of time steps
J = 100; % Steps between output
L = 20; % Space period
q = -1;
mu = 2/3;
v = 1;
h = L/N; % Space step
n = (-N/2:1:N/2-1)'; % Indices
x = n*h; % Grid points
u0 = (12 - 48*x.^2 + 16*x.^4)./...
(8.*sqrt(6).*exp(x.*x/2.)*pi.^0.25);
u = u0; % Initial Condition Phi_4 of harmonic potential
k = 2*n.*pi/L; % Wavenumbers
t = 0;
for m = 1:1:M % Start time loop
u = exp(dt/2*1j*q*x.^2).*u; % Solve non-constant part of LSE
c = fftshift(fft(u)); % Take Fourier transform
c = exp(dt/2*1j*-k.^2).*c; % Advance in Fourier Space
u = ifft(fftshift(c)); % Return to Physical Space
end
plot(x,u0,'k',  'LineWidth', 2,'DisplayName', 'initial wavefunction')
hold on
plot(x,real(u),'r',  'LineWidth', 2,'DisplayName', 'final wavefunction real part')
plot(x,imag(u),'m',  'LineWidth', 2,'DisplayName', 'final wavefunction imag part')
legend('show');

u_an=exp(-1j*9/2*20).*u0;%*20 is not important, because it always in the eigen state and the modulus is constant
%just the real and imag of wavefunction is changed.
figure;
plot(x, abs(u), 'r.', 'MarkerSize', 12, 'DisplayName', '|u|');
hold on;
plot(x, abs(u_an), 'b-', 'LineWidth', 1.5, 'DisplayName', '|u_{an}|');
hold off;
legend('show')

[Sulaymon Eshkabilov]

In code 1, you are shifting '0' freq component twice. Doing it once is sufficient. Therefore, using fftshit() once in code 2 is just good :).

[Daniel Niu]

Dear Sulaymon, why the code doesn't get the right result?

for m = 1:1:M % Start time loop
u = exp(dt/2*1j*q*x.^2).*u; % Solve non-constant part of LSE
c = (fft(fftshift(u))); % Take Fourier transform
c = exp(dt/2*1j*-k.^2).*c; % Advance in Fourier Space
u = (ifft(fftshift(c))); % Return to Physical Space
end

thank you!

Answer 1

[Paul]

Let's take a simplified view of the loops in code1 and code2, and renaming the variables appropriately. Also, because of the way that n is defined, it's correct to use ifftshift instead of fftshift in a few places (and also makes the code much more clear IMO) as shown below

% code1
u1 = exp(dt/2*1j*q*x.^2).*u1; % Solve non-constant part of LSE
%c1 = fftshift(fft(fftshift(u1))); % Take Fourier transform
c1 = fftshift(fft(ifftshift(u1))); % Take Fourier transform
c1 = exp(dt/2*1j*-k.^2).*c1; % Advance in Fourier Space
%u1 = fftshift(ifft(fftshift(c1))); % Return to Physical Space
u1 = fftshift(ifft(ifftshift(c1))); % Return to Physical Space


% code2
u2 = exp(dt/2*1j*q*x.^2).*u2; % Solve non-constant part of LSE
c2 = fftshift(fft(u2)); % Take Fourier transform
c2 = exp(dt/2*1j*-k.^2).*c2; % Advance in Fourier Space
%u2 = ifft(fftshift(c2)); % Return to Physical Space
u2 = ifft(ifftshift(c2)); % Return to Physical Space

If we ignore those exp() multiplications for a moment, the code would look like this:

% code1
u1 = u;
c1 = fftshift(fft(ifftshift(u1))); % Take Fourier transform
u1 = fftshift(ifft(ifftshift(c1))); % Return to Physical Space

% code2
u2 = u;
c2 = fftshift(fft(u2)); % Take Fourier transform
u2 = ifft(ifftshift(c2)); % Return to Physical Space

With these simplifications, it should be apparent that u1 == u == u2 by direct substitution.

For example with u2

u2 = ifft(ifftshift(c2));
u2 = ifft(ifftshift(fftshift(fft(u2))));
u2 = ifft(ifftshift(fftshift(fft(u))));
u2 = ifft(((fft(u))));   % ifftshift and fftshift cancel
u2 = u

and similarly for u1. Perhaps this gives an intuitive understanding of why code1 and code2 yield the same result.

But the code does have those exp() terms, particularly on c2, so it may be useful to have a more mathematical explanation.

In this problem, the input u is defined on an "fft-symmetric" set of discrete-time points, i..e, n = N/2:N/2-1 with N being even (if N were odd, we'd need n = -(N-1)/2:(N-1)/2).

As in code1 and code2, let U1 and U2 respectively be

U1 = fft(ifftshift(u1));
U2 = fft(u2);

U1 and U2 have the same magnitude but differ by a linear phase factor of the form exp(-1j*w*m) where w = (0:N-1)/N*2*pi and m = n(1) (example here), which is the DFT shift theorem. Then next step in the code is to multiply U1 and U2 by another exp() term, but the linear phase factor is still the difference between U1 and U2. If U1 and U2 differ by a linear phase factor, the DFT shift theorem (special case for fft-symmetric case) also tells us that

fftshift(ifft(U1)) == ifft(U2)

In short, because you're going from spatial to frequency to spatial, the shifting on input to fft and output from ifft in code1 don't really matter. But of course you have to do the shifting on c both ways in both code1 and code2 to correctly apply the exp() factor to c because k is defined as proportional to n.

Personally, I think code1 is much clearer to see what's happening, even if it's marignally less efficient. I'd also consider modifying code1 to look like

% code1
u1 = exp(dt/2*1j*q*x.^2).*u1; % Solve non-constant part of LSE
c1 = fftshift(fft(ifftshift(u1))); % Take Fourier transform symmetric in frequency
c1 = exp(dt/2*1j*-k.^2).*c1; % Advance in Fourier Space
c1 = ifftshift(c1); % convert back to DFT frequency
u1 = fftshift(ifft(c1)); % Return to Physical Space symmetric in time

[Nordine]

I like to apply the split step fourrier method to coupled non linear differential equations ( partial differential in time and space). that needs to use runge kutta ( or ode ) in loop .

is there someone code it?

[Paul]

You should consider opening a new question for more detail on what you want to do and what you've tried so far.