%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% see how a cusp propagates on an elastic string and %
% how we can simulate that using finite difference.  %
% y_tt = y_xx - B * y_t; x periodic.                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

B = 0;
xmin = 0; xmax = 2*pi; mesh = 256; 
xdel = (xmax-xmin)/mesh; xdel2 = xdel*xdel;
tmin = 0; tmax = 1; 
r = 0.9; r2 = r*r;
tdel = r*xdel; tdel2 = tdel*tdel;
n = 0:mesh-1;
x = n*xdel+xmin;
yi = (x<pi).*x + (x>=pi).*(2*pi-x);
ydot = x*0; ylast = yi - ydot * tdel;
clf; plot (x,yi); hold on; 
axis([xmin xmax min(yi) max(yi)]);

% finite difference
tnow = tmin; y=yi;
while (tnow < tmax)
  yp = y;
  y = ( ([y(2:mesh) y(1)] + [y(mesh) y(1:mesh-1)] - 2*y)*r2 ...
      - ylast * ( 1 - B*tdel/2 ) + 2*y ) / ( 1 + B*tdel/2 );
  ylast = yp;
  tnow = tnow + tdel;
end;
plot (x, y, 'g--');

% spectral solution -- exact
yw = fft(yi); ydotw = fft(ydot);
k = 2*pi/(xmax-xmin)*...
    (n.*(n<mesh/2)-(mesh-n).*(n>=mesh/2));
channel = 10;
filter = abs(k)<channel*2*pi/(xmax-xmin);
% for each k: exp(ikx)*(a*cos(k*t)+b*sin(k*t))
a = filter.*yw; 
b = filter.*ydotw./k; b(1)=0;
% final profile
yw = a.*cos(k*tnow)+b.*sin(k*tnow);
yf = ifft(yw);
plot (x, yf, 'r-.');

title (['y_{tt} = y_{xx} - By_{t}, B = ' num2str(B) ...
    ', \Delta x = ' num2str(xdel) ', \Delta t = ' num2str(tdel) ...
    ', t_f = ', num2str(tnow)]);
xlabel ('x'); ylabel ('y');
text (pi+0.3,3.,'t=0');
text (pi,2,'t=t_f');
print -depsc2 cusp.eps;
