riccati_rb_gamma.m

Go to the documentation of this file.

function gamma = riccati_rb_gamma(model, reduced_data, sim_data)
%[gamma] = riccati_rb_gamma(model, reduced_data, sim_data)
%
% Calculate the norm of the inverse Lyapunov operator for the RB
% approximation. This can be done in several ways, depending on the field
% reduced_data.estim.gamma.mode which is directly copied from
% model.RB_gamma_mode. 
%
%    interpol: use an interpolation based calculation of the constant
%    kernel_interpolation: kernel interpolation based approximation
%    exponential_bound: use an exponential bound to calculate the constant
%
% Andreas Schmidt, 2015
if ~isfield(reduced_data.estim, 'gamma')
    gamma = 1;
    return
end

if strcmp(reduced_data.estim.gamma.mode, 'interpol')
    gamma = reduced_data.estim.gamma.interpolant(get_mu(model)');
elseif strcmp(reduced_data.estim.gamma.mode, 'kernel_interpolation')
    % Get the kernel weights and the kernel
    alpha = reduced_data.estim.gamma.alpha;
    k = riccati_kernel_functions(model, reduced_data.estim.gamma.kernel);
    mus = reduced_data.estim.gamma.interpol_mus;
    
    gamma = 0;
    mu = get_mu(model);
    
    for j = 1:length(alpha)
        gamma = gamma + alpha(j)*k(mu(:)', mus(j,:));
    end
elseif strcmp(reduced_data.estim.gamma.mode, 'exponential_bound')
    % TODO: Make me fast!!!
    model_data = gen_model_data(model);
    [E,A,B,~] = riccati_assemble_system_matrices(model, model_data);
    
    nu = eigs((A+A'), E, 1, 'sm')/2;
    n = normest(E\(B*B'*reduced_data.RB*sim_data.PN*reduced_data.RB'*E));

    gamma = -reduced_data.estim.normEInv^2/2/(nu+n);
else
    gamma = -1;
    warning('The gamma approximation mode is not known')
end

end