
>> N = 2^10;
>> dx = 2*pi/N;
>> x = -0:dx:2*pi-dx;
>> eps = 1;
>> u = sqrt(sin(x).^2 + eps);
>> plot(x,u), figure(1)
>> hold on
>> eps = .1;
>> u = sqrt(sin(x).^2 + eps);
>> plot(x,u), figure(1)
>> eps = .01;
>> u = sqrt(sin(x).^2 + eps);
>> plot(x,u), figure(1)
% as epsilon -> 0, the profile develops corners at x=0,pi.  Now look
% at how this affects the power spectrum

>> eps = 1;
>> u = sqrt(sin(x).^2 + eps);
>> [k,amp] = find_spec(u);
>> hold off
>> plot(k,amp); figure(1)
>> eps = .1;
>> u = sqrt(sin(x).^2 + eps);
>> [k,amp] = find_spec(u);
>> plot(k,amp); figure(1)
>> eps = .01;
>> u = sqrt(sin(x).^2 + eps);
>> [k,amp] = find_spec(u);
>> plot(k,amp); figure(1)

>> help heat4

  [u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N)
 
  this solves the periodic heat equation u_t = u_xx with initial data u_0 =
  cos(x) with periodic boundary conditions using spectral derivatives in
  space and explicit time-stepping.  t_0 is the initial time, t_f is the
  final time, N is the number of mesh-points, and M is the number of
  time steps.  err is the error. 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.)

% now we solve the heat equation with explicit time-stepping, having
% chosen a time step that won't to lead to trouble. 
>> [u,err,x,t,kk,amp] = heat4(0,1,100,8);
r =
    0.0162
r =
    0.9100

% plotting the power spectrum, we see that noise stay at the level of
% round-off
>> plot(kk,amp(:,1)); figure(1)
>> hold off;
>> plot(kk,amp(:,1)); figure(1)
>> hold on
>> plot(kk,amp(:,21)); figure(1)
>> plot(kk,amp(:,61)); figure(1)
>> plot(kk,amp(:,101)); figure(1)

% now we solve the heat equation with explicit time-stepping, having
% chosen a time step that is sure to lead to trouble.  (did this by
% taking dx much smaller)
>> [u,err,x,t,kk,amp] = heat4(0,1,100,128);
r =
    4.1501
r =
  -38.6900

>> clf; hold on;
>> plot(kk,amp(:,1)); figure(1)
>> plot(kk,amp(:,2)); figure(1)
>> plot(kk,amp(:,3)); figure(1)
>> plot(kk,amp(:,4)); figure(1)

% from the power spectrum, we can see the noise is growing.  now we look
% for this in the error.
>> plot(x,err(:,1))
>> figure(1)
>> plot(x,err(:,2)),figure(1)
>> plot(x,err(:,3)),figure(1)
>> plot(x,err(:,4)),figure(1)
>> plot(x,err(:,5)),figure(1)
>> plot(x,err(:,6)),figure(1)
>> plot(x,err(:,10)),figure(1)
>> plot(x,err(:,6)),figure(1)
>> plot(x,err(:,7)),figure(1)
% up until now, the error was around 1e-5.  Now we see that the error
% has grown hugely
>> plot(x,err(:,8)),figure(1)
% looking at the power spectrum, we see that the high modes have grown
% to be on the order of 1.e-5, which is why they're now visible in the
% error.  (They were always present, not just visible.)
>> plot(kk,amp(:,8)); figure(1)

% now solve with the same dt and dx, but with the implicit method:
>> [u,err,x,t,kk,amp] = heat5(0,1,100,128);
r =
    4.1501

% at this time the power spectrum is still fine.
>> plot(kk,amp(:,8)); figure(1)
% as is the error:
>> plot(x,err(:,8)),figure(1)

% now look at the power-spectrum for burger's equation.
>> load data.mat
>> who
Your variables are:

N         amp_512   k         kk_512    u_256     x_256     
amp       dx        kk        t         u_512     x_512     
amp_128   eps       kk_128    u         x         
amp_256   err       kk_256    u_128     x_128     

% plot the solution and watch it become singular in finite time:
>> plot(x_512,u_512(:,1)); figure(1); hold on;
>> plot(x_512,u_512(:,21)); figure(1);
>> plot(x_512,u_512(:,41)); figure(1);
>> plot(x_512,u_512(:,61)); figure(1);
>> plot(x_512,u_512(:,81)); figure(1);
>> plot(x_512,u_512(:,101)); figure(1);
>> hold off
% plot the power spectrum and watch it expand to fill all modes. 
>> plot(kk_512,amp_512(:,1)); figure(1);
>> hold on
>> plot(kk_512,amp_512(:,11)); figure(1);
>> plot(kk_512,amp_512(:,21)); figure(1);
>> plot(kk_512,amp_512(:,41)); figure(1);
>> plot(kk_512,amp_512(:,61)); figure(1);
>> plot(kk_512,amp_512(:,81)); figure(1);
>> diary off
