scm_demo.m

Go to the documentation of this file.

% scm_demo.m
% a simple demo script which produces scm_offline_data for the inf-sup
% constant of the scm_minimal_model. Then for a fine set of parameters in
% [0,1] the exact constant beta(mu), the lower bound beta_{LB}(mu) and
% the exact cercivity constant alpha(mu) are computet and the results are
% plottet.
% This shows that the scm_minimal_model is coercive up to mu = 0.5 and from
% there on it is only inf-sup stable. For mu = 0.5 it is neither cercive
% nor inf-sup stable. The main point is to show, that the SCM is working
% fine for little 'lin_stat' inf-sup stable problems.

% Dominik Garmatter 20.09 2012

clear;
% generate model and model_data
size = 100;
model = scm_minimal_model(size);
model_data = gen_model_data(model);
K = model.get_inner_product_matrix(model_data);
% add components to data
old_decomp_mode = model.decomp_mode;
model.decomp_mode = 1;
A_comp = model.operators(model, model_data);
model.decomp_mode = old_decomp_mode;
model_data.A_comp = A_comp;
% scm settings
model.scm_eps_tol = eps;
model.scm_size_C = 20;
model.scm_M_alpha = 20;
model.scm_M_plus = 200;
% generate offline_data for beta
model.scm_desired_constant = 2;
scm_offline_data_beta = scm_offline(model, model_data);
% write into reduced data
reduced_data_beta = [];
reduced_data_beta.scm_offline_data = scm_offline_data_beta;
% for every mu get the direct constant, and the lower bound
mu_fineness = 99;
gather_alpha = zeros(mu_fineness+1,1);
gather_beta = zeros(mu_fineness+1,1);
gather_beta_LB = zeros(mu_fineness+1,1);
for i = 1 : mu_fineness+1
    fprintf('.');
        model = set_mu(model, 0 + (i-1)/mu_fineness); % realises discrete [0,1]
    model.decomp_mode = 2;
    A_coeff = model.operators(model, model_data);
    A = lincomb_sequence(A_comp, A_coeff);
    % also solve the EWP for the coercivity constant
    LM_alpha = eigs(0.5*(A+A'),K,1,'LM');
    if LM_alpha >= 0
        gather_alpha(i) = eigs(0.5*(A+A') - (LM_alpha+1e-10)*K,K,1,'LM') + (LM_alpha+1e-10);
    else
        gather_alpha(i) = LM_alpha;
    end
    LM_beta = sqrt(eigs(A*K*A',K,1,'LM'));
    if LM_beta >= 0
        gather_beta(i) = sqrt(eigs(A*K*A' - (LM_beta+1e-10)*K,K,1,'LM') + (LM_beta+1e-10)); % inv(K) == K (normally it would be A*inv(K)*A')
    else
        gather_beta(i) = LM_beta;
    end
    gather_beta_LB(i) = scm_lower_bound(model, reduced_data_beta);
end
fprintf('\n');
disp('plotting...')
% plot
f = figure;
subplot(1,2,1)
plot(0:1/mu_fineness:1,gather_alpha,'r')
xlabel('\mu')
leg = legend('\alpha','\alpha_{LB}');
set(leg,'Location','Best')
subplot(1,2,2)
plot(0:1/mu_fineness:1,gather_beta,'r:',0:1/mu_fineness:1,gather_beta_LB,'g--')
xlabel('\mu')
leg = legend('\beta','\beta_{LB}');
set(leg,'Location','Best')
axes('position',[0,0,1,1],'visible','off');
text(0.15,0.97,['This plot shows that the scm\_minimal\_model is coercive up to \mu = ' mat2str(0.5) ' and always inf-sup stable (except for \mu = ' mat2str(0.5) '). The main point is that the scm is working fine for small lin\_stat inf-sup stable problems.'])
set(gcf,'Color','white')