% p31f.m - Solve KdV eq. u_t + uu_x + u_xxx = 0 on [-pi,pi] by
%         FFT with integrating factor v = exp(-ik^3t)*u-hat.

clf; clear;

function [x tdata udata] = KdVsolu(uex, L, N, tmax, dt, nplot) 
  x = (2*pi/N)*(-N/2:N/2-1)';
  u = feval(uex, x, 0);    v = fft(u); 
  k = [0:N/2-1 0 -N/2+1:-1]'/L; ik3 = 1i*k.^3;
  
  nplt = floor((tmax/nplot)/dt); nmax = round(tmax/dt);
  udata = u; tdata = 0; 
  for n = 1:nmax
    t = n*dt; g = -.5i*dt*k;
    E = exp(dt*ik3/2); E2 = E.^2;
    a = g.*fft(real( ifft(     v    ) ).^2);
    b = g.*fft(real( ifft(E.*(v+a/2)) ).^2);     % 4th-order
    c = g.*fft(real( ifft(E.*v + b/2) ).^2);     % Runge-Kutta
    d = g.*fft(real( ifft(E2.*v+E.*c) ).^2);
    v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
    if mod(n,nplt) == 0 
      u = real(ifft(v)); 
      udata = [udata u]; tdata = [tdata t];
    end
  end
end

A = .3; L=15; y0=-20/L;
u =inline('12*A^2*sech(A*L*(x-y0)-4*A^3*t).^2', 'x', 't'); 

tmax=60; dt=0.01; nplot=25;  Nvec=32:8:100;
err0 = []; err1 = []; 
for N=Nvec
    [x, tdata, udata] = KdVsolu(u, L, N, tmax, dt, nplot);  
    ue = feval(u, x, tdata(nplot+1));
    uerr=fft(udata(:,nplot+1)-ue,N)/N;
    err0 = [ err0 norm(abs(uerr), 2) ];                % L^2 norm
    k = [0:N/2 -N/2+1:-1]';
    err1 = [ err1 norm(abs(uerr).*sqrt(1+k.^2), 2) ];  % H^1 norm      
end
[tt, xx]=meshgrid(tdata,x);  plot3(tt,xx, udata);
xlabel t, ylabel x

figure;
semilogy(Nvec, err0, 'r+', 'markersize', 14), hold on
semilogy(Nvec, err1, 'go', 'markersize', 14),
grid on, xlabel N, ylabel error
title('The L^2(+) and H^1(o) convergence rates');



