% This code uses Matlab's syntax.
% This file is best viewed with Notepad or Wordpad,
% but not with Notepad++ (some of the lines will appear as shifted there).
% This is a logical flow-chart for implementing the averaging equation, Eq. (28),
% of the paper "Long-time simulations of nonlinear Schr¨odinger-type equations
% using step size exceeding threshold of numerical instability," by T.I. Lakoba,
% J. Sci. Comput. DOI 10.1007/s10915-016-0346-y
%
% Recall from the line below Eq. (25) that
% \hat{v}(n) = \hat{u}(n) * exp(i*k^2*dt*n). (1)
% Here \hat{u}(n) is the Fourier transform of the solution u at time step n,
% and dt is the time step.
% For the purpose of coding, it is convenient to write (28) as:
% vhat_new = c(1)*vhat(1) + c(2)*vhat(2) + ... c(7)*vhat(7), (2)
% where the averaging coefficients c(1) -- c(7) are defined below, and
% vhat(1) stands for \hat{v}(n-3), vhat(2) stands for \hat{v}(n-2), etc..
% Note: The variables \hat{v}(n) depend not only on time (=dt*n) but also on wavenumber, k.
% However, since this is text is just the logical flow-chart rather than an actual code,
% we will omit the k-dependence (and also the x-dependence later).
% Finally, note that equation (2) above can be rewritten in terms of variables \hat{u}:
% uhat_new = c(1)*uhat(1)*exp(-3*i*k^2*dt) + c(2)*uhat(2)*exp(-2*i*k^2*dt) + ... ,
% (3)
% where we have used equation (1) relating vhat and uhat.
% In the beginning of the code we specify how often we apply stabilization step:
tstabilize = 5; % time interval over which we stabilize the solution
nn_tstabilize = round(tstabilize/dt);
% the corresponding number of time steps
% Coefficients needed for averaging, as per Eq. (28):
c(1) = 1/64;
c(2) = -3/32;
c(3) = 15/64;
c(4) = 11/16;
c(5) = c(3);
c(6) = c(2);
c(7) = c(1);
% We also define exponential factors to relate \hat{v} with \hat{u}, as stated above:
ef(1) = exp(-3*i*k^2*dt);
ef(2) = exp(-2*i*k^2*dt);
ef(3) = exp(-1*i*k^2*dt);
ef(4) = 1;
ef(5) = exp(1*i*k^2*dt);
ef(6) = exp(2*i*k^2*dt);
ef(7) = exp(3*i*k^2*dt);
% As with vhat(n), these exponential factors depend on k. However, we do not indicate
% such a dependence explicitly in this logical flow-chart.
% -------------------------
% Calculations begin here:
for n = 1 : nmax
u = ... % this is solution u(n) computed by Leap-frog
u1 = ... % this is solution u(n+1) computed by Leap-frog
if n > 4 % we simply want the stabilization step being applied at very first steps
% We collect 8 = (7+1) values of u to use in the averaging formula (28).
% The reason why we need 1 more value than in (28) will become apparent soon:
% we need to obtain a new value not only for u(n) but also for u(n+1).
for j = 1 : 8
if mod(n + 3, nn_tstabilize) == j-1
ustab(j) = u; % this is the collected value of the solution;
% recall that u also depends on x, but we do not indicate
% such a dependence in this flow-chart
end
end
if mod(n, nn_tstabilize) == 4
% In order to avoid the accummulation of the parasitic error, we
% constract both u(n) and u(n+1) from scratch.
u = 0;
u1 = 0;
% "Build" the new solutions u(n) and u(n+1) by adding one term at a time
% as per Eq. (28).
for j = 1 : 7
u = u + c(j)*ifft( ef(j).* fft(ustab(j)) );
u1 = u1 + c(j)*ifft( ef(j).* fft(ustab(j+1)) );
end
% Note 1:
% For those unfamiliar with Matlab's syntax, fft and ifft stand for the
% Fast Fourier Transform and its inverse.
% Note 2:
% The occurrence of the exponential factors is clarified by equation (3) above.
end
% Strictly speaking, above we have reset the u,u1-values 3 time steps back and so
% need to also reset the time (i.e. the step counter n) 3 steps back.
% However, not doing so will introduce a negligible error (proportional to
% 3/nn_tstabilize << 1) in the total simulated time, and so we won't do it
% in order not to make the logic of the code more complex.
end % this ends the loop verifying if n > 4
end % this ends the loop of computational steps