rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
rbf_elliptic.m
Go to the documentation of this file.
1 function res = rbf_elliptic(step,params);
2 %function res = rbf_elliptic(step,params);
3 % function for meshless collocation of fem-problem. step can be
4 % used to choose several cases. step=1 does not need to be called
5 % directly but is used by the other steps. So call step=2, 3 and
6 % change options there.
7 %
8 % step 1: single solve. params.gamma can be set,
9 % params.colpoints_choice can be set to 1(random) or 2(structured).
10 % params.model_choice = 1 (poisson) or 2 (singular elliptic)
11 % If params.no_plots is set, then no plots
12 % and no screen output are performed. Output fields are
13 % res.u_h and the l2-error res.err and condition number res.cond
14 %
15 % step 2 % working example for poisson
16 %
17 % step 3 % working example for difficult ellipic model
18 %
19 % step 4 % determine errors and matrix conditions over a range of gammas
20 % and plot of error by repeated calls of step 1. Finds optimal gamma
21 
22 % B. Haasdonk
23 
24 res = [];
25 
26 if (nargin < 1) | (isempty(step))
27  step = 2; % default: show working example of poisson model
28 end;
29 
30 if (nargin < 2)
31  params = [];
32 end;
33 
34 switch step
35  case 1
36  % check if plots are required
37  if isfield(params,'no_plots')
38  no_plots = params.no_plots;
39  else % default: plotting on
40  no_plots = 0;
41  end;
42 
43  % set gamma of rbf-kernel
44  if isfield(params,'gamma')
45  gamma = params.gamma;
46  else % default
47  gamma = 50;
48 % gamma = 0.01;
49  end;
50 
51  % set gamma of rbf-kernel
52  if isfield(params,'colpoints_choice')
53  colpoints_choice = params.colpoints_choice;
54  else
55  colpoints_choice = 2;
56  end;
57 
58  % set model type
59  if isfield(params,'model_choice')
60  model_choice = params.model_choice;
61  else
62  model_choice = 1; % poisson
63  %model_choice = 2; % elliptic
64  end;
65 
66  % set refinement steps
67  if isfield(params,'refinements')
68  refinements = params.refinements;
69  else
70  refinements = 2;
71  end;
72 
73  if model_choice == 1
74  model = praktikum_poisson_model(params);
75  elseif model_choice == 2
76  model = praktikum_ellipt_model(params);
77  end;
78 
79  % scalar argument and scalar valued rbf kernel and derivatives
80  % function can be called in many points s simultaneous.
81  % s is scalar (vectorial if many points)
82  kernel = @(s) exp(-gamma * s.^2);
83  kernel_derivative = @(s) -2*gamma * s.*exp(-gamma * s.^2);
84  %kernel_2nd_derivative = @(s) -2*gamma*exp(-gamma * s.^2).* ...
85  % (- 2 *gamma * s + 1);
86  % numerical gradient determination by central difference of a function f
87  % x is 1x2-vec (or matrix n x 2 if many points)
88  %ds = (1e-3)*gamma;
89 % num_gradient = @(f,x) [(f([x(:,1),x(:,2)+ds/2]) ...
90 % -f([x(:,1),x(:,2)-ds/2]))/ds, ...
91 % (f([x(:,1)+ds/2,x(:,2)]) ...
92 % -f([x(:,1)-ds/2,x(:,2)]))/ds];
93 
94 % collocation centers for unit square
95 switch colpoints_choice
96  case 1 % random points for unit square
97  if model_choice ==2
98  error('random points not implemented for l-shape domain')
99  end;
100  rand('seed',101); % initialize random number generator
101  n_int = 15 * 2^refinements;
102  n_dir = 8 * 2^refinements; % dividable by 4
103  n_neu = 0;
104  X_int = rand(n_int,2);
105  x_dir = rand(n_dir/4,1);
106  x_zero = zeros(size(x_dir,1),1);
107  x_ones = ones(size(x_dir,1),1);
108 
109  X_dir = [x_dir, x_zero;
110  x_zero, x_dir;
111  x_dir,x_ones;
112  x_ones,x_dir];
113  %X_neu = rand_uniform(n_neu,{[0,1]})';
114  X_neu = zeros(0,2);
115  N_neu = model.normals(X_neu);
116  case 2 % choose points from grid for unit square and l-shape
117  if model_choice == 1
118  p.xnumintervals = 2*2^refinements;
119  p.ynumintervals = 2*2^refinements;
120  p.xrange = [0,1];
121  p.yrange = [0,1];
122  else
123  p.xnumintervals = 5*2^refinements;
124  p.ynumintervals = 5*2^refinements;
125  p.xrange = [-1,1];
126  p.yrange = [-1,1];
127  end;
128  grid = triagrid(p);
129  Xall = [grid.X(:), grid.Y(:)];
130  if model_choice == 2
131  keep_points = ones(size(Xall,1),1);
132  Xremove = find( (Xall(:,1)>1e-10) & (Xall(:,2)>1e-10));
133  Xremove0 = find( sum(Xall.^2,2)<0.26^2);% remove origin ball
134  % due
135  % to f singularity
136  Xremove = [Xremove;Xremove0];
137  keep_points(Xremove) = 0;
138  Xall = Xall(find(keep_points),:);
139  end;
140  bt = model.boundary_type(Xall);
141 % keyboard;
142  i_int = find(bt == 0);
143  i_dir = find(bt == -1);
144  i_neu = find(bt == -2);
145 
146  n_int = length(i_int);
147  n_dir = length(i_dir);
148  n_neu = length(i_neu);
149  X_int = Xall(i_int,:);
150  X_dir = Xall(i_dir,:);
151  X_neu = Xall(i_neu,:);
152  N_neu = model.normals(X_neu);
153 
154  otherwise('colpoint_choice unknown')
155  end;
156 
157  X = [X_int;X_dir;X_neu];
158  n = size(X,1);
159  if n~=(n_int+n_dir+n_neu)
160  error('error in system dimension!');
161  end;
162 
163  % definition of elementary functions
164  % evaluation of all basis functions in single point:
165  % (phi_i(x),...,phi_n(x))
166  phi_vec = @(x) ...
167  kernel(sqrt(sum((repmat(x(:)',n,1)-X).^2,2)));
168  % auxiliary function grad_phi_vec_aux defined at end of file
169  grad_phi_vec = @(x) ...
170  grad_phi_vec_aux(x,kernel_derivative,n,X);
171  % d_11 phi + d_22 phi = kernel''(|x-xi|)
172  % +
173  % kernel'(|x-xi|)*
174  % (2-(x1+x2-xi1-xi2)|x-xi|)/|x-xi|^2
175  %div_grad_phi_vec = @(x) ...
176  % kernel_2nd_derivative(sqrt(sum((repmat(x(:)',n,1)-X).^2,2))) ...
177  % + kernel_derivative(sqrt(sum((repmat(x(:)',n,1)-X).^2,2))) .* ...
178  % ((-sum(X,2)+x(1)+x(2)).*sqrt(sum((repmat(x(:)',n,1)-X).^2,2)) ...
179  % + 2) ...
180  % .* sqrt(sum((repmat(x(:)',n,1)-X).^2,2)).^(-2);
181  % evaluation of all (a * grad phi_j) in single point numerical derivative
182  %a_grad_phi_vec = @(x) ...
183  % model.diffusivity(x(:)') * (num_gradient(phi_vec,x(:)'));
184  % evaluation of all (a * grad phi_j) in single point
185  a_grad_phi_vec = @(x) ...
186  model.diffusivity(x(:)') * grad_phi_vec(x(:)');
187  a_grad_phi_n_vec = @(x,n) ...
188  a_grad_phi_vec(x) * n(:);
189 
190  % evaluation of all (a * grad phi_j * n) in single point with
191  % normal n
192  %grad_phi_vec = @(x) ...
193  % kernel(sum((repmat(x(:)',n,1)-X).^2,2));
194 
195  %evaluation of all - a div grad phi_i(x) in single point x
196  %minus_a_div_grad_phi_vec = @(x) ...
197  % - model.diffusivity(x(:)') * ...
198  % div_grad_phi_vec(x);
199  %evaluation of all - (grad a) (grad phi_i(x)) in single point x
200  %minus_grad_a_grad_phi_vec = @(x) ...
201  % - sum(repmat(model.diffusivity_gradient(x(:)'),size(X,1),1)...
202  % .*grad_phi_vec(x),2);
203  %evaluation of all (- div (a grad phi_i(x))) in single point x
204  %minus_div_a_grad_phi_vec = @(x) ...
205  % minus_grad_a_grad_phi_vec(x)+ minus_a_div_grad_phi_vec(x);
206  %evaluation of all (- div (a grad phi_i(x))) in single point x
207  %numerical derivation:
208  ds = (1e-3)*gamma;
209  minus_div_a_grad_phi_vec = @(x) ...
210  - ...
211  ((a_grad_phi_vec(x+[ds/2,0])-...
212  a_grad_phi_vec(x-[ds/2,0]))/ds) * [1;0] ...
213  - ...
214  ((a_grad_phi_vec(x+[0,ds/2])-...
215  a_grad_phi_vec(x-[0,ds/2]))/ds) * [0;1];
216 
217  % assemble system, sorry one loop, can be removed later
218  A_int = zeros(n_int,n);
219  for i = 1:n_int;
220  A_int(i,:) = minus_div_a_grad_phi_vec(X_int(i,:))';
221  end;
222  A_dir = zeros(n_dir,n);
223  for i = 1:n_dir;
224  A_dir(i,:) = phi_vec(X_dir(i,:))';
225  end;
226  A_neu = zeros(n_neu,n);
227  for i = 1:n_neu;
228  A_neu(i,:) = a_grad_phi_n_vec(X_neu(i,:),N_neu(i,:))';
229  end;
230 % keyboard;
231 
232  A = [A_int; A_dir; A_neu];
233 
234  b_int = model.source(X_int);
235  b_dir = model.dirichlet_values(X_dir);
236  b_neu = model.neumann_values(X_neu);
237  b = [b_int;b_dir;b_neu];
238 % keyboard;
239 
240  % solve system
241  lambda = A\b;
242 % lambda = bicgstab(A,b);
243 
244  % discrete function
245  %u(x) = sum_i lambda_i phi_i(x)
246  u_h = @(x) lambda' * phi_vec(x);
247 
248  % evaluation of discrete function in collocation points
249  u = zeros(n,1);
250  % sorry, one loop, can be removed later...
251  for i = 1:n;
252  u(i) = u_h(X(i,:));
253  end;
254 
255  % error to exact solution
256  u_exact = model.solution(X);
257  max_err = max(abs(u-u_exact));
258 
259  % plots if required
260  if ~no_plots
261  % right hand side
262  figure, plot3(X(:,1),X(:,2),b,'*');
263  title('right hand side of eqn system');
264  % matrix
265  figure, pcolor(A);
266  colorbar;
267  title('matrix of eqn system');
268  % discrete solution
269  bt = model.boundary_type(X);
270  dir_ind = find(bt==-1);
271  neu_ind = find(bt==-2);
272  legends = {'int','dir','neu'};
273  if isempty(neu_ind)
274  legends = {'int','dir'};
275  end;
276  figure, plot3(X(:,1),X(:,2),u,'*');
277  hold on, plot3(X(dir_ind,1),X(dir_ind,2),u(dir_ind),'ro');
278  hold on, plot3(X(neu_ind,1),X(neu_ind,2),u(neu_ind),'go');
279  title('discrete solution')
280  legend(legends);
281  % exact solution
282  figure, plot3(X(:,1),X(:,2),u_exact,'*');
283  hold on, plot3(X(dir_ind,1),X(dir_ind,2),u_exact(dir_ind),'ro');
284  hold on, plot3(X(neu_ind,1),X(neu_ind,2),u_exact(neu_ind),'go');
285  title('exact solution');
286  legend(legends);
287  % error
288  figure, plot3(X(:,1),X(:,2),u-u_exact,'*');
289  title('error');
290  disp(['L-infty error: ',num2str(max_err)])
291  end;
292 
293  % return results
294  res.u_h = u_h;
295  res.err = max_err;
296  res.cond = cond(A);
297 
298  case 2 % working example for poisson on unit square
299  params = [];
300  params.model_choice = 1; % poisson
301 % params.gamma = 0.0063; % for ref = 0;
302 % params.refinements = 0;
303 % params.gamma = 0.063; % for ref = 1;
304 % params.refinements = 1;
305  params.gamma = 0.794; % for ref = 2;
306  params.refinements = 2;
307  params.all_dirichlet_boundary = 0;
308 % params.all_dirichlet_boundary = 1;
309  % call collocation
310  res = rbf_elliptic(1,params);
311 % res = rbf_elliptic(4,params);
312 
313  case 3 % working example for difficult l-shaped domain with ellipic model
314  params = [];
315  params.model_choice = 2; % elliptic l-shaped
316 % params.gamma = 0.001; % optimal for refinements 0
317  params.gamma = 31.62; % optimal for refinements 1
318 % params.gamma = 80; % optimal for refinements 2
319  params.refinements = 1;
320  params.all_dirichlet_boundary = 1;
321 % params.all_dirichlet_boundary = 0;
322  % call collocation
323  res = rbf_elliptic(1,params);
324 % res = rbf_elliptic(4,params);
325 
326  case 4 % determine errors and matrix conditions over a range of gammas
327  % range of gamma between 0.001 and 10, log-equidistant
328  gamma_log = log(10)*(-3:(1/10):3);
329  gammas = exp(gamma_log);
330  errs = zeros(size(gammas));
331  conds = zeros(size(gammas));
332  for gamma_index = 1:length(gammas);
333  params.gamma = gammas(gamma_index);
334  params.no_plots = 1;
335  params.colpoints_choice = 2;
336  % solve for given parameter
337  res = rbf_elliptic(1,params);
338  errs(gamma_index) = res.err;
339  conds(gamma_index) = res.cond;
340  end;
341  figure,plot(gammas,errs);
342  set(gca,'Xscale','log');
343  set(gca,'Yscale','log');
344  xlabel('gamma');
345  ylabel('l_\infty-error');
346  title('error over gamma');
347  figure,plot(gammas,conds);
348  set(gca,'Xscale','log');
349  set(gca,'Yscale','log');
350  xlabel('gamma');
351  ylabel('cond(A)');
352  title('condition number over gamma');
353  res.gammas = gammas;
354  res.errs = errs;
355  res.conds = conds;
356  [e,i] = min(errs);
357  res.gamma_opt = gammas(i(1));
358  disp(['optimal gamma:',num2str(res.gamma_opt),...
359  ', linfty-error : ',num2str(e)]);
360  otherwise
361  error('step unknown')
362 end;
363 
364 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% auxiliary functions %%%%%%%%%%%%%%%%%%%%%%%%%
365 
366 % small auxiliary function used for computing grad_phi_vec
367 function res = grad_phi_vec_aux(x,kernel_derivative,n,X)
368 xdiffs = repmat(x(:)',n,1)-X;
369 xdiffnorms = sqrt(sum(xdiffs.^2,2));
370 res = repmat(...
371  kernel_derivative(xdiffnorms).* (xdiffnorms.^(-1)), ...
372  1,2)...
373  .*xdiffs;
374 % by factor 0/0 (division by zero) a nan can appear, det to 0:
375 i = find(isnan(res));
376 res(i) = 0;
function Udirichlet = dirichlet_values(model, X, Y)
UDIRICHLET = DIRICHLET_VALUES([X],[Y], MODEL) Examples dirichlet_values([0,1,2],[1,1,1],struct(name_dirichlet_values, homogeneous, ... c_dir, 1)) dirichlet_values([0:0.1:1],[0],struct(name_dirichlet_values, xstripes, ... c_dir, [0 1 2], ... dir_borders, [0.3 0.6])) dirichlet_values([0:0.1:1],[0],struct(name_dirichlet_values, box, ... c_dir, 1, ... dir_box_xrange, [0.3 0.6], ... dir_box_yrange, [-0.1 0.1]))
function res = rbf_elliptic(step, params)
function for meshless collocation of fem-problem. step can be used to choose several cases...
Definition: rbf_elliptic.m:17