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.
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
15 % step 2 % working example
for poisson
17 % step 3 % working example
for difficult ellipic model
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
26 if (nargin < 1) | (isempty(step))
27 step = 2; %
default: show working example of poisson model
36 % check
if plots are required
37 if isfield(params,
'no_plots')
38 no_plots = params.no_plots;
39 else % default: plotting on
43 % set gamma of rbf-kernel
44 if isfield(params,'gamma')
51 % set gamma of rbf-kernel
52 if isfield(params,'colpoints_choice')
53 colpoints_choice = params.colpoints_choice;
59 if isfield(params,'model_choice')
60 model_choice = params.model_choice;
62 model_choice = 1; % poisson
63 %model_choice = 2; % elliptic
66 % set refinement steps
67 if isfield(params,'refinements')
68 refinements = params.refinements;
74 model = praktikum_poisson_model(params);
75 elseif model_choice == 2
76 model = praktikum_ellipt_model(params);
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)
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];
94 % collocation centers for unit square
95 switch colpoints_choice
96 case 1 % random points for unit square
98 error('random points not implemented for l-shape domain')
100 rand('seed',101); % initialize random number generator
101 n_int = 15 * 2^refinements;
102 n_dir = 8 * 2^refinements; % dividable by 4
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);
109 X_dir = [x_dir, x_zero;
113 %X_neu = rand_uniform(n_neu,{[0,1]})
';
115 N_neu = model.normals(X_neu);
116 case 2 % choose points from grid for unit square and l-shape
118 p.xnumintervals = 2*2^refinements;
119 p.ynumintervals = 2*2^refinements;
123 p.xnumintervals = 5*2^refinements;
124 p.ynumintervals = 5*2^refinements;
129 Xall = [grid.X(:), grid.Y(:)];
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
136 Xremove = [Xremove;Xremove0];
137 keep_points(Xremove) = 0;
138 Xall = Xall(find(keep_points),:);
140 bt = model.boundary_type(Xall);
142 i_int = find(bt == 0);
143 i_dir = find(bt == -1);
144 i_neu = find(bt == -2);
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);
154 otherwise('colpoint_choice unknown
')
157 X = [X_int;X_dir;X_neu];
159 if n~=(n_int+n_dir+n_neu)
160 error('error in system dimension!
');
163 % definition of elementary functions
164 % evaluation of all basis functions in single point:
165 % (phi_i(x),...,phi_n(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|)
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)) ...
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(:);
190 % evaluation of all (a * grad phi_j * n) in single point with
192 %grad_phi_vec = @(x) ...
193 % kernel(sum((repmat(x(:)',n,1)-X).^2,2));
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:
209 minus_div_a_grad_phi_vec = @(x) ...
211 ((a_grad_phi_vec(x+[ds/2,0])-...
212 a_grad_phi_vec(x-[ds/2,0]))/ds) * [1;0] ...
214 ((a_grad_phi_vec(x+[0,ds/2])-...
215 a_grad_phi_vec(x-[0,ds/2]))/ds) * [0;1];
217 % assemble system, sorry one loop, can be removed later
218 A_int = zeros(n_int,n);
220 A_int(i,:) = minus_div_a_grad_phi_vec(X_int(i,:))';
222 A_dir = zeros(n_dir,n);
224 A_dir(i,:) = phi_vec(X_dir(i,:))';
226 A_neu = zeros(n_neu,n);
228 A_neu(i,:) = a_grad_phi_n_vec(X_neu(i,:),N_neu(i,:))';
232 A = [A_int; A_dir; A_neu];
234 b_int = model.source(X_int);
236 b_neu = model.neumann_values(X_neu);
237 b = [b_int;b_dir;b_neu];
242 % lambda = bicgstab(A,b);
245 %u(x) = sum_i lambda_i phi_i(x)
246 u_h = @(x) lambda' * phi_vec(x);
248 % evaluation of discrete function in collocation points
250 % sorry, one loop, can be removed later...
255 % error to exact solution
256 u_exact = model.solution(X);
257 max_err = max(abs(u-u_exact));
262 figure, plot3(X(:,1),X(:,2),b,'*');
263 title('right hand side of eqn system');
267 title('matrix of eqn system');
269 bt = model.boundary_type(X);
270 dir_ind = find(bt==-1);
271 neu_ind = find(bt==-2);
272 legends = {
'int',
'dir',
'neu'};
274 legends = {
'int',
'dir'};
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')
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');
288 figure, plot3(X(:,1),X(:,2),u-u_exact,'*');
290 disp(['L-infty error: ',num2str(max_err)])
298 case 2 % working example for poisson on unit square
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;
311 % res = rbf_elliptic(4,params);
313 case 3 % working example for difficult l-shaped domain with ellipic model
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;
323 res = rbf_elliptic(1,params);
324 % res = rbf_elliptic(4,params);
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);
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;
341 figure,plot(gammas,errs);
342 set(gca,'Xscale','log');
343 set(gca,'Yscale','log');
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');
352 title('condition number over gamma');
357 res.gamma_opt = gammas(i(1));
358 disp(['optimal gamma:',num2str(res.gamma_opt),...
359 ', linfty-error : ',num2str(e)]);
361 error('step unknown')
364 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% auxiliary functions %%%%%%%%%%%%%%%%%%%%%%%%%
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));
371 kernel_derivative(xdiffnorms).* (xdiffnorms.^(-1)), ...
374 % by factor 0/0 (division by zero) a nan can appear, det to 0:
375 i = find(isnan(res));
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...