burgers_1d_model.m

Go to the documentation of this file.

function model = burgers_1d_model(n_xi,T,x_left, x_right, v)
% simple model for Burgers PDE
% $d/dt x  + d/d_xi (1/2 * v * x^2) = 0$
% on the unit square xi in [0,1]
% with initial value
% x(.,0) = x0(.) = piecewise constant x_left and x_right on left/right half
%               of domain
% and dirichlet boundary values x_left and x_right and 
% discretization via Lax-friedrichs (central differences with
% numerical diffusion)
% The v is the velocity. the parameter vector is [x_left, x_right, v]
%
%
% The corresponding script showing several aspects: call 
%    burgers_examples(step), 
% where step can vary from 1 to 7. Please run these steps to get an 
% insight into the model.  
%
%The model provides the function f(.,.) for the ode system
%   d/dt x(t)  = f(x(t),u(t))  with x(0) = x0
%which represents a nonlinear non-viscous Burgers equation discretized with 
%upwind finite differences. This means, we have n_xi spatial intervals, 
%hence n_xi+1 points on the grid and a solution variable in each gridpoint, 
%hence the state x(t) is n_xi + 1 dimensional.
%The input u is assumed to be 2-dimensional: first entry left boundary 
%value x_left, second entry right boundary value 
%x_right, so that value will be "entering" the domain over time.
%x0 is arbitrary, but the model provides a "default" value in the sense of 
%having half the domain filled with both boundary values, a so called 
%"Riemann-problem" for the hyperbolic problem. Step 1-6 use the default 
%initial value, step 7 uses an arbitrary initial value.
%
%It is a model, that is scalable, i.e. one can arbitrarily 
%set the n_xi variable, that defines the number of spatial intervals 
%of the discretized Burgers equation. 

% B .Haasdonk 8.6.2016

model = [];
% number of x intervals
model.n_xi = n_xi; 
model.delta_xi = 1/n_xi;
model.grid = 0:1/n_xi:1;
model.x_left = x_left;
model.x_right = x_right;
model.v = v;
model.T = T;
model.lambda = 0.5; % for Lax-Friedrichs-flux: numerical diffusion
model.x0 = my_init_values(model);
model.f = @(x,u) burgers_nonlinearity(model,x,u);
model.explicit_euler_time_integration = @(model,x0,u,nt) ...
    explicit_euler_time_integration(model, x0, u, nt);
model.plot_results = @(model,X) plot_results(model,X);

function x0 = my_init_values(model);
% piecewise constant initial value for Riemann problem
x0 = zeros(model.n_xi+1,1);
mid_ind = floor((model.n_xi+1)/2);
x0(1:mid_ind) = model.x_left;
x0(mid_ind+1:end) = model.x_right;

function f = burgers_nonlinearity(model,x,u);
% For simplicity: Lax-Friedrichs flux
% g(x,w) = 1/2 (f(x) + f(w)) + 1/(2 lambda) (x-w)
% where (x,w) are neighbouring values of the state
% then
% 1/delta_t (x(t+delta_t)-x(t)) can be approximated by 
% -1/delta_x * (g(x,xshift_plus1 - g(x_shift_min1,x))

%% The following would be correct if boundary values instead of
%% input u were to be used:
%x_min1 = model.x_left; % left dirichlet value
%x_plus1 = model.x_right; % right dirichlet value
%xshift_plus1 = [x(2:end); x_plus1];
%xshift_min1 = [x_min1; x(1:end-1)];

v = model.v; % velocity
xshift_plus1 = [x(2:end); u(2)];
xshift_min1 = [u(1); x(1:end-1)];
%%%% Lax Friedrichs Flux
%f = + 1/ (2*model.delta_xi) * ...
%    ( v * ( xshift_plus1.^2 - xshift_min1.^2) ...
%      - 1/model.lambda * (2*x -xshift_plus1 - xshift_min1));    
%%%% Engquist Osher Flux (cmp: Kroener '97)
g_x_xshift_plus1 = fplus(model,x) + fminus(model,xshift_plus1);
g_xshift_min1_x = fplus(model,xshift_min1) + fminus(model,x);
f = - 1/ (model.delta_xi) * ...
    ( g_x_xshift_plus1 - g_xshift_min1_x );

function f = fplus(model,x) 
f = max(2*model.v*x,0) .* x;

function f = fminus(model,x) 
f = min(2*model.v*x,0) .* x;

function X = explicit_euler_time_integration(model, x0, u, nt)
% simple explicit Euler scheme
nx = length(x0);
delta_t = model.T/nt;
X = zeros(nx,nt+1);
x = x0(:);
X(:,1) = x;
for i=1:nt
  ut = u((i-1)*delta_t);
  xnew = x + delta_t* model.f(x,ut);
%  keyboard;
  X(:,i+1) = xnew;
  x = xnew;
end;

function p = plot_results(model,X)
ts = [0,floor((size(X,2)-1)/2), size(X,2)-1];
for i = 1:length(ts);
  subplot(length(ts),1,i)
  hold on;
  plot(model.grid,X(:,ts(i)+1)); 
  title(['state at t = ',num2str(ts(i)),' * \Delta t']); 
end;