advection_Nadapt.m

Go to the documentation of this file.

% small script demonstrating RB-approach with adaptive N choice 
% for a simple advection problem. Only dirichlet-data
% parametrization, should be made more complex later.
% Demonstrates the unequal error estimator increase
% is a motivation for model adaptivity.
%
% The results of the MATHMOD-Paper on Nadaptivity are generated
% with this script.

% Bernard Haasdonk 28.10.2008

%step = 0; % check continuity of dirichlet data
%step = 1; % single simulation
%step = 2; % generate reduced basis, detailed_data is stored
%step = 3 % load detailed_data and run reduced simulation with error plot
step = 4 % run n-adaptive reduced simulation with errot plot
%step = 5 % generate figures of error estimator and N evolution for paper
%step = 6 % perform time measurement for paper and generate error
%         % estimator figure

% default: load 'better' basis of step 3:
% Basis with large error 1e-2:
detailedfname = 'advection_Nadapt_detailed.mat';
%detailedfname = 'advection_adaptivity_detailed_new.mat';

params = [];
% load demo_lin_evol_params;
% ===>
%                            xrange: [0 1.0000e-003]
%                            yrange: [0 2.0000e-004]
%                                 T: 0.5000
%                                nt: 200
%                  bnd_rect_corner1: [2x5 double]
%                  bnd_rect_corner2: [2x5 double]
%                    bnd_rect_index: [-1 -2 -2 -1 -2]
%             name_dirichlet_values: 'weighted_boxes'
%                    dir_box_xrange: {[-2.2204e-016 5.0000e-004]  [5.0000e-004 1.0000e-003]}
%                    dir_box_yrange: {[2.0000e-004 2.0000e-004]  [2.0000e-004 2.0000e-004]}
%                             c_dir: 1
%                              beta: 1
%                name_neuman_values: 'rightflow'
%                  name_diffusivity: 'homogeneous'
%           name_diffusive_num_flux: 'gradient'
%                                 k: 0
%                         name_flux: 'gdl2'
%           use_velocity_matrixfile: 1
%                     divclean_mode: 'none'
%                            lambda: 2.0000e-007
%               velocity_matrixfile: 'vel_gdl2_200x40_l2e-07.mat'
%          name_convective_num_flux: 'lax-friedrichs'
%                  name_init_values: 'wave'
%                            freq_x: 3.1416e+004
%                            freq_y: 0
%                            c_init: 1
%                           verbose: 9
%                     show_colorbar: 1
%                     colorbar_mode: 'EastOutside'
%                          mu_names: {'c_init'  'beta'  'k'}
%                         mu_ranges: {[0 1]  [0 1]  [0 5.0000e-008]}
%                L_I_inv_norm_bound: 1
%                    L_E_norm_bound: 1
%     detailed_simulation_algorithm: 'detailed_lin_evol_simulation'
%               operators_algorithm: 'fv_operators_implicit_explicit'
%             init_values_algorithm: 'fv_init_values'
%                        lxf_lambda: 1.0194e+003
%                data_const_in_time: 1
%                     xnumintervals: 200
%                     ynumintervals: 40
%                          gridtype: 'rectgrid'
%                        error_norm: 'l2'
%             coercivity_bound_name: 'fv_coercivity_bound'
%                      alpha_over_k: 1.9915e+007
%                 energy_norm_gamma: 0.5000
%                        axis_equal: 1
%                   rb_problem_type: 'lin_evol'
%    inner_product_matrix_algorithm: 'fv_inner_product_matrix'
%                   error_algorithm: 'fv_error'
%            RB_extension_algorithm: 'RB_extension_PCA_fixspace'
%                   RB_stop_timeout: 3600
%                   RB_stop_epsilon: 1.0000e-008
%                      RB_stop_Nmax: 150
%                RB_error_indicator: 'estimator'
%                RB_generation_mode: 'uniform_fixed'
%                   RB_numintervals: [4 4 4]
%                       flux_linear: 1


% the following is required for step 0 and 1:

params.T = 2;

