eop_C_L_I.m

Go to the documentation of this file.

function model = eop_C_L_I(model)
%model = eop_C_L_I(model)
%
% function to compute and empirical continuousity constant of the Matrix 
% L_I^-1 of a linear evolution scheme. Required for the nomal L2-error
% estimator. 

% computes the largest magnitude of ||L_I^-1|| for all \mu\in M_train and 
% takes the maximum of those as empirical continuousity constant.

% Dominik Garmatter 11.03 2012

model_data = gen_model_data(model);

% get the components of L_I
model.decomp_mode = 1;
[L_I_comp, ~, ~] = model.operators_ptr(model, model_data);

% formt he Training set Mtrain as a discretisation of the parameter space
par.numintervals = model.RB_numintervals-2;
par.range = model.mu_ranges;
diskret_par_space = cubegrid(par);
Mtrain = get(diskret_par_space, 'vertex')';

spectral_build = zeros(size(Mtrain,2),1);
n = size(L_I_comp{1},1);

% build vector
for i = 1:size(Mtrain,2)
    disp(i)
    % assemble L_I for the current mu
    model = model.set_mu(model,Mtrain(:,i));
    model.decomp_mode = 2;
    [L_I_coeff, ~, ~] = model.operators_ptr(model, model_data);
    L_I = lincomb_sequence(L_I_comp,L_I_coeff);
    % compute the largest magnitude
    spectral_build(i) = eigs(@(x)my_mult_for_L_I(x,L_I,model_data.W),n,1,'LM');
    try %first try to solve the EWP with default tolerance eps
      spectral_build(i) = eigs(@(x)my_mult_for_L_I(x, L_I, model_data.W), n, model_data.W, 1,'LM');
    catch err 
      if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) %if no eigenvalues to sufficient accuracy was found reduce the tolerance and try it again
          opts.tol = 1e-14;
          no_solution = 1;
          while no_solution
            disp(['I reduced the tolerance of an eigenvalueproblem to ' mat2str(opts.tol)]);
            try
              spectral_build(i) = eigs(@(x)my_mult_for_L_I(x, L_I, model_data.W), n, model_data.W, 1,'LM',opts);
              no_solution  = 0; %end when there is no error
            catch err2 %when it still fails reduce the tolerance and try again!
              no_solution = 1;
              if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14'))  
                opts.tol = 100 * opts.tol;
              end
              if opts.tol > 1e-3 %if the tolerance is reduced until a certain value, give an error since the tolerance is too soft!
                  error(['I reduced the tolerance to ' mat2str(opts.tol) ' but I still could not find a solution. Reconsider your Problem']);
              end
            end
          end
      else %if there is another error rethrow
          rethrow(err)
      end
    end
end
model.L_I_inv_norm_bound = max(spectral_build);


%==========================================================================


function y = my_mult_for_L_I(x,L,K)
% AFUN for eigs to compute ||(L_I)^-1(\mu)||_2.
% Does the operationy = K^-1/2 * L^-1 * K * L^-1 * K^-1/2 * x.
Kinv = sqrt(inv(K));
z1 = Kinv * x;
z2 = L\z1;
z3 = K*z2;
z4 = L\z3;
z5 = Kinv * z4;
y = z5;