% turn on the diary
>> diary 8_17_02

% I know I have a file "burgers.m" but I can't remember much about it.
% so I type "help burgers".  This shows me whatever commented-out lines
% left to myself as instructions.
>> help burgers

  [u,x,t,kk,amp] = burgers(t_0,t_f,M,N)
 
  solves the periodic Burgers equation without viscosity:
                 u_t + u u_x = u_t + (u^2/2)_x = 0 
  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.)
 
  Note: for burgers equation with initial data cos(x), we know that the 
  solution blows up at time: t = 1

% I define the input parameters..
>> t_0 = 0; t_f = 1; M = 1000; N = 256;

% I now execute the program...
>> [u,x,t,kk,amp] = burgers(t_0,t_f,M,N);

% I want to see what objects I have around and what their dimensions
% are
>> whos
  Name      Size           Bytes  Class

  M         1x1                8  double array
  N         1x1                8  double array
  amp     128x1001       1025024  double array
  kk        1x128           1024  double array
  t         1x1001          8008  double array
  t_0       1x1                8  double array
  t_f       1x1                8  double array
  u       256x1001       2050048  double array
  x       256x1             2048  double array

Grand total is 385773 elements using 3086184 bytes

% I see that u has dimensions of x by t.  I check this by
% freezing a t (the first time, index 1) and checking that
% u(:,1) has the same dimensions as x...
>> size(u(:,1))

ans =

   256     1

% I now want to watch a movie of my solutions...  I don't feel
% like watching all 1001 samples, so I only wach some of them.
% A slightly better-written code would print out only a subset
% of the solutions.
>> for i=1:41; plot(x,u(:,25*(i-1)+1)); hold on; plot(x+2*pi,u(:,25*(i-1)+1)); title('cos initial data, shocks at t = 1'); xlabel(t(25*(i-1)+1)); hold off; figure(1); axis([0,4*pi,-1.5,1.5]); pause(1); end

% I want to look at the power spectrum.  Since I have 256 
% mesh-points, I can resolve up to 128 modes.
>> for i=1:41; plot(kk,amp(:,25*(i-1)+1)); xlabel(t(25*(i-1)+1)); figure(1); axis([0,128,-20,3]); pause(1); end

% I see that the solution becomes numerically unresolved when the
% highest modes (kk=128) are no longer at the level of roundoff
% error.  Also, I see that as time passes, the high modes like
% kk = 128 have fourier coefficients that are on the order of the
% low modes like kk=1.  This is why you see the oscillations in 
% the late-time solutions...  For more about what this means, look
% at Rob's very nice notes which you can find linked off of my
% page (go one up) on computing the NLS

% turn off the diary
>> diary off