params.gridtype = 'triagrid'; % required for gui!
params.grid_initfile = '2dsquaretriang';
eps = 1e-6;
params.bnd_rect_corner1= [-eps, -eps ;  -eps,-eps; 1-eps, -eps; -eps ,1-eps]';
params.bnd_rect_corner2= [1+eps, +eps; 1+eps,+eps; 1+eps,1+eps; 1+eps, 1+eps]';
%bnd_rect_corner1= [-1, -1; eps,eps];
%bnd_rect_corner2= [2, 2; 2,2]';
params.bnd_rect_index= [-1; -1; -2; -2];
params.nt = 500;
grid = triagrid(params);

%plot_element_data_sequence(grid,grid.NBI.*(grid.NBI<0),params);
%keyboard;

%params.grid_initfile = 'circle_grid';
%grid = triagrid(params);
%plot(grid);
%params.nt = 300;

params.name_dirichlet_values = 'decomp_function_ptr';
params.dirichlet_values_components_ptr = @advection_Nadapt.my_udir_components;
params.dirichlet_values_coefficients_ptr = @advection_Nadapt.my_udir_coefficients;
params.name_init_values = 'decomp_function_ptr';
params.init_values_components_ptr = @advection_Nadapt.my_u0_components;
params.init_values_coefficients_ptr = @advection_Nadapt.my_u0_coefficients;
params.c_dir_1 = 1;
params.c_dir_2 = 1;
%params.c_dir_3 = 1;
params.name_neuman_values = 'outflow';

params.transport_x = 1;
params.transport_y = 1;
params.name_flux = 'transport';
params.flux_linear = 1;
params.name_convective_num_flux = 'engquist-osher';

params.name_diffusivity = 'homogeneous';
params.name_diffusive_num_flux = 'gradient';
params.k = 0.00;

params.detailed_simulation_algorithm = 'detailed_lin_evol_simulation';
params.operators_algorithm = 'fv_operators_implicit_explicit';
params.init_values_algorithm = 'fv_init_values';
params.data_const_in_time = 0;
params.rb_problem_type= 'lin_evol';
params.verbose = 9;
%params.verbose = 10;

% the following is required for step 2 

%params.mu_names= {'c_dir_1'  'c_dir_2'  'c_dir_3'};
params.mu_names= {'c_dir_1'  'c_dir_2'  'k'};
params.mu_ranges= {[0 1]  [0 1]  [0 1]};
params.inner_product_matrix_algorithm= 'fv_inner_product_matrix';
params.L_I_inv_norm_bound = 1;
params.L_E_norm_bound =1;
params.error_norm = 'l2';

params.show_colorbar = 1;
params.colorbar_mode = 'EastOutside';

% the following is required for incremental basis generation

%             coercivity_bound_name: 'fv_coercivity_bound'
%                      alpha_over_k: 1.9915e+007
%                 energy_norm_gamma: 0.5000
%                        axis_equal: 1
%                   error_algorithm: 'fv_error'

params.RB_extension_algorithm= 'RB_extension_PCA_fixspace';
params.RB_stop_timeout= 4*3600;
params.RB_stop_epsilon= 1.0000e-008;
params.RB_stop_Nmax= 150;
params.RB_error_indicator= 'estimator';                
params.RB_generation_mode= 'uniform_fixed';
params.RB_numintervals= [2 2 6];

detailed_data.grid = grid;

switch step
 case 0 % check continuity of dirichlet data 
  x = 0.0;
  y = 0;
  params.decomp_mode = 0; % complete
%  params.affine_decomp_mode = 'complete';
%  params.affine_decomp_mode = 'coefficients';
  udir = zeros(params.nt+1,1);
  dt = params.T/params.nt;
  for ti = 0:params.nt
    params.t = ti*dt;
    udir(ti+1) = dirichlet_values(x,y,params);
%    coeff = dirichlet_values(x,y,params);
%    udir(ti+1) = coeff(3);
  end;
  plot(udir);
 case 1
  disp('performing single detailed simulation')
  U = detailed_simulation(detailed_data.grid, params);
  plot_element_data_sequence(detailed_data.grid,U,params);
 case 2 % dummy reduced basis by single trajectory and simulate
  % turn on matrix inversion for faster basis generation:
  params.skip_matrix_inversion = 0; 
  % single trajectory:

