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?