%
%  w4a.m -- spectral diffusion w/ nonlinearity
%           -- bifurcation & instability
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  parameters

global n2 N D a b

choice = 2;

switch choice
  case {1,2}
    D =  1.0;  a =  1.2;   b = 1.0;
    N =   64;  L = pi/2;  dx = 2*L/N;  
    T =   24; dt =  1.2;
end

n2 = (( (0:1:N-1)*(pi/2/L) ).^2)';

%  set grids & coefficients

xg = 0:dx:2*L-dx;  x = [xg 2*L];
tg = 0:dt:T;

%  IV initialization (column), FFT & plot

switch choice
  case 1
    IVs  = (4/pi^2) * (xg.*(2*L-xg))';
  case 2
    IVs  = ((xg.*(2*L-xg) .* (sin(3*xg)+cos(5*xg)))/4)';
end

u = IVs;
usave = [u ; u(1)];

u2 = [u ; 0 ; -flipud(u(2:end))];
uf = -imag(fftshift(fft(u2))/N);
c  = uf(N+1:2*N);

figure(1);  clf;  subplot(2,1,1) 
plot(x,usave(:,1),'g');  hold on;  plot(x,usave(:,1),'k.');
axis([0 2*L -1.5 1.5]);

%  timestep PDE using matlab ODE solver

Xode = 'X01';  method = 'ode15s';
options = odeset('reltol',1e-3,'abstol',1e-6,'refine','1','stats','off');

disp(' ')
disp([' spectral diffusion: ' Xode])
disp(' ')

jj = [];  tj = [];
for tt = 0:dt:T-dt

   IVs = c;
   [t y] = feval(method,Xode,[tt tt+dt],IVs,options);

%  reconstruct u 

   c  = y(end,:)';
   uf = [0 ; -flipud(c(2:end)) ; c];
   u2 = ifft(fftshift(-i*uf))*N;
   u  = real(u2(1:N));  u = [u ; u(1)];

%  variational energy

   ud = diff(u);
   u_x = [ud(1)/dx  ; (ud(1:end-1) + ud(2:end)) / (2*dx) ; ud(end)/dx];

   jj = [jj ; trapz((D*u_x.^2 - a*u.^2 + (b/2)*u.^4)/2)*dx];
   tj = [tj ; t(end)];

% interval plotting
   usave = [usave u];
   plot(x,u);   axis([0 2*L -2 2]);  pause(0.1)
end

plot(x,u,'r');  hold on;  plot(x,u,'k.');

xlabel('\bf x-axis');  ylabel('\bf u-axis')
title('\bf diffusion & nonlinearity')

subplot(2,1,2);  plot(tj,jj,'x')
xlabel('\bf t-axis');  ylabel('\bf J[u]-axis')
title('\bf energy dissipation')