%  disp('detailed simulations');
%  params = params.set_mu([1,1,1],params);
%  U1 =  filecache_function(@detailed_simulation,detailed_data.grid,params);
%  params = params.set_mu([0,1,0],params);
%  U2 =  filecache_function(@detailed_simulation,detailed_data.grid,params);
%  params = params.set_mu([0,0,1],params);
%  U3 =  filecache_function(@detailed_simulation, ...
%                          detailed_data.grid,params);

%A = feval(params.inner_product_matrix_algorithm,detailed_data.grid,params);
%  disp('generating basis');
%  disp('orthonormalization 1:');
%  UON1 = orthonormalize_qr(U1,A);
%  detailed_data.RB= UON1;

  %  disp(['found ',num2str(size(UON1,2)),' basis vectors, only taking ', ...
%       '25']);
%  UON1 = UON1(:,1:25); % only take 30 for each trajectory
%  disp('orthonormalization 2:');
%  % orthonormalize UON2 with respect to UON1:
%  U2P = U2 - UON1 * (UON1' * A * U2);  
%  UON2 = orthonormalize_qr(U2P,A);
%  disp(['found ',num2str(size(UON2,2)),' basis vectors, only taking ', ...
%       '25']);
%  UON2 = UON2(:,1:25); % only take 20 for each trajectory
%  disp('orthonormalization 3:');
%  % orthonormalize UON3 with respect to [UON1, UON2]
%  U3P = U3 - [UON1 UON2] * ([UON1'; UON2'] * A * U3);  
%  UON3 = orthonormalize_qr(U3P,A);
%  disp(['found ',num2str(size(UON3,2)),' basis vectors, only taking ', ...
%       ' 25']);
%  UON3 = UON3(:,1:25); % only take 20 for each trajectory
%  detailed_data.RB = [UON1,UON2,UON3];
%  disp(['found total ',num2str(size(detailed_data.RB,2)),' basis vectors.']);

% Gram Schmidt is much too sensitive for large trajectories
%  UON = orthonormalize_gram_schmidt(U,A);

  detailed_data.grid = grid;
%  detailed_data = rb_basis_generation(detailed_data,params);
%  A = feval(params.inner_product_matrix_algorithm,detailed_data.grid,params);

%  params.set_mu([1,1,1],params);
%  U =  detailed_simulation(detailed_data.grid,params);
%  UON = orthonormalize_qr(U,A);
  
%  detailed_data.RB = UON;
  % adaptive basis generation:
  detailed_data = rb_basis_generation(detailed_data,params);
  disp('saving detailed_data')
  save(detailedfname,'detailed_data','params');
%  params.title = 'detailed simulation result'; 
%  plot_element_data_sequence(detailed_data.grid,U,params);
 case 3
  disp('loading detailed data and parameters from file')
  load(detailedfname);
  disp('reduced simulation:')
  offline_data = rb_offline_prep(detailed_data,params);
%  params.N = size(detailed_data.RB,2)/2;
  params.N = round(size(detailed_data.RB,2));
  reduced_data = rb_online_prep(offline_data, params);
  simulation_data = rb_simulation(reduced_data,params);
  Uappr = rb_reconstruction(detailed_data,simulation_data);  
  params.title = 'reduced simulation result';   
  plot_element_data_sequence(detailed_data.grid,Uappr,params);
  figure, plot(simulation_data.Delta);
  title('error estimator');
 case 4
  disp('loading detailed data and parameters from file')
  load(detailedfname);
  disp('reduced simulation:')
  offline_data = rb_offline_prep(detailed_data,params);
  params.N_mode = 'adapt';
  params.N_adapt_look_ahead_nt = 40;
%  params.N_adapt_slope = 0.5;
  params.N_adapt_slope = 10;
  params.N = size(detailed_data.RB,2); % set max!
%  params.N_adapt_epsilon = 7e-2;
%  params.N_adapt_epsilon = 5e-2;
%  params.N_adapt_epsilon = 3e-2;
%  params.N_adapt_epsilon = 1e-1;
%  params.N_adapt_epsilon = 7e-1;
%  params.N_adapt_epsilon = 1e-3; => everything over 50% N

% schoen niedriges N:
  params.N_mode = 'adapt';
  params = params.set_mu([1 0 0],params);
  params.N_adapt_epsilon = 5e-2;
% schoen adaptives N:
  params.N_mode = 'adapt';
  params = params.set_mu([1 1 0],params);
  params.N_adapt_epsilon = 5e-2;
% schoen adaptives N mit Diffusion:
  params.N_mode = 'adapt';
  params = params.set_mu([1 1 0.25],params);
  params.N_adapt_epsilon = 7e-2;
% N-"switch" zu bestimmten Zeiten: verringerung des Knicks: 
  params.N_mode = 'switch';
  params.N_switch_times =  [0,1];
  params.N_switch_values = [40,64];
% schoen adaptives N mit grosser Diffusion:
  params.N_mode = 'adapt';
  params = params.set_mu([1 1 1],params);
  params.N_adapt_epsilon = 5e-2;
  
  reduced_data = rb_online_prep(offline_data, params);
  simulation_data = rb_simulation(reduced_data,params);
%  Uappr = rb_reconstruction(detailed_data,simulation_data);  
%  params.title = 'reduced simulation result';  
%  plot_element_data_sequence(detailed_data.grid,Uappr,params);


  Nscaled = simulation_data.N /max(simulation_data.N)*...
            params.N_adapt_epsilon;
  figure, plot([simulation_data.Delta;...
                0:(params.N_adapt_epsilon)/params.nt:params.N_adapt_epsilon;...
                Nscaled]');
  title('error estimator, N scaled, target line');
%  figure, plot(simulation_data.N);
%  title('reduced dimension N');
 
%  disp('time measurement: epsilon = 1e-2:')
%  params.N_adapt_epsilon = 1e-2; 
%  tic, simulation_data = rb_simulation(reduced_data,params), toc
%  disp(['average N: ', num2str(mean(simulation_data.N))]);
%  disp('time measurement: epsilon = 1e-5:')
%  params.N_adapt_epsilon = 1e-5; 
%  tic, simulation_data = rb_simulation(reduced_data,params), toc
%  disp(['average N: ', num2str(mean(simulation_data.N))]);
  
 case 5
  disp('loading detailed data and parameters from file')
  load(detailedfname);
  offline_data = rb_offline_prep(detailed_data,params);

  params = params.set_mu([1 1 1],params); % => non-monotonic error curve for Nfixed!
  
  linestyles = {'-.','-','-',':'};  
  linecolors = {[0,0,1],[0,0.5,0],[1,0,0],[0,0,0]};  
  linewidths = [1,1,2,2];
  
  params = params.set_mu([1,1,1],params); % convection 
  params.N_mode = 'adapt';
%  params.N_adapt_epsilon = 7e-1;
  params.N_adapt_epsilon = 7e-2;
  params.N_adapt_look_ahead_nt = 40;
  params.N_adapt_slope = 10;
  disp('reduced_simulation:');
  params.N = 64;
  reduced_data = rb_online_prep(offline_data, params);
  simulation_data = rb_simulation(reduced_data,params);

  N1 = 30;
  N2 = 50;
  params.N_mode = 'fixed';

  params.N = N1;
  reduced_data = rb_online_prep(offline_data, params);
  disp('reduced_simulation:');
  simulation_data_N1 = rb_simulation(reduced_data,params);
  params.N = N2;
  reduced_data = rb_online_prep(offline_data, params);
  disp('reduced_simulation:');
  simulation_data_N2 = rb_simulation(reduced_data,params);
  %Uappr = rb_reconstruction(detailed_data,simulation_data); 
  if isfield(simulation_data,'N'); 
    Ns = simulation_data.N;
  else
    Ns = params.N * ones(1,params.nt+1);
  end;
  figure; % subplot(2,1,2); 
  p = plot([0:params.nt],[ones(size(Ns))*N1;ones(size(Ns))*N2;Ns]');
  for i=1:length(p)
    set(p(i),'Color',linecolors{i},...
             'LineWidth',linewidths(i),...
             'LineStyle',linestyles{i})
  end;
  %set(p(3),'LineWidth',2);
  %set(p(1),'LineStyle','-.');
  title('model dimension N');
  legend(['N-fixed N=',num2str(N1)],...
         ['N-fixed N=',num2str(N2)],'N-adaptive')
  set(gca,'XLim',[0,params.nt]);
  set(gca,'YLim',[0,100]);
  xlabel('time index k');
  saveas(gcf,'Nadaptive_error2_large','epsc');
  set(gcf,'Position',[472   404   564   238]);
  %subplot(2,1,1), 
  figure, p = plot([0:params.nt],[simulation_data_N1.Delta; ...
                    simulation_data_N2.Delta;...
                    simulation_data.Delta;...
                    0:(params.N_adapt_epsilon)/params.nt:params.N_adapt_epsilon
                   ]');
  %keyboard;
  for i=1:length(p)
    set(p(i),'Color',linecolors{i},...
             'LineWidth',linewidths(i),...
             'LineStyle',linestyles{i})
  end;
  %keyboard;
  %set(p(1),'LineStyle','-.');
  %set(p(4),'LineStyle',':');
  %set(p(3),'LineWidth',2);
  title('error estimator evolution');
  l = legend({['N-fixed, N=',num2str(N1)],...
              ['N-fixed, N=',num2str(N2)], ...
              'N-adaptive', ...
              'N-adaptive-target'},'location','northwest');
  set(gca,'XLim',[0,params.nt]);
  xlabel('time index k');
  saveas(gcf,'Nadaptive_Ns2_large','epsc');
  set(gcf,'Position',[472   404   564   238]);
  po = get(l,'Position');
  set(l,'Position',po);
  disp('please save figures as eps (automatic saving does not get aspectratio)')

 case 6
  % time measurement:

  disp('loading detailed data and parameters from file')
  load(detailedfname);
  offline_data = rb_offline_prep(detailed_data,params);

  params = params.set_mu([1 1 1],params); % => non-monotonic error curve for Nfixed!
%  params = params.set_mu([1 1 0.666666],params); % => non-monotonic error curve for Nfixed!
%  params = params.set_mu([1 1 0.2],params); => fixed fastest
%  params = params.set_mu([1 1 0.5],params); % => fixed fastest
%  params = params.set_mu([1 1 0],params); % => almost no time difference!
                                   % diffusivity required!
  params.N = size(detailed_data.RB,2); % set max! 
                                       
  nrep = 10;
  t_exact = zeros(1,nrep);
  disp('detailed_simulation:');
  disp('(skipping detailed computation loop)')
  %  for r = 1:nrep
  for r = 1:1
    tic, U = detailed_simulation(detailed_data.grid,params);
    t_exact(r) = toc;
  end;
  disp('                      time      av(N)     Delta       err ');
  disp(['Exact simulation: ', ...
        num2str(mean(t_exact)), '+-',num2str(std(t_exact)),...
        '   ', num2str(detailed_data.grid.nelements),' 0 0 ']);
  
  disp('reduced_simulations, fixed Ns:');
%  Ns = [100,85,70,55,40,25,10];
%  Ns = [64,56,48,40,32,24,16,8];
  Ns = [32,28,24,20,16,12,8];
  t_Nfix = zeros(length(Ns),nrep);
  Delta_Nfix  = zeros(size(Ns));
  % prevent matrix inversion for fair comparison
  params.skip_matrix_inversion = 1; 
  params.N_mode = 'fixed';
  
  %dummy simulation
  dummy_reduced_data = rb_online_prep(offline_data, params); 
  dummy = rb_simulation(dummy_reduced_data,params);
  
  for r = 1:nrep
    for n = 1:length(Ns) 
      simulation_data = [];
      params.N = Ns(n);
      reduced_data = rb_online_prep(offline_data, params); 
      tic;
      simulation_data = rb_simulation(reduced_data,params);
      t_Nfix(n,r)= toc;
      Uappr = rb_reconstruction(detailed_data,simulation_data); 
      err = fv_error(Uappr,U,detailed_data.grid,params);  
      err_Nfix(n) = max(err);
      Delta_Nfix(n) = simulation_data.Delta(end);
      %  disp(['N = ',num2str(Ns(n)),'  ',num2str(t_Nfix(n,r)),'    ',...
      % num2str(Ns(n)),'    ',...
      % num2str(Delta_Nfix(n)),'    ',num2str(err_Nfix(n))]);
    end;
  end;
  disp('overall:')
  for n = 1:length(Ns) 
    disp(['N = ',num2str(Ns(n)),'  ',num2str(mean(t_Nfix(n,:))),...
          '+-',num2str(std(t_Nfix(n,:))),'   ',...
          num2str(Ns(n)),'    ',...
          num2str(Delta_Nfix(n)),'    ',num2str(err_Nfix(n))]);
  end;
  
  disp('reduced_simulations, adaptive Ns:');
%  epsilons = [5e-6,1e-5,5e-5,1e-4,5e-4,1e-3];
%  epsilons = [1e-2,2e-2,5e-2,1e-1,2e-1,5e-1,1,5];
  epsilons = [1e-1,2e-1,5e-1,7e-1,1,2,5];
  Delta_Nadaptive  = zeros(size(epsilons));
  t_Nadaptive = zeros(length(epsilons),nrep);
  Nav_Nadaptive = zeros(size(epsilons));
  params.N = size(detailed_data.RB,2);
  reduced_data = rb_online_prep(offline_data, params); 
  params.N_mode = 'adapt';
%  params.N_adapt_look_ahead_nt = 20;
  params.N_adapt_look_ahead_nt = 40;
  %  params.N_adapt_slope = 0.5;
%  params.N_adapt_slope = 5;
  params.N_adapt_slope = 10;
  
  for r = 1:nrep
    for e = 1:length(epsilons) 
      simulation_data =[];
      params.N_adapt_epsilon = epsilons(e);
      tic;
      simulation_data = rb_simulation(reduced_data,params);
      t_Nadaptive(e,r)= toc;
      Uappr = rb_reconstruction(detailed_data,simulation_data); 
      err = fv_error(Uappr,U,detailed_data.grid,params);  
%      err_Nadaptive(e) = err(end);
      err_Nadaptive(e) = max(err);
      Delta_Nadaptive(e) = simulation_data.Delta(end);
      Nav_Nadaptive(e) = mean(simulation_data.N);
      %  disp(['epsilon = ',num2str(epsilons(e)),...
      % '  ',num2str(t_Nadaptive(e,r)),'    ',num2str(Nav_Nadaptive(e)),...
      % '    ',...
      % num2str(Delta_Nadaptive(e)),'    ',num2str(err_Nadaptive(e))]);
    end;
  end;
  disp('overall:')
  for e = 1:length(epsilons) 
    disp(['epsilon = ',num2str(epsilons(e)),...
          '  ',num2str(mean(t_Nadaptive(e,:))),'+-',...
          num2str(std(t_Nadaptive(e,:))),'    ',num2str(Nav_Nadaptive(e)),...
          '    ',...
          num2str(Delta_Nadaptive(e)),'    ',num2str(err_Nadaptive(e))]);
  end;
  disp(['above results averaged over ',num2str(nrep),' runs']);

  
%    disp('reduced_simulations, switching Ns:');
%%  epsilons = [5e-6,1e-5,5e-5,1e-4,5e-4,1e-3];
%  switch_times = 0:0.25:2;
%  Delta_Nswitch  = zeros(size(switch_times));
%  t_Nswitch = zeros(length(switch_times),nrep);
%  Nav_Nswitch = zeros(size(switch_times));
%  params.N = size(detailed_data.RB,2);
%  reduced_data = rb_online_prep(offline_data, params); 
%  params.N_mode = 'switch';
%  params.N_switch_values = [40,64];
%  params.N_switch_times = [0,1];
%  
%  for r = 1:nrep
%    for e = 1:length(switch_times) 
%      simulation_data =[];
%      params.N_switch_times(2) = switch_times(e);
%      tic;
%      simulation_data = rb_simulation(reduced_data,params);
%      t_Nswitch(e,r)= toc;
%      Uappr = rb_reconstruction(detailed_data,simulation_data); 
%      err = fv_error(Uappr,U,detailed_data.grid,params);  
%%      err_Nadaptive(e) = err(end);
%      err_Nswitch(e) = max(err);
%      Delta_Nswitch(e) = simulation_data.Delta(end);
%      Nav_Nswitch(e) = mean(simulation_data.N);
%      %  disp(['epsilon = ',num2str(epsilons(e)),...
%      %        '  ',num2str(t_Nadaptive(e,r)),'    ',num2str(Nav_Nadaptive(e)),...
%      %        '    ',...
%      %        num2str(Delta_Nadaptive(e)),'    ',num2str(err_Nadaptive(e))]);
%    end;
%  end;
%  disp('overall:')
%  for e = 1:length(switch_times) 
%    disp(['switch_time = ',num2str(switch_times(e)),...
%         '  ',num2str(mean(t_Nswitch(e,:))),'+-',...
%         num2str(std(t_Nswitch(e,:))),'    ',num2str(Nav_Nswitch(e)),...
%         '    ',...
%         num2str(Delta_Nswitch(e)),'    ',num2str(err_Nswitch(e))]);
%  end;
%  disp(['above results averaged over ',num2str(nrep),' runs']);

  linestyles = {'-.','-','-',':'};  
  linecolors = {[0,0,1],[0,0.5,0],[1,0,0],[0,0,0]};  
  linewidths = [1,1,2,2];
  markers = {'','x','d',''};
  
  figure
  plot(mean(t_Nadaptive,2),Delta_Nadaptive',...
       'Color',linecolors{3},...
       'LineStyle',linestyles{3},'Marker',markers{3},'LineWidth',linewidths(3));
  hold on;
  plot(mean(t_Nfix,2),Delta_Nfix',...
              'Color',linecolors{2},...
       'LineStyle',linestyles{2},'Marker',markers{2},'LineWidth',linewidths(2));
  hold on;
%  plot(mean(t_Nswitch,2),Delta_Nswitch','r-');
%  hold on;
  xlabel('runtime')
  ylabel('error estimator')
  set(gca,'Yscale','log')
  title('error estimator over runtime')
%  legend('N-adaptive','N-fixed','N-switch')
  legend('N-adaptive','N-fixed')
  
  figure
  plot(mean(t_Nadaptive,2),err_Nadaptive',...
              'Color',linecolors{3},...
       'LineStyle',linestyles{3},'Marker',markers{3},'LineWidth',linewidths(3));
  hold on;
  plot(mean(t_Nfix,2),err_Nfix',....
                     'Color',linecolors{2},...
       'LineStyle',linestyles{2},'Marker',markers{2},'LineWidth',linewidths(2));
%  hold on;
%  plot(mean(t_Nswitch,2),err_Nswitch','r-');
  xlabel('runtime')
  ylabel('error')
  set(gca,'Yscale','log')
  title('error over runtime')
%  legend('N-adaptive','N-fixed','N-switch')
  legend('N-adaptive','N-fixed')

%  figure
%  plot(Nav_Nadaptive,err_Nadaptive',...
%              'Color',linecolors{3},...
%       'LineStyle',linestyles{3},'Marker',markers{3});
%  hold on;
%  plot(Ns,err_Nfix',...
%                     'Color',linecolors{2},...
%       'LineStyle',linestyles{2},'Marker',markers{2});
%%  hold on;
%%  plot(Nav_Nswitch,err_Nswitch','r-');
%  xlabel('average N')
%  ylabel('error')
%  set(gca,'Yscale','log')
%  title('error over average N')
%%  legend('N-adaptive','N-fixed','N-switch')
%  legend('N-adaptive','N-fixed')

  keyboard;
  
  disp('please save figures as eps and fig by hand.')
  
 otherwise
  error('step number unknown!');  
end;

return;

%| \docupdate