2
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • 3
    $\begingroup$ If every method you're testing isn't converging to more than 1e-3, your analytical solution may be off by 1e-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:48
  • $\begingroup$ how do I know whether there are enough terms in the FFT? the initial pulse sech(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
  • $\begingroup$ another problem is that when I implement some of the higher-order methods, I sometimes get Nan from the FFT for small number of steps, but this problem disappears for step number larger than around 1000. I'm still figuring why, and I don't know if this has an effect on the convergence (the simple code I posted above doesn't seem to suffer from this problem) $\endgroup$ – Physicist Aug 23 '17 at 15:15
  • 1
    $\begingroup$ You know if you have enough terms by seeing how it converges. It the convergence when changing time is falling off, then you know you don't have enough spatial resolution. It's very problem and accuracy dependent how to balance these. But the second issue could be CFL stability problems. The spatial resolution and the timestepping method put constraints on the allowed size for the timesteps. $\endgroup$ – Chris Rackauckas Aug 23 '17 at 15:25
  • $\begingroup$ Can you write down the equations for your problem? $\endgroup$ – nicoguaro♦ Aug 23 '17 at 16:03

Your Answer

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

Browse other questions tagged or ask your own question.