burgers_examples.m

function res = burgers_examples(step) 
%function res = burgers_examples(step) 
%
% Burgers equation with 1D finite difference discretization
% demonstrating rarefaction waves and shocks moving
%
% The steps demonstrate the movement of a shock, the development of 
% a rarefaction wave for different initial values and velocity directions, 
% different ODE solvers, different inputs and initial conditions. 
% 
% As example solver I provide an implicit Euler solver and a visualization 
% routine with the model. In step=6 I demonstrate that also matlab ODE 
% solvers can be used (but are slower) than the implicit Euler scheme. 
%
% The model is "tough" from model reduction perspective, as it has a 
% moving discontinuity if a shock develops, which cannot very well be 
% modelled by low-dimensional spaces. 
%
% You can define the input signal u arbitrarily, also time-varying. In 
% the examples I set it to constant-in-time values. 
% However the implicit Euler solver has a time-step width set to be stable 
% for u in [0,1]^2. Similar the inital value can be chosen arbitrarily in 
% x0 \in [0,1]^(n_xi+1) guaranteeing a stable implicit Euler discretization.
% So you should prevent input or initial value going beyond this range.

% B. Haasdonk 8.6.2016

if nargin < 1
  step = 1;
end;

%%%% This parameter scales the spatial and time-domain resolution:
n_xi = 200; % choose between 10 and 1000 points

%%%% do not change the following values, as these nicely result in
%%%% stable numerical scheme and give shock at boundary at end time:
nt = 2* n_xi;
T = 0.25;

switch step
 case 1 % rarefaction wave
  
  % model mainly makes sense with constant u corresponding to
  % initial values
  
  disp('example of shock moving to the right');
  
  x_left = 1;
  x_right = 0;  
  v = 1; 
  model = burgers_1d_model(n_xi, T, x_left, x_right, v);
  u = @(t) [x_left, x_right]';
  x0 = model.x0;
  %  nt = model.n_xi;
  %  nt = 200;
  X = model.explicit_euler_time_integration(model,x0,u,nt);
  model.plot_results(model,X);

 case 2 % shock moving
  
  disp('example of rarefaction wave to the right');  
  
  %%%% change the following values between [-1,1]:
  x_left = 0;
  x_right = 1;  
  v = 1; 
  model = burgers_1d_model(n_xi, T, x_left, x_right, v);
  u = @(t) [x_left, x_right]';
  x0 = model.x0;
  X = model.explicit_euler_time_integration(model,x0,u,nt);
  model.plot_results(model,X);
  
 case 3 % same as 2 but reverted velocity

  disp('example of shock wave to the left');  

  x_left = 0;
  x_right = 1;  
  v = -1; 
  model = burgers_1d_model(n_xi, T, x_left, x_right, v);
  u = @(t) [x_left, x_right]';
  x0 = model.x0;
  X = model.explicit_euler_time_integration(model,x0,u,nt);
  model.plot_results(model,X);
  
 case 4 

    disp('example of rarefaction wave to the left');  

  x_left = 1;
  x_right = 0;  
  v = -1; 
  model = burgers_1d_model(n_xi, T, x_left, x_right, v);
  u = @(t) [x_left, x_right]';
  x0 = model.x0;
  X = model.explicit_euler_time_integration(model,x0,u,nt);
  model.plot_results(model,X);
  
 case 5 % use arbitrary input u

  disp('arbitrary input (boundary values) 0.5:');
  disp('result is shock and rarefaction wave');
  x_left = 0;
  x_right = 1;  
  v = 1; 
  model = burgers_1d_model(n_xi, T, x_left, x_right, v);
  u = @(t) [0.5,0.5]';
  x0 = model.x0;
  tic
  X = model.explicit_euler_time_integration(model,x0,u,nt);
  t = toc;
  disp(['computation took ',num2str(t),' sec.']);
  model.plot_results(model,X);
  
  case 6 % use matlab integrator for case 5

   disp('same as step 5 (rarefaction wave + shock) but with ode45:')
   x_left = 0;
   x_right = 1;  
   v = 1; 
   model = burgers_1d_model(n_xi, T, x_left, x_right, v);
   u = @(t) [0.5,0.5]';
   x0 = model.x0;
   
   odefun = @(t,x) model.f(x,u(t)); 
   tspan = [0,T/2,T];
   tic
   [tout, Xtransp] = ode45(odefun, tspan, x0);
   t = toc;
   disp(['computation took ',num2str(t),' sec.']);
   model.plot_results(model,Xtransp');
   
 case 7 % use other initial value generating a shock
  
  disp('development of shock');
  x_left = NaN; % not used
  x_right = NaN; % not used if having external x0 
  v = 1; 
  model = burgers_1d_model(n_xi, T, x_left, x_right, v);
  u = @(t) [0.5,0.5]';
%  x0 = model.x0;
  x0 = zeros(n_xi+1,1);
  x0(20:40) = (0:20)/20;
  x0(41:60) = (19:-1:0)/20;
  tic
  X = model.explicit_euler_time_integration(model,x0,u,nt);
  t = toc;
  disp(['computation took ',num2str(t),' sec.']);
  model.plot_results(model,X);

    
 otherwise
  error('step unknown');
end;