% [u,x,t,kk,amp] = int_fact_KdV(t_0,t_f,dt,NPrint,N)
%
% this solves the periodic KdV equation 
%
%         u_t + 6 u u_x + u_xxx = 0
%
% with periodic boundary conditions using spectral derivatives 
% and an integrating factor for the time-stepping.
% t_0 is the initial time, t_f is the
% final time, dt is the time-step size, N is the number of mesh-points, 
% and M is the number of printouts.  It also returns the power spectrum of the
% solution, where kk is the collection of wave numbers and amp = log(f_k^2)
% is the magnitude of the kth Fourier mode.
%
% N SHOULD BE A POWER OF 2!  
% (It will work no matter what, but will go much faster if N is a power
% of 2.)
%
% The initial data is set within the code, search for "SET INITIAL
% DATA".  
%
% To turn off the nonlinear effects, look within the tstep subroutine

function [u,x,t,kk,amp] = int_fact_KdV(t_0,t_f,dt,NPrint,N)

% define the mesh in space
dx = 2*pi/N;
x = 0:dx:2*pi-dx;
x = x';
% the wave-numbers...
kk = (-N/2+1):1:(N/2-1);

% total number of steps to be taken:
Ntot = (t_f-t_0)/dt;
% test to see if Ntot is an integer, to see if your choice of dt
% is compatible with your choice of t_0 and t_f
if abs(Ntot - round(Ntot)) > 100000*eps
  disp(' ')
  disp('dt is not consistent with your choice of t_0 and t_f.  Choose another dt!')
  disp(' ')
  return
end
Ntot = round(Ntot);
% number of steps to be taken before a (screen) printout
NRun = Ntot/NPrint;
% test to see if NRun is an integer to see if your choice of dt is
% consistent with NPrint
if abs(NRun - round(NRun)) > 1000*eps
  disp(' ')
  disp('NPrint is not consistent with dt.  Choose another NPrint!')
  disp(' ')
  return
end
NRun = round(NRun);
% define the vector t that is to be printed out...
t = t_0:NRun*dt:t_f;

% SET INITIAL DATA
% sigma = 20;
% sigma = 100;
% u_now = sigma/2*sech(sqrt(sigma)/2*(x-pi/2)).^2;
sigma1 = 75;
sigma2 = 250;
u_now = sigma2/2*sech(sqrt(sigma2)/2*(x-pi/2)).^2;
% u_now = zeros(size(u_now));
% u_now = u_now + sigma1/2*sech(sqrt(sigma1)/2*(x-pi)).^2;

% We're computing solutions of
%
%    u_t = - u_xxx + Non   where Non = - 6 u u_x
%
% In Fourier space, this is
%
%   d/dt u_k = i*k^3 u_k + Non_k
%
% which has the solution
%
% u_k(t) = u_k(0) exp(ik^3 t) + int_0^t exp(ik^3 (t-s)) Non_k(s) ds
%
% Approximating Non_k(s) with Non_k(0),
%
% u_k(t) = u_k(0) exp(ik^3 t) + i/k^3 Non_k(0) (1-exp(ik^3 t))
%
% and so in Fourier space,
%
% u_k_new = u_k_now exp(ik^3 dt) + i/k^3 Non_k_now (1-exp(ik^3 dt))
%
% u_0_new = u_0_now.
%
% To get a method with higher accuracy, I would've chosen a better
% approximation of the integral than the left-hand rule.

% save first printout:
u(:,1) = u_now;
v = fft(u_now); 
amp(N/2:N-1,1) = log10( abs(v(1:N/2)) + 1.e-16);
amp(1:N/2-1,1) = log10( abs(v(N/2+2:N)) + 1.e-16);

for jj=1:NPrint

  % take NRun time-steps
  for j=1:NRun
      
   % compute solution with one step of size dt
   u_coarse = tstep(u_now,dt,N);
   % compute solution with two steps of size dt/2
   u_fine = tstep(u_now,dt/2,N);
   u_fine = tstep(u_fine,dt/2,N);
   % use Richardson extrapolation to combine u_coarse and
   % u_fine to get a solution with local truncation
   % error of O(dt^3)...
   u_now = 2*u_fine - u_coarse;

  end
  % save printout
  u(:,jj+1) = u_now;
  v = fft(u_now); 
  amp(N/2:N-1,jj+1) = log10( abs(v(1:N/2)) + 1.e-16);
  amp(1:N/2-1,jj+1) = log10( abs(v(N/2+2:N)) + 1.e-16);
end
x(N+1) = 2*pi;
u(N+1,:) = u(1,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u_now = tstep(u_now,dt,N);

% We're computing solutions of
%
%    u_t = - u_xxx + Non   where Non = - 6 u u_x
%
% In Fourier space, this is
%
%   d/dt u_k = i*k^3 u_k + Non_k
%
% which has the solution
%
% u_k(t) = u_k(0) exp(ik^3 t) + int_0^t exp(ik^3 (t-s)) Non_k(s) ds
%
% Approximating Non_k(s) with Non_k(0),
%
% u_k(t) = u_k(0) exp(ik^3 t) + i/k^3 Non_k(0) (1-exp(ik^3 t))
%
% and so in Fourier space,
%
% u_k_new = u_k_now exp(ik^3 dt) + i/k^3 Non_k_now (1-exp(ik^3 dt))
%
% u_0_new = u_0_now.

% Non = - 6 u u_x = - 3 (u^2)_x
% Non_k = - 3 i k (u^2)_k

    Non = -3*u_now.^2;
    % now differentiate Non.
    Non = fft(Non);
    Non(1) = 0;
    for k=1:(N/2-1)
      Non(1+k) = i*k*Non(1+k);
      Non(N-k+1) = -i*k*Non(N-k+1);
    end
% at this point, I have Non_k_now computed.    If I want to remove the 
% nonlinear effects, then I set Non to zero as below
%    Non = zeros(size(Non));


% start by computing Fourier transform of u_now...
    v = fft(u_now);
    for k=1:(N/2-1)
      v(1+k) = v(1+k)*exp(i*k^3*dt) + i/k^3*Non(1+k)*(1-exp(i*k^3*dt));
      v(N-k+1) = v(N-k+1)*exp(i*(-k)^3*dt) + i/(-k)^3*Non(N-k+1)*(1-exp(i*(-k)^3*dt));
    end
% at this point, v is v_new_k    

% take the inverse fourier transform.  
    u_now = ifft(v);
    % for good measure, impose reality.
    u_now = real(u_now);