eop_fd_functionals.m

Go to the documentation of this file.

function [v, l] = eop_fd_functionals(model, model_data)
%[v, l] = eop_fd_functionals(model, model_data)
%
% function which computes the riesz-representant v of the functionals, so
% that one can evaluate them via
% l(P^M;\mu)=v'*sim_data.U
% (see eop_theta_functional for the 'theta'-case) 
% if only v is desired as output, the function is affinely decomposed.
%
% Also the continuousity constant l of the functionals can be computed and
% given as output. If v and l isare desired as output no affine decomposition
% is available (and not neccessary).

% In the following \Omega_0 is always a subset \Omega which has to match
% the discretization.

% possible functionals (i.e. strings for model.name_output_functional = '')
% - 'average price': l(P^M;\mu) = 1/|\Omega_0|*\int_{\Omega_0}P^M(x)dx
% -  'theta': l(P^M;\mu) = -1/|\Omega_0|*int_{\Omega}\frac{\partial
%           P^M}{\partial t}(x) dx
% - 'average dP/dS1': l(P^M;\mu) =
%                           1/|\Omega_0|*\int_{\Omega_0}\frac{\parial P^M}{\partial S_1}(x)dx
% - 'average dP/dS2': l(P^M;\mu) =
%                           1/|\Omega_0|*\int_{\Omega_0}\frac{\parial P^M}{\partial S_2}(x)dx

% Dominik Garmatter 25.07 2012


% the evaluation of all functionals is done on a subset of Omega. Here we
% check if the subset defined by
% [model.functional_subset_S1]x[model.functional_subset_S2] is conform with
% the discretization and Omega.
x1 = model.functional_subset_S1(1)/model.h1;
x1_check = floor(model.functional_subset_S1(1)/model.h1);
x2 = model.functional_subset_S1(2)/model.h1;
x2_check = floor(model.functional_subset_S1(2)/model.h1);
y1 = model.functional_subset_S2(1)/model.h2;
y1_check = floor(model.functional_subset_S2(1)/model.h2);
y2 = model.functional_subset_S2(2)/model.h2;
y2_check = floor(model.functional_subset_S2(2)/model.h2);
if x1 ~= x1_check || x2 ~= x2_check || y1 ~= y1_check || y2 ~= y2_check
  error('subset input inconsistent with discretisation!');
end
if x1*model.h1 < model.xrange(1) || x2*model.h1 > model.xrange(2)...
   || y1*model.h2 < model.yrange(1) || y2*model.h2 > model.yrange(2)
  error('subset does not fit the grid!');
end

if nargout < 2 % only v is computed, affinely decomposed
  if ~isfield(model,'decomp_mode')
    model.decomp_mode = 0;
  end  
    
  if model.decomp_mode == 2; %no high-dimensional data available
    v = 1; % coefficient of the riesz-representant is always 1
    
  else % high-dimensional data is available, i.e. if decomp-mode is 1 or 0
    % define the subset and the number of gridpoints in it
    H = model_data.grid.nvertices;
    subset_S1 = model.functional_subset_S1;
    subset_S2 = model.functional_subset_S2;
    N1_subset = (subset_S1(2) - subset_S1(1))/model.h1 + 1;
    N2_subset = (subset_S2(2) - subset_S2(1))/model.h2 + 1;
    H_subset  = N1_subset*N2_subset;
    % the normal set is [0xS1bar]x[0xS2bar] but the in indexing in matlab
    % starts with 1, so we move the subset to match this
    new_subset_S1 = [subset_S1(1)/model.h1 + 1 subset_S1(2)/model.h1 + 1];
    new_subset_S2 = [subset_S2(1)/model.h2 + 1 subset_S2(2)/model.h2 + 1];

  switch model.name_output_functional

    case 'average price'

        % riesz-representatn is 1 in the subset and 0 else
        if model.functional_subset_S1(1) == model.xrange(1) || model.functional_subset_S1(2) == model.xrange(2)...
            || model.functional_subset_S2(1) == model.yrange(1) || model.functional_subset_S2(2) == model.yrange(2)
            v_build = ones(H,1);
        else
            v_build = zeros(H,1);
            for i = 1:N2_subset
                v_build((new_subset_S2(1)-1+(i-1))*model.N1+new_subset_S1(1) : (new_subset_S2(1)-1+(i-1))*model.N1+new_subset_S1(2)) = 1;
            end
        end
        v = 1/H_subset*v_build; % scale with the respective factors
        if model.decomp_mode == 1 % cell output in component case
          v = {v};
        end
        
    case 'theta'
        
        % riesz-representatn is 1 in the subset and 0 else
        if model.functional_subset_S1(1) == model.xrange(1) || model.functional_subset_S1(2) == model.xrange(2)...
                || model.functional_subset_S2(1) == model.yrange(1) || model.functional_subset_S2(2) == model.yrange(2)
            v_build = ones(H,1);
        else
            v_build = zeros(H,1);
            for i = 1:N2_subset
                v_build((new_subset_S2(1)-1+(i-1))*model.N1+new_subset_S1(1) : (new_subset_S2(1)-1+(i-1))*model.N1+new_subset_S1(2)) = 1;
            end
        end
        v = -1/(H_subset*model.deltaT)*v_build; % scale with the respective factors
        if model.decomp_mode == 1 % cell output in component case
          v = {v};
        end
        
    case 'average dP/dS1'

        % in case of dP/dS1 the central difference quotient is used for the
        % inner points and the left (or right) difference quotient is used
        % at the left (or right) boundary of the subset. Therefore many
        % terms are 0 and this remains
        v_build = zeros(model.N1,1);
        if N1_subset == 1
            error('subset too small for a partial derivative in S1 direction!');
        elseif N1_subset == 2
            v_build(new_subset_S1(1) : new_subset_S1(2)) = [-2, 2];
        elseif N1_subset == 3
            v_build(new_subset_S1(1) : new_subset_S1(2)) = [-1.5, 0, 1.5];
        else
            v_build(new_subset_S1(1))   = -1.5;
            v_build(new_subset_S1(1)+1) = 0.5;
            v_build(new_subset_S1(2)-1) = -0.5;
            v_build(new_subset_S1(2))   = 1.5;
        end
        v_build_2 = zeros(H,1);
        % write those entrys at the positions of the left (or right) boundary of the subset
        v_build_2((new_subset_S2(1)-1)*model.N1+1:new_subset_S2(2)*model.N1) = repmat(v_build,1,N2_subset);
        v = 1/(H_subset*model.h1)*v_build_2; % scale with the respective factors
        if model.decomp_mode == 1 %cell output in component case
          v = {v};
        end
        
    case 'average dP/dS2'

        % in case of dP/dS2 the central difference quotient is used for the
        % inner points and the left (or right) difference quotient is used
        % at the lower (or upper) boundary of the subset. Therefore many
        % terms are 0 and this remains
        v_build = zeros(H,1);
        if N2_subset == 1
            error('subset too small for a partial derivative in S2 direction!');
        elseif N2_subset == 2
            v_build((new_subset_S2(1)-1)*model.N1+new_subset_S1(1) : (new_subset_S2(1)-1)*model.N1+new_subset_S1(2)) = -2;
            v_build(new_subset_S2(1)*model.N1+new_subset_S1(1) : new_subset_S2(1)*model.N1+new_subset_S1(2))         = 2;
        elseif N2_subset == 3
            v_build((new_subset_S2(1)-1)*model.N1+new_subset_S1(1) : (new_subset_S2(1)-1)*model.N1+new_subset_S1(2)) = -1.5;
            v_build((new_subset_S2(1)+1)*model.N1+new_subset_S1(1) : (new_subset_S2(1)+1)*model.N1+new_subset_S1(2)) = 1.5;
        else
            v_build((new_subset_S2(1)-1)*model.N1+new_subset_S1(1) : (new_subset_S2(1)-1)*model.N1+new_subset_S1(2)) = -1.5;
            v_build(new_subset_S2(1)*model.N1+new_subset_S1(1) : new_subset_S2(1)*model.N1+new_subset_S1(2))         = 0.5;
            v_build((new_subset_S2(2)-2)*model.N1+new_subset_S1(1) : (new_subset_S2(2)-2)*model.N1+new_subset_S1(2)) = -0.5;
            v_build((new_subset_S2(2)-1)*model.N1+new_subset_S1(1) : (new_subset_S2(2)-1)*model.N1+new_subset_S1(2)) = 1.5;
        end
        v = 1/(H_subset*model.h2)*v_build; % scale with the respective factors
        if model.decomp_mode == 1 %cell output in component case
          v = {v};
        end
        
    otherwise
        error('type of output functional is unknown!');
  end %end of cases if nargout < 2
  end %end of decomp_mode

else %nargout = 2 case, i.e. compute both v and l - NOT affinely decomposed
    % define the subset and the number of gridpoints in it
    H = model_data.grid.nvertices;
    subset_S1 = model.functional_subset_S1;
    subset_S2 = model.functional_subset_S2;
    N1_subset = (subset_S1(2) - subset_S1(1))/model.h1 + 1;
    N2_subset = (subset_S2(2) - subset_S2(1))/model.h2 + 1;
    H_subset  = N1_subset*N2_subset;
    % the normal set is [0xS1bar]x[0xS2bar] but the in indexing in matlab
    % starts with 1, so we move the subset to match this
    new_subset_S1 = [subset_S1(1)/model.h1 + 1 subset_S1(2)/model.h1 + 1];
    new_subset_S2 = [subset_S2(1)/model.h2 + 1 subset_S2(2)/model.h2 + 1];

  switch model.name_output_functional

    case 'average price'

        % riesz-representatn is 1 in the subset and 0 else
        if model.functional_subset_S1(1) == model.xrange(1) || model.functional_subset_S1(2) == model.xrange(2)...
            || model.functional_subset_S2(1) == model.yrange(1) || model.functional_subset_S2(2) == model.yrange(2)
            v_build = ones(H,1);
        else
            v_build = zeros(H,1);
            for i = 1:N2_subset
                v_build((new_subset_S2(1)-1+(i-1))*model.N1+new_subset_S1(1) : (new_subset_S2(1)-1+(i-1))*model.N1+new_subset_S1(2)) = 1;
            end
        end
        v = 1/H_subset*v_build; % scale with the respective factors
        l = sqrt(H/H_subset); % continuousity constant
        
    case 'theta'
        
        % riesz-representatn is 1 in the subset and 0 else
        if model.functional_subset_S1(1) == model.xrange(1) || model.functional_subset_S1(2) == model.xrange(2)...
                || model.functional_subset_S2(1) == model.yrange(1) || model.functional_subset_S2(2) == model.yrange(2)
            v_build = ones(H,1);
        else
            v_build = zeros(H,1);
            for i = 1:N2_subset
                v_build((new_subset_S2(1)-1+(i-1))*model.N1+new_subset_S1(1) : (new_subset_S2(1)-1+(i-1))*model.N1+new_subset_S1(2)) = 1;
            end
        end
        v = -1/(H_subset*model.deltaT)*v_build; % scale with the respective factors
        l = 1/(model.deltaT)*sqrt(H/H_subset); % continuousity constant
        
    case 'average dP/dS1'

        % in case of dP/dS1 the central difference quotient is used for the
        % inner points and the left (or right) difference quotient is used
        % at the left (or right) boundary of the subset. Therefore many
        % terms are 0 and this remains
        v_build = zeros(model.N1,1);
        if N1_subset == 1
            error('subset too small for a partial derivative in S1 direction!');
        elseif N1_subset == 2
            v_build(new_subset_S1(1) : new_subset_S1(2)) = [-2, 2];
            l_build = 8; % equals v_build'*v_build
        elseif N1_subset == 3
            v_build(new_subset_S1(1) : new_subset_S1(2)) = [-1.5, 0, 1.5];
            l_build = 4.5; % equals v_build'*v_build
        else
            v_build(new_subset_S1(1))   = -1.5;
            v_build(new_subset_S1(1)+1) = 0.5;
            v_build(new_subset_S1(2)-1) = -0.5;
            v_build(new_subset_S1(2))   = 1.5;
            l_build = 5; % equals v_build'*v_build
        end
        v_build_2 = zeros(H,1);
        % write those entrys at the positions of the left (or right) boundary of the subset
        v_build_2((new_subset_S2(1)-1)*model.N1+1:new_subset_S2(2)*model.N1) = repmat(v_build,1,N2_subset);
        v = 1/(H_subset*model.h1)*v_build_2; % scale with the respective factors
        l = H/(model.h1*H_subset)*sqrt(N2_subset*1/H*l_build); % continuousity constant
        
    case 'average dP/dS2'

        % in case of dP/dS2 the central difference quotient is used for the
        % inner points and the left (or right) difference quotient is used
        % at the lower (or upper) boundary of the subset. Therefore many
        % terms are 0 and this remains
        v_build = zeros(H,1);
        if N2_subset == 1
            error('subset too small for a partial derivative in S2 direction!');
        elseif N2_subset == 2
            v_build((new_subset_S2(1)-1)*model.N1+new_subset_S1(1) : (new_subset_S2(1)-1)*model.N1+new_subset_S1(2)) = -2;
            v_build(new_subset_S2(1)*model.N1+new_subset_S1(1) : new_subset_S2(1)*model.N1+new_subset_S1(2))         = 2;
            l_build = 8; % equals v_build'*v_build
        elseif N2_subset == 3
            v_build((new_subset_S2(1)-1)*model.N1+new_subset_S1(1) : (new_subset_S2(1)-1)*model.N1+new_subset_S1(2)) = -1.5;
            v_build((new_subset_S2(1)+1)*model.N1+new_subset_S1(1) : (new_subset_S2(1)+1)*model.N1+new_subset_S1(2)) = 1.5;
            l_build = 4.5; % equals v_build'*v_build
        else
            v_build((new_subset_S2(1)-1)*model.N1+new_subset_S1(1) : (new_subset_S2(1)-1)*model.N1+new_subset_S1(2)) = -1.5;
            v_build(new_subset_S2(1)*model.N1+new_subset_S1(1) : new_subset_S2(1)*model.N1+new_subset_S1(2))         = 0.5;
            v_build((new_subset_S2(2)-2)*model.N1+new_subset_S1(1) : (new_subset_S2(2)-2)*model.N1+new_subset_S1(2)) = -0.5;
            v_build((new_subset_S2(2)-1)*model.N1+new_subset_S1(1) : (new_subset_S2(2)-1)*model.N1+new_subset_S1(2)) = 1.5;
            l_build = 5; % equals v_build'*v_build
        end
        v = 1/(H_subset*model.h2)*v_build; % scale with the respective factors
        l = H/(model.h2*H_subset)*sqrt(N1_subset*1/H*l_build); % continuousity constant
        
    otherwise
        error('type of output functional is unknown!');
  end %end of nargout = 2 case
      
end %end of nargout if-cases