Actual source code: laplacian.m
petsc-3.12.2 2019-11-22
1: function [lambda, V, A] = laplacian(varargin)
3: % LAPLACIAN Sparse Negative Laplacian in 1D, 2D, or 3D
4: %
5: % [~,~,A]=LAPLACIAN(N) generates a sparse negative 3D Laplacian matrix
6: % with Dirichlet boundary conditions, from a rectangular cuboid regular
7: % grid with j x k x l interior grid points if N = [j k l], using the
8: % standard 7-point finite-difference scheme, The grid size is always
9: % one in all directions.
10: %
11: % [~,~,A]=LAPLACIAN(N,B) specifies boundary conditions with a cell array
12: % B. For example, B = {'DD' 'DN' 'P'} will Dirichlet boundary conditions
13: % ('DD') in the x-direction, Dirichlet-Neumann conditions ('DN') in the
14: % y-direction and period conditions ('P') in the z-direction. Possible
15: % values for the elements of B are 'DD', 'DN', 'ND', 'NN' and 'P'.
16: %
17: % LAMBDA = LAPLACIAN(N,B,M) or LAPLACIAN(N,M) outputs the m smallest
18: % eigenvalues of the matrix, computed by an exact known formula, see
19: % http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_of_the_second_derivative
20: % It will produce a warning if the mth eigenvalue is equal to the
21: % (m+1)th eigenvalue. If m is absebt or zero, lambda will be empty.
22: %
23: % [LAMBDA,V] = LAPLACIAN(N,B,M) also outputs orthonormal eigenvectors
24: % associated with the corresponding m smallest eigenvalues.
25: %
26: % [LAMBDA,V,A] = LAPLACIAN(N,B,M) produces a 2D or 1D negative
27: % Laplacian matrix if the length of N and B are 2 or 1 respectively.
28: % It uses the standard 5-point scheme for 2D, and 3-point scheme for 1D.
29: %
30: % % Examples:
31: % [lambda,V,A] = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
32: % % Everything for 3D negative Laplacian with mixed boundary conditions.
33: % laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
34: % % or
35: % lambda = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
36: % % computes the eigenvalues only
37: %
38: % [~,V,~] = laplacian([200 200],{'DD' 'DN'},30);
39: % % Eigenvectors of 2D negative Laplacian with mixed boundary conditions.
40: %
41: % [~,~,A] = laplacian(200,{'DN'},30);
42: % % 1D negative Laplacian matrix A with mixed boundary conditions.
43: %
44: % % Example to test if outputs correct eigenvalues and vectors:
45: % [lambda,V,A] = laplacian([13,10,6],{'DD' 'DN' 'P'},30);
46: % [Veig D] = eig(full(A)); lambdaeig = diag(D(1:30,1:30));
47: % max(abs(lambda-lambdaeig)) %checking eigenvalues
48: % subspace(V,Veig(:,1:30)) %checking the invariant subspace
49: % subspace(V(:,1),Veig(:,1)) %checking selected eigenvectors
50: % subspace(V(:,29:30),Veig(:,29:30)) %a multiple eigenvalue
51: %
52: % % Example showing equivalence between laplacian.m and built-in MATLAB
53: % % DELSQ for the 2D case. The output of the last command shall be 0.
54: % A1 = delsq(numgrid('S',32)); % input 'S' specifies square grid.
55: % [~,~,A2] = laplacian([30,30]);
56: % norm(A1-A2,inf)
57: %
58: % Class support for inputs:
59: % N - row vector float double
60: % B - cell array
61: % M - scalar float double
62: %
63: % Class support for outputs:
64: % lambda and V - full float double, A - sparse float double.
65: %
66: % Note: the actual numerical entries of A fit int8 format, but only
67: % double data class is currently (2010) supported for sparse matrices.
68: %
69: % This program is designed to efficiently compute eigenvalues,
70: % eigenvectors, and the sparse matrix of the (1-3)D negative Laplacian
71: % on a rectangular grid for Dirichlet, Neumann, and Periodic boundary
72: % conditions using tensor sums of 1D Laplacians. For more information on
73: % tensor products, see
74: % http://en.wikipedia.org/wiki/Kronecker_sum_of_discrete_Laplacians
75: % For 2D case in MATLAB, see
76: % http://www.mathworks.com/access/helpdesk/help/techdoc/ref/kron.html.
77: %
78: % This code is a part of the BLOPEX package:
79: % http://en.wikipedia.org/wiki/BLOPEX or directly
80: % http://code.google.com/p/blopex/
82: % Revision 1.1 changes: rearranged the output variables, always compute
83: % the eigenvalues, compute eigenvectors and/or the matrix on demand only.
85: % License: BSD
86: % Copyright 2010-2011 Bryan C. Smith, Andrew V. Knyazev
87: % $Revision: 1.1 $ $Date: 1-Aug-2011
88: % Tested in MATLAB 7.11.0 (R2010b) and Octave 3.4.0.
90: tic
92: % Input/Output handling.
93: if nargin > 3
94: error('BLOPEX:laplacian:TooManyInputs',...
95: '%s','Too many input arguments.');
96: elseif nargin == 0
97: error('BLOPEX:laplacian:NoInputArguments',...
98: '%s','Must have at least one input argument.');
99: end
101: if nargout > 3
102: error('BLOPEX:laplacian:TooManyOutputs',...
103: '%s','Maximum number of outputs is 3.');
104: end
106: u = varargin{1};
107: dim2 = size(u);
108: if dim2(1) ~= 1
109: error('BLOPEX:laplacian:WrongVectorOfGridPoints',...
110: '%s','Number of grid points must be in a row vector.')
111: end
112: if dim2(2) > 3
113: error('BLOPEX:laplacian:WrongNumberOfGridPoints',...
114: '%s','Number of grid points must be a 1, 2, or 3')
115: end
116: dim=dim2(2); clear dim2;
118: uint = round(u);
119: if max(uint~=u)
120: warning('BLOPEX:laplacian:NonIntegerGridSize',...
121: '%s','Grid sizes must be integers. Rounding...');
122: u = uint; clear uint
123: end
124: if max(u<=0 )
125: error('BLOPEX:laplacian:NonIntegerGridSize',...
126: '%s','Grid sizes must be positive.');
127: end
129: if nargin == 3
130: m = varargin{3};
131: B = varargin{2};
132: elseif nargin == 2
133: f = varargin{2};
134: a = whos('regep','f');
135: if sum(a.class(1:4)=='cell') == 4
136: B = f;
137: m = 0;
138: elseif sum(a.class(1:4)=='doub') == 4
139: if dim == 1
140: B = {'DD'};
141: elseif dim == 2
142: B = {'DD' 'DD'};
143: else
144: B = {'DD' 'DD' 'DD'};
145: end
146: m = f;
147: else
148: error('BLOPEX:laplacian:InvalidClass',...
149: '%s','Second input must be either class double or a cell array.');
150: end
151: else
152: if dim == 1
153: B = {'DD'};
154: elseif dim == 2
155: B = {'DD' 'DD'};
156: else
157: B = {'DD' 'DD' 'DD'};
158: end
159: m = 0;
160: end
162: if max(size(m) - [1 1]) ~= 0
163: error('BLOPEX:laplacian:WrongNumberOfEigenvalues',...
164: '%s','The requested number of eigenvalues must be a scalar.');
165: end
167: maxeigs = prod(u);
168: mint = round(m);
169: if mint ~= m || mint > maxeigs
170: error('BLOPEX:laplacian:InvalidNumberOfEigs',...
171: '%s','Number of eigenvalues output must be a nonnegative ',...
172: 'integer no bigger than number of grid points.');
173: end
174: m = mint;
176: bdryerr = 0;
177: a = whos('regep','B');
178: if sum(a.class(1:4)=='cell') ~= 4 || sum(a.size == [1 dim]) ~= 2
179: bdryerr = 1;
180: else
181: BB = zeros(1, 2*dim);
182: for i = 1:dim
183: if (length(B{i}) == 1)
184: if B{i} == 'P'
185: BB(i) = 3;
186: BB(i + dim) = 3;
187: else
188: bdryerr = 1;
189: end
190: elseif (length(B{i}) == 2)
191: if B{i}(1) == 'D'
192: BB(i) = 1;
193: elseif B{i}(1) == 'N'
194: BB(i) = 2;
195: else
196: bdryerr = 1;
197: end
198: if B{i}(2) == 'D'
199: BB(i + dim) = 1;
200: elseif B{i}(2) == 'N'
201: BB(i + dim) = 2;
202: else
203: bdryerr = 1;
204: end
205: else
206: bdryerr = 1;
207: end
208: end
209: end
211: if bdryerr == 1
212: error('BLOPEX:laplacian:InvalidBdryConds',...
213: '%s','Boundary conditions must be a cell of length 3 for 3D, 2',...
214: ' for 2D, 1 for 1D, with values ''DD'', ''DN'', ''ND'', ''NN''',...
215: ', or ''P''.');
216: end
218: % Set the component matrices. SPDIAGS converts int8 into double anyway.
219: e1 = ones(u(1),1); %e1 = ones(u(1),1,'int8');
220: if dim > 1
221: e2 = ones(u(2),1);
222: end
223: if dim > 2
224: e3 = ones(u(3),1);
225: end
227: % Calculate m smallest exact eigenvalues.
228: if m > 0
229: if (BB(1) == 1) && (BB(1+dim) == 1)
230: a1 = pi/2/(u(1)+1);
231: N = (1:u(1))';
232: elseif (BB(1) == 2) && (BB(1+dim) == 2)
233: a1 = pi/2/u(1);
234: N = (0:(u(1)-1))';
235: elseif ((BB(1) == 1) && (BB(1+dim) == 2)) || ((BB(1) == 2)...
236: && (BB(1+dim) == 1))
237: a1 = pi/4/(u(1)+0.5);
238: N = 2*(1:u(1))' - 1;
239: else
240: a1 = pi/u(1);
241: N = floor((1:u(1))/2)';
242: end
243:
244: lambda1 = 4*sin(a1*N).^2;
245:
246: if dim > 1
247: if (BB(2) == 1) && (BB(2+dim) == 1)
248: a2 = pi/2/(u(2)+1);
249: N = (1:u(2))';
250: elseif (BB(2) == 2) && (BB(2+dim) == 2)
251: a2 = pi/2/u(2);
252: N = (0:(u(2)-1))';
253: elseif ((BB(2) == 1) && (BB(2+dim) == 2)) || ((BB(2) == 2)...
254: && (BB(2+dim) == 1))
255: a2 = pi/4/(u(2)+0.5);
256: N = 2*(1:u(2))' - 1;
257: else
258: a2 = pi/u(2);
259: N = floor((1:u(2))/2)';
260: end
261: lambda2 = 4*sin(a2*N).^2;
262: else
263: lambda2 = 0;
264: end
265:
266: if dim > 2
267: if (BB(3) == 1) && (BB(6) == 1)
268: a3 = pi/2/(u(3)+1);
269: N = (1:u(3))';
270: elseif (BB(3) == 2) && (BB(6) == 2)
271: a3 = pi/2/u(3);
272: N = (0:(u(3)-1))';
273: elseif ((BB(3) == 1) && (BB(6) == 2)) || ((BB(3) == 2)...
274: && (BB(6) == 1))
275: a3 = pi/4/(u(3)+0.5);
276: N = 2*(1:u(3))' - 1;
277: else
278: a3 = pi/u(3);
279: N = floor((1:u(3))/2)';
280: end
281: lambda3 = 4*sin(a3*N).^2;
282: else
283: lambda3 = 0;
284: end
285:
286: if dim == 1
287: lambda = lambda1;
288: elseif dim == 2
289: lambda = kron(e2,lambda1) + kron(lambda2, e1);
290: else
291: lambda = kron(e3,kron(e2,lambda1)) + kron(e3,kron(lambda2,e1))...
292: + kron(lambda3,kron(e2,e1));
293: end
294: [lambda, p] = sort(lambda);
295: if m < maxeigs - 0.1
296: w = lambda(m+1);
297: else
298: w = inf;
299: end
300: lambda = lambda(1:m);
301: p = p(1:m)';
302: else
303: lambda = [];
304: end
306: V = [];
307: if nargout > 1 && m > 0 % Calculate eigenvectors if specified in output.
308:
309: p1 = mod(p-1,u(1))+1;
310:
311: if (BB(1) == 1) && (BB(1+dim) == 1)
312: V1 = sin(kron((1:u(1))'*(pi/(u(1)+1)),p1))*(2/(u(1)+1))^0.5;
313: elseif (BB(1) == 2) && (BB(1+dim) == 2)
314: V1 = cos(kron((0.5:1:u(1)-0.5)'*(pi/u(1)),p1-1))*(2/u(1))^0.5;
315: V1(:,p1==1) = 1/u(1)^0.5;
316: elseif ((BB(1) == 1) && (BB(1+dim) == 2))
317: V1 = sin(kron((1:u(1))'*(pi/2/(u(1)+0.5)),2*p1 - 1))...
318: *(2/(u(1)+0.5))^0.5;
319: elseif ((BB(1) == 2) && (BB(1+dim) == 1))
320: V1 = cos(kron((0.5:1:u(1)-0.5)'*(pi/2/(u(1)+0.5)),2*p1 - 1))...
321: *(2/(u(1)+0.5))^0.5;
322: else
323: V1 = zeros(u(1),m);
324: a = (0.5:1:u(1)-0.5)';
325: V1(:,mod(p1,2)==1) = cos(a*(pi/u(1)*(p1(mod(p1,2)==1)-1)))...
326: *(2/u(1))^0.5;
327: pp = p1(mod(p1,2)==0);
328: if ~isempty(pp)
329: V1(:,mod(p1,2)==0) = sin(a*(pi/u(1)*p1(mod(p1,2)==0)))...
330: *(2/u(1))^0.5;
331: end
332: V1(:,p1==1) = 1/u(1)^0.5;
333: if mod(u(1),2) == 0
334: V1(:,p1==u(1)) = V1(:,p1==u(1))/2^0.5;
335: end
336: end
337:
338:
339: if dim > 1
340: p2 = mod(p-p1,u(1)*u(2));
341: p3 = (p - p2 - p1)/(u(1)*u(2)) + 1;
342: p2 = p2/u(1) + 1;
343: if (BB(2) == 1) && (BB(2+dim) == 1)
344: V2 = sin(kron((1:u(2))'*(pi/(u(2)+1)),p2))*(2/(u(2)+1))^0.5;
345: elseif (BB(2) == 2) && (BB(2+dim) == 2)
346: V2 = cos(kron((0.5:1:u(2)-0.5)'*(pi/u(2)),p2-1))*(2/u(2))^0.5;
347: V2(:,p2==1) = 1/u(2)^0.5;
348: elseif ((BB(2) == 1) && (BB(2+dim) == 2))
349: V2 = sin(kron((1:u(2))'*(pi/2/(u(2)+0.5)),2*p2 - 1))...
350: *(2/(u(2)+0.5))^0.5;
351: elseif ((BB(2) == 2) && (BB(2+dim) == 1))
352: V2 = cos(kron((0.5:1:u(2)-0.5)'*(pi/2/(u(2)+0.5)),2*p2 - 1))...
353: *(2/(u(2)+0.5))^0.5;
354: else
355: V2 = zeros(u(2),m);
356: a = (0.5:1:u(2)-0.5)';
357: V2(:,mod(p2,2)==1) = cos(a*(pi/u(2)*(p2(mod(p2,2)==1)-1)))...
358: *(2/u(2))^0.5;
359: pp = p2(mod(p2,2)==0);
360: if ~isempty(pp)
361: V2(:,mod(p2,2)==0) = sin(a*(pi/u(2)*p2(mod(p2,2)==0)))...
362: *(2/u(2))^0.5;
363: end
364: V2(:,p2==1) = 1/u(2)^0.5;
365: if mod(u(2),2) == 0
366: V2(:,p2==u(2)) = V2(:,p2==u(2))/2^0.5;
367: end
368: end
369: else
370: V2 = ones(1,m);
371: end
372:
373: if dim > 2
374: if (BB(3) == 1) && (BB(3+dim) == 1)
375: V3 = sin(kron((1:u(3))'*(pi/(u(3)+1)),p3))*(2/(u(3)+1))^0.5;
376: elseif (BB(3) == 2) && (BB(3+dim) == 2)
377: V3 = cos(kron((0.5:1:u(3)-0.5)'*(pi/u(3)),p3-1))*(2/u(3))^0.5;
378: V3(:,p3==1) = 1/u(3)^0.5;
379: elseif ((BB(3) == 1) && (BB(3+dim) == 2))
380: V3 = sin(kron((1:u(3))'*(pi/2/(u(3)+0.5)),2*p3 - 1))...
381: *(2/(u(3)+0.5))^0.5;
382: elseif ((BB(3) == 2) && (BB(3+dim) == 1))
383: V3 = cos(kron((0.5:1:u(3)-0.5)'*(pi/2/(u(3)+0.5)),2*p3 - 1))...
384: *(2/(u(3)+0.5))^0.5;
385: else
386: V3 = zeros(u(3),m);
387: a = (0.5:1:u(3)-0.5)';
388: V3(:,mod(p3,2)==1) = cos(a*(pi/u(3)*(p3(mod(p3,2)==1)-1)))...
389: *(2/u(3))^0.5;
390: pp = p1(mod(p3,2)==0);
391: if ~isempty(pp)
392: V3(:,mod(p3,2)==0) = sin(a*(pi/u(3)*p3(mod(p3,2)==0)))...
393: *(2/u(3))^0.5;
394: end
395: V3(:,p3==1) = 1/u(3)^0.5;
396: if mod(u(3),2) == 0
397: V3(:,p3==u(3)) = V3(:,p3==u(3))/2^0.5;
398: end
399:
400: end
401: else
402: V3 = ones(1,m);
403: end
404:
405: if dim == 1
406: V = V1;
407: elseif dim == 2
408: V = kron(e2,V1).*kron(V2,e1);
409: else
410: V = kron(e3, kron(e2, V1)).*kron(e3, kron(V2, e1))...
411: .*kron(kron(V3,e2),e1);
412: end
413:
414: if m ~= 0
415: if abs(lambda(m) - w) < maxeigs*eps('double')
416: sprintf('\n%s','Warning: (m+1)th eigenvalue is nearly equal',...
417: ' to mth.')
418:
419: end
420: end
421:
422: end
424: A = [];
425: if nargout > 2 %also calulate the matrix if specified in the output
426:
427: % Set the component matrices. SPDIAGS converts int8 into double anyway.
428: % e1 = ones(u(1),1); %e1 = ones(u(1),1,'int8');
429: D1x = spdiags([-e1 2*e1 -e1], [-1 0 1], u(1),u(1));
430: if dim > 1
431: % e2 = ones(u(2),1);
432: D1y = spdiags([-e2 2*e2 -e2], [-1 0 1], u(2),u(2));
433: end
434: if dim > 2
435: % e3 = ones(u(3),1);
436: D1z = spdiags([-e3 2*e3 -e3], [-1 0 1], u(3),u(3));
437: end
438:
439:
440: % Set boundary conditions if other than Dirichlet.
441: for i = 1:dim
442: if BB(i) == 2
443: eval(['D1' char(119 + i) '(1,1) = 1;'])
444: elseif BB(i) == 3
445: eval(['D1' char(119 + i) '(1,' num2str(u(i)) ') = D1'...
446: char(119 + i) '(1,' num2str(u(i)) ') -1;']);
447: eval(['D1' char(119 + i) '(' num2str(u(i)) ',1) = D1'...
448: char(119 + i) '(' num2str(u(i)) ',1) -1;']);
449: end
450:
451: if BB(i+dim) == 2
452: eval(['D1' char(119 + i)...
453: '(',num2str(u(i)),',',num2str(u(i)),') = 1;'])
454: end
455: end
456:
457: % Form A using tensor products of lower dimensional Laplacians
458: Ix = speye(u(1));
459: if dim == 1
460: A = D1x;
461: elseif dim == 2
462: Iy = speye(u(2));
463: A = kron(Iy,D1x) + kron(D1y,Ix);
464: elseif dim == 3
465: Iy = speye(u(2));
466: Iz = speye(u(3));
467: A = kron(Iz, kron(Iy, D1x)) + kron(Iz, kron(D1y, Ix))...
468: + kron(kron(D1z,Iy),Ix);
469: end
470: end
472: disp(' ')
473: toc
474: if ~isempty(V)
475: a = whos('regep','V');
476: disp(['The eigenvectors take ' num2str(a.bytes) ' bytes'])
477: end
478: if ~isempty(A)
479: a = whos('regexp','A');
480: disp(['The Laplacian matrix takes ' num2str(a.bytes) ' bytes'])
481: end
482: disp(' ')