scm_coercive_demo.m

Go to the documentation of this file.

function scm_coercive_demo()
%scm_coercive_demo()
% 
% Demo of SCM (coercive constant)
% 

% David Kreplin 05.01.2015

clear;

model = thermalblock_model_hp();
model.gen_model_data=@lin_stat_gen_model_data;
model_data = model.gen_model_data(model);

operator=@ (model) model.operators(model,model_data);
innerproduct_matrix = model.get_inner_product_matrix(model_data);
% SCM properties
params.Ma = inf;
params.Mp = 0;
params.NumberOfParameters = 1000;
params.eps=0.15;
params.max_iterations=100;
params.parallel=0;
params.dirichlet_gids = model_data.df_info.dirichlet_gids;
params.type='coercive';



%OFFLINE
SCM_data=SCM(model,params,operator,innerproduct_matrix);
%Alternative (see documentation): SCM_data=SCM(model,params,model_data);

rmserror=0;
maderror = 0;
innerproduct_matrix = model.get_inner_product_matrix(model_data);
innerproduct_matrix = remove_dirichlet(innerproduct_matrix,0,params.dirichlet_gids);
fprintf('Results:\ncoercive constant      upper bound       lower bound (scm coercive constant)\n');
for i = 1:20
    mu = rand_uniform(1,model.mu_ranges);
    model = set_mu(model,mu);
    model.decomp_mode=0;
    [A,~] =  model.operators(model,model_data);
     A=remove_dirichlet(A,0,params.dirichlet_gids);
    coer=eigs(A,innerproduct_matrix,1,'sm');
    alpha_UB =  SCMonline(mu,model,operator,SCM_data,'ub');
    %ONLINE:
    alpha_LB = SCMonline(mu,model,operator,SCM_data,'lb');
    %alpha_LB = SCMonline(mu,model,model_data,SCM_data);
    
    fprintf('%f               %f          %f\n',coer,alpha_UB,alpha_LB);
    
    rmserror = rmserror+(coer-alpha_LB)^2;
    maderror = maderror + abs(coer-alpha_LB);
end
rmserror = sqrt(rmserror/20)
maderror = maderror/20



keyboard;
end


function A=remove_dirichlet(A,mode,I)
if mode == 0
    A(I,:) = [];
    A(:,I)= [];
else
    for i = 1:length(A)
        A{1,i}(I,:) = [];
        A{1,i}(:,I) = [];
    end
end
end