Actual source code: laplacian.m

petsc-3.14.5 2021-03-03
Report Typos and Errors
  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('  ')