rbf_elliptic.m

Go to the documentation of this file.

function res = rbf_elliptic(step,params);
%function res = rbf_elliptic(step,params);
% function for meshless collocation of fem-problem. step can be
% used to choose several cases. step=1 does not need to be called
% directly but is used by the other steps. So call step=2, 3 and
% change options there.
%
% step 1: single solve. params.gamma can be set,
%         params.colpoints_choice can be set to 1(random) or 2(structured).  
%         params.model_choice = 1 (poisson) or 2 (singular elliptic)
%         If params.no_plots is set, then no plots
%         and no screen output are performed. Output fields are 
%         res.u_h and the l2-error res.err and condition number res.cond
%
% step 2 % working example for poisson
% 
% step 3 % working example for difficult ellipic model
%  
% step 4 % determine errors and matrix conditions over a range of gammas
%          and plot of error by repeated calls of step 1. Finds optimal gamma

res = [];

if (nargin < 1) | (isempty(step))
  step = 2; % default: show working example of poisson model
end;

if (nargin < 2) 
  params = [];
end;

switch step
 case 1
   % check if plots are required
  if isfield(params,'no_plots')
    no_plots = params.no_plots;
  else % default: plotting on
    no_plots = 0;
  end;
  
  % set gamma of rbf-kernel
  if isfield(params,'gamma')
    gamma = params.gamma;
  else % default
    gamma = 50;
%    gamma = 0.01;
  end;

  % set gamma of rbf-kernel
  if isfield(params,'colpoints_choice')
    colpoints_choice = params.colpoints_choice;
  else
    colpoints_choice = 2;
  end;

  % set model type 
  if isfield(params,'model_choice')
    model_choice = params.model_choice;
  else
    model_choice = 1; % poisson
    %model_choice = 2; % elliptic
  end;
  
  % set refinement steps 
  if isfield(params,'refinements')
    refinements = params.refinements;
  else
    refinements = 2; 
  end;

  if model_choice == 1
    model = praktikum_poisson_model(params);
  elseif model_choice == 2
    model = praktikum_ellipt_model(params);
  end;
  
  % scalar argument and scalar valued rbf kernel and derivatives
  % function can be called in many points s simultaneous.
  % s is scalar (vectorial if many points)
  kernel = @(s) exp(-gamma * s.^2);
  kernel_derivative = @(s) -2*gamma * s.*exp(-gamma * s.^2);
  %kernel_2nd_derivative = @(s) -2*gamma*exp(-gamma * s.^2).* ...
  %    (- 2 *gamma * s + 1);
  % numerical gradient determination by central difference of a function f
  % x is 1x2-vec (or matrix n x 2 if many points)
  %ds = (1e-3)*gamma;
%  num_gradient = @(f,x) [(f([x(:,1),x(:,2)+ds/2]) ...
%                         -f([x(:,1),x(:,2)-ds/2]))/ds, ...
%                   (f([x(:,1)+ds/2,x(:,2)]) ...
%                    -f([x(:,1)-ds/2,x(:,2)]))/ds];

% collocation centers for unit square
switch colpoints_choice
  case 1 % random points for unit square
    if model_choice ==2
      error('random points not implemented for l-shape domain')
    end;
    rand('seed',101); % initialize random number generator
    n_int = 15 * 2^refinements;
    n_dir = 8 * 2^refinements; % dividable by 4
    n_neu = 0;
    X_int = rand(n_int,2);
    x_dir = rand(n_dir/4,1);
    x_zero = zeros(size(x_dir,1),1);
    x_ones = ones(size(x_dir,1),1);
    
    X_dir = [x_dir, x_zero;
             x_zero, x_dir;
             x_dir,x_ones;
             x_ones,x_dir]; 
    %X_neu = rand_uniform(n_neu,{[0,1]})';
    X_neu = zeros(0,2);
    N_neu = model.normals(X_neu);
   case 2 % choose points from grid for unit square and l-shape
    if model_choice == 1
      p.xnumintervals = 2*2^refinements;
      p.ynumintervals = 2*2^refinements;
      p.xrange = [0,1];
      p.yrange = [0,1];
    else
      p.xnumintervals = 5*2^refinements;
      p.ynumintervals = 5*2^refinements;
      p.xrange = [-1,1];
      p.yrange = [-1,1];
    end;
    grid = triagrid(p);
    Xall = [grid.X(:), grid.Y(:)];
    if model_choice == 2
      keep_points = ones(size(Xall,1),1);
      Xremove = find( (Xall(:,1)>1e-10) & (Xall(:,2)>1e-10));
      Xremove0 = find( sum(Xall.^2,2)<0.26^2);% remove origin ball
                                               % due
                                              % to f singularity
      Xremove = [Xremove;Xremove0];
      keep_points(Xremove) = 0;
      Xall = Xall(find(keep_points),:);
    end;
    bt = model.boundary_type(Xall);
%    keyboard;
    i_int = find(bt == 0);
    i_dir = find(bt == -1);
    i_neu = find(bt == -2);
    
    n_int = length(i_int);
    n_dir = length(i_dir);
    n_neu = length(i_neu);
    X_int = Xall(i_int,:);
    X_dir = Xall(i_dir,:);
    X_neu = Xall(i_neu,:);
    N_neu = model.normals(X_neu);
    
   otherwise('colpoint_choice unknown')
  end;
  
  X = [X_int;X_dir;X_neu];
  n = size(X,1);
  if n~=(n_int+n_dir+n_neu) 
    error('error in system dimension!');
  end;
  
  % definition of elementary functions
  % evaluation of all basis functions in single point:
  % (phi_i(x),...,phi_n(x))
  phi_vec = @(x) ...
            kernel(sqrt(sum((repmat(x(:)',n,1)-X).^2,2)));
  % auxiliary function grad_phi_vec_aux defined at end of file
  grad_phi_vec = @(x) ...
      grad_phi_vec_aux(x,kernel_derivative,n,X);
    % d_11 phi + d_22 phi = kernel''(|x-xi|)
  %                       +
  %                       kernel'(|x-xi|)*
  %                         (2-(x1+x2-xi1-xi2)|x-xi|)/|x-xi|^2
  %div_grad_phi_vec = @(x) ...
  %    kernel_2nd_derivative(sqrt(sum((repmat(x(:)',n,1)-X).^2,2))) ...
  %    + kernel_derivative(sqrt(sum((repmat(x(:)',n,1)-X).^2,2))) .* ...
  %    ((-sum(X,2)+x(1)+x(2)).*sqrt(sum((repmat(x(:)',n,1)-X).^2,2)) ...
  %     + 2) ...
  %    .* sqrt(sum((repmat(x(:)',n,1)-X).^2,2)).^(-2);
  % evaluation of all (a * grad phi_j) in single point numerical derivative
  %a_grad_phi_vec = @(x) ...
  %    model.diffusivity(x(:)') * (num_gradient(phi_vec,x(:)'));
  % evaluation of all (a * grad phi_j) in single point
  a_grad_phi_vec = @(x) ...
      model.diffusivity(x(:)') * grad_phi_vec(x(:)');
  a_grad_phi_n_vec = @(x,n) ...
      a_grad_phi_vec(x) * n(:);
  
  % evaluation of all (a * grad phi_j * n) in single point with
  % normal n
  %grad_phi_vec = @(x) ...
  %       kernel(sum((repmat(x(:)',n,1)-X).^2,2));  
  
  %evaluation of all   - a div grad phi_i(x) in single point x
  %minus_a_div_grad_phi_vec = @(x) ...
  %    - model.diffusivity(x(:)') * ...
  %    div_grad_phi_vec(x);      
  %evaluation of all   - (grad a) (grad phi_i(x)) in single point x
  %minus_grad_a_grad_phi_vec = @(x) ...
  %    - sum(repmat(model.diffusivity_gradient(x(:)'),size(X,1),1)...
  %    .*grad_phi_vec(x),2);
  %evaluation of all (- div (a grad phi_i(x))) in single point x
  %minus_div_a_grad_phi_vec = @(x) ...
  %    minus_grad_a_grad_phi_vec(x)+ minus_a_div_grad_phi_vec(x);
  %evaluation of all (- div (a grad phi_i(x))) in single point x
  %numerical derivation:
  ds = (1e-3)*gamma;
  minus_div_a_grad_phi_vec = @(x) ...
      - ...
      ((a_grad_phi_vec(x+[ds/2,0])-...
       a_grad_phi_vec(x-[ds/2,0]))/ds) * [1;0] ...
      - ...
      ((a_grad_phi_vec(x+[0,ds/2])-...
       a_grad_phi_vec(x-[0,ds/2]))/ds) * [0;1];
  
  % assemble system, sorry one loop, can be removed later
  A_int = zeros(n_int,n);
  for i = 1:n_int;
    A_int(i,:) = minus_div_a_grad_phi_vec(X_int(i,:))';
  end;
  A_dir = zeros(n_dir,n);
  for i = 1:n_dir;
    A_dir(i,:) = phi_vec(X_dir(i,:))';
  end;
  A_neu = zeros(n_neu,n);
  for i = 1:n_neu;
    A_neu(i,:) = a_grad_phi_n_vec(X_neu(i,:),N_neu(i,:))';
  end;
%  keyboard;
  
  A = [A_int; A_dir; A_neu];
  
  b_int = model.source(X_int);
  b_dir = model.dirichlet_values(X_dir);
  b_neu = model.neumann_values(X_neu);
  b = [b_int;b_dir;b_neu];
%  keyboard;
  
  % solve system
  lambda = A\b;
%  lambda = bicgstab(A,b);
  
  % discrete function
  %u(x) = sum_i lambda_i phi_i(x)
  u_h = @(x) lambda' * phi_vec(x); 
  
  % evaluation of discrete function in collocation points
  u = zeros(n,1);
  % sorry, one loop, can be removed later...
  for i = 1:n;
    u(i) = u_h(X(i,:));
  end;
  
  % error to exact solution
  u_exact = model.solution(X);
  max_err = max(abs(u-u_exact));

  % plots if required
  if ~no_plots    
    % right hand side
    figure, plot3(X(:,1),X(:,2),b,'*');
    title('right hand side of eqn system');
    % matrix
    figure, pcolor(A);
    colorbar;
    title('matrix of eqn system');
    % discrete solution
    bt = model.boundary_type(X);
    dir_ind = find(bt==-1);
    neu_ind = find(bt==-2);
    legends = {'int','dir','neu'};
    if isempty(neu_ind)
      legends = {'int','dir'};
    end;
    figure, plot3(X(:,1),X(:,2),u,'*');
    hold on, plot3(X(dir_ind,1),X(dir_ind,2),u(dir_ind),'ro');
    hold on, plot3(X(neu_ind,1),X(neu_ind,2),u(neu_ind),'go');
    title('discrete solution')
    legend(legends);
    % exact solution
    figure, plot3(X(:,1),X(:,2),u_exact,'*');
    hold on, plot3(X(dir_ind,1),X(dir_ind,2),u_exact(dir_ind),'ro');
    hold on, plot3(X(neu_ind,1),X(neu_ind,2),u_exact(neu_ind),'go');
    title('exact solution');
    legend(legends);
    % error
    figure, plot3(X(:,1),X(:,2),u-u_exact,'*');
    title('error');
    disp(['L-infty error: ',num2str(max_err)])
  end;  

  % return results
  res.u_h = u_h;
  res.err = max_err;
  res.cond = cond(A);

 case 2 % working example for poisson on unit square
  params = [];
  params.model_choice = 1; % poisson
%  params.gamma = 0.0063; % for ref = 0;
%  params.refinements = 0;
%  params.gamma = 0.063; % for ref = 1;
%  params.refinements = 1;
  params.gamma = 0.794; % for ref = 2;
  params.refinements = 2;
  params.all_dirichlet_boundary = 0;
%  params.all_dirichlet_boundary = 1;
  % call collocation
  res = rbf_elliptic(1,params);
%  res = rbf_elliptic(4,params);
  
 case 3 % working example for difficult l-shaped domain with ellipic model
  params = [];
  params.model_choice = 2; % elliptic l-shaped
%  params.gamma = 0.001; % optimal for refinements 0
  params.gamma = 31.62; % optimal for refinements 1
%  params.gamma = 80; % optimal for refinements 2
  params.refinements = 1; 
  params.all_dirichlet_boundary = 1;
%  params.all_dirichlet_boundary = 0;
  % call collocation
  res = rbf_elliptic(1,params);
%  res = rbf_elliptic(4,params);
    
 case 4 % determine errors and matrix conditions over a range of gammas
  % range of gamma between 0.001 and 10, log-equidistant
  gamma_log = log(10)*(-3:(1/10):3);
  gammas = exp(gamma_log);
  errs = zeros(size(gammas));
  conds = zeros(size(gammas));
  for gamma_index = 1:length(gammas);
    params.gamma = gammas(gamma_index);
    params.no_plots = 1;
    params.colpoints_choice = 2;
    % solve for given parameter
    res = rbf_elliptic(1,params);
    errs(gamma_index) = res.err;
    conds(gamma_index) = res.cond;
  end;
  figure,plot(gammas,errs);
  set(gca,'Xscale','log');
  set(gca,'Yscale','log');
  xlabel('gamma');
  ylabel('l_\infty-error');
  title('error over gamma');
  figure,plot(gammas,conds);
  set(gca,'Xscale','log');
  set(gca,'Yscale','log');
  xlabel('gamma');
  ylabel('cond(A)');
  title('condition number over gamma');
  res.gammas = gammas;
  res.errs = errs;
  res.conds = conds;  
  [e,i] = min(errs);
  res.gamma_opt = gammas(i(1));
  disp(['optimal gamma:',num2str(res.gamma_opt),...
        ', linfty-error : ',num2str(e)]);
 otherwise 
  error('step unknown')
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% auxiliary functions %%%%%%%%%%%%%%%%%%%%%%%%%

% small auxiliary function used for computing grad_phi_vec
function res = grad_phi_vec_aux(x,kernel_derivative,n,X)
xdiffs = repmat(x(:)',n,1)-X;
xdiffnorms = sqrt(sum(xdiffs.^2,2));
res = repmat(...
    kernel_derivative(xdiffnorms).* (xdiffnorms.^(-1)), ...
    1,2)...
      .*xdiffs;
% by factor 0/0 (division by zero) a nan can appear, det to 0:
i = find(isnan(res));
res(i) = 0;