I tried to simulate soliton propagation by solving the nonlinear Schrödinger equation using the split-step Fourier method. The following is an example of the Matlab code copied from a textbook.
alpha = 0;
beta_2 = 1;
gamma = 1;
A = N*sech(T); % initial pulse shape
L = max(size(A));
delta_omega = 1/L/delta_T*2*pi;
omega = (-L/2:1:L/2-1)*delta_omega;
solution = A; % analytical solution (same as original it is periodic)
step_num_hist = round(10.^linspace(1,3.5,20)); % total number of steps
avg_err = zeros(3,max(size(step_num_hist)));
end_pt = 0.5*pi; % the end of simulation, one soliton period
for step_num = step_num_hist
h = end_pt/step_num; % step size
A_t = A;
for n=1:step_num
A_f = fftshift(fft(A_t));
A_f = A_f.*exp(-alpha*(h/2)-1i*beta_2/2*omega.^2*(h/2));
A_t = ifft(ifftshift(A_f));
A_t = A_t.*exp(1i*gamma*(abs(A_t).^2*h));
A_f = fftshift(fft(A_t));
A_f = A_f.*exp(-alpha*(h/2)-1i*beta_2/2*omega.^2*(h/2));
A_t = ifft(ifftshift(A_f));
end
comp_result = abs(A_t);
avg_err(1,m) = sum(abs(comp_result.^2 - abs(solution).^2))/max(size(solution))/max(abs(solution).^2); % avg relative intensity error, found in a paper
m = m+1
end
The code seems to produce something that makes sense, so I don't think it is wrong. The problem is that when the step number increases (step size h
decreases) beyond a certain number (in my case around 1e-3), the error between the computed solution and the analytical solution doesn't decrease any further. This happens in both Matlab and Python. I also tried implementing a few higher-order methods that I found in the literature. They do converge faster, but the 'average relative intensity error' does not go below somewhere about 1e-3, while a paper I found which compares various algorithms can go down to like 1e-7 or even 1e-10 (total number of steps also of the order of 10e4).
So my question is why is that? Why does the error decrease with step size initially but then stops decreasing after the step size gets to a certain point? What can possibly go wrong?
1e-3
, your analytical solution may be off by1e-3
. Otherwise you should look for other sources of numerical error. There's spatial and timestepping error: do you have enough terms in your FFT? Is the problem numerically stable in each of its calculations? $\endgroup$ – Chris Rackauckas Aug 23 '17 at 14:48sech(x)
is a bit similar to a Gaussian pulse, and I included about 2^12 or 2^13 terms from -5 to +5. $\endgroup$ – Physicist Aug 23 '17 at 15:13