1 ## Copyright (c) 2010-2011 Andrew V. Knyazev <andrew.knyazev@ucdenver.edu>
2 ## Copyright (c) 2010-2011 Bryan C. Smith <bryan.c.smith@ucdenver.edu>
3 ## All rights reserved.
5 ## Redistribution and use in source and binary forms, with or without
6 ## modification, are permitted provided that the following conditions are met:
7 ## * Redistributions of source code must retain the above copyright
8 ## notice, this list of conditions and the following disclaimer.
9 ## * Redistributions in binary form must reproduce the above copyright
10 ## notice, this list of conditions and the following disclaimer in the
11 ## documentation and/or other materials provided with the distribution.
12 ## * Neither the name of the <organization> nor the
13 ## names of its contributors may be used to endorse or promote products
14 ## derived from this software without specific prior written permission.
16 ## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17 ## ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18 ## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19 ## DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
20 ## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21 ## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22 ## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23 ## ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 ## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25 ## SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 % LAPLACIAN Sparse Negative Laplacian in 1D, 2D, or 3D
29 % [~,~,A]=LAPLACIAN(N) generates a sparse negative 3D Laplacian matrix
30 % with Dirichlet boundary conditions, from a rectangular cuboid regular
31 % grid with j x k x l interior grid points if N = [j k l], using the
32 % standard 7-point finite-difference scheme, The grid size is always
33 % one in all directions.
35 % [~,~,A]=LAPLACIAN(N,B) specifies boundary conditions with a cell array
36 % B. For example, B = {'DD' 'DN' 'P'} will Dirichlet boundary conditions
37 % ('DD') in the x-direction, Dirichlet-Neumann conditions ('DN') in the
38 % y-direction and period conditions ('P') in the z-direction. Possible
39 % values for the elements of B are 'DD', 'DN', 'ND', 'NN' and 'P'.
41 % LAMBDA = LAPLACIAN(N,B,M) or LAPLACIAN(N,M) outputs the m smallest
42 % eigenvalues of the matrix, computed by an exact known formula, see
43 % http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_of_the_second_derivative
44 % It will produce a warning if the mth eigenvalue is equal to the
45 % (m+1)th eigenvalue. If m is absebt or zero, lambda will be empty.
47 % [LAMBDA,V] = LAPLACIAN(N,B,M) also outputs orthonormal eigenvectors
48 % associated with the corresponding m smallest eigenvalues.
50 % [LAMBDA,V,A] = LAPLACIAN(N,B,M) produces a 2D or 1D negative
51 % Laplacian matrix if the length of N and B are 2 or 1 respectively.
52 % It uses the standard 5-point scheme for 2D, and 3-point scheme for 1D.
55 % [lambda,V,A] = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
56 % % Everything for 3D negative Laplacian with mixed boundary conditions.
57 % laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
59 % lambda = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
60 % % computes the eigenvalues only
62 % [~,V,~] = laplacian([200 200],{'DD' 'DN'},30);
63 % % Eigenvectors of 2D negative Laplacian with mixed boundary conditions.
65 % [~,~,A] = laplacian(200,{'DN'},30);
66 % % 1D negative Laplacian matrix A with mixed boundary conditions.
68 % % Example to test if outputs correct eigenvalues and vectors:
69 % [lambda,V,A] = laplacian([13,10,6],{'DD' 'DN' 'P'},30);
70 % [Veig D] = eig(full(A)); lambdaeig = diag(D(1:30,1:30));
71 % max(abs(lambda-lambdaeig)) %checking eigenvalues
72 % subspace(V,Veig(:,1:30)) %checking the invariant subspace
73 % subspace(V(:,1),Veig(:,1)) %checking selected eigenvectors
74 % subspace(V(:,29:30),Veig(:,29:30)) %a multiple eigenvalue
76 % % Example showing equivalence between laplacian.m and built-in MATLAB
77 % % DELSQ for the 2D case. The output of the last command shall be 0.
78 % A1 = delsq(numgrid('S',32)); % input 'S' specifies square grid.
79 % [~,~,A2] = laplacian([30,30]);
82 % Class support for inputs:
83 % N - row vector float double
85 % M - scalar float double
87 % Class support for outputs:
88 % lambda and V - full float double, A - sparse float double.
90 % Note: the actual numerical entries of A fit int8 format, but only
91 % double data class is currently (2010) supported for sparse matrices.
93 % This program is designed to efficiently compute eigenvalues,
94 % eigenvectors, and the sparse matrix of the (1-3)D negative Laplacian
95 % on a rectangular grid for Dirichlet, Neumann, and Periodic boundary
96 % conditions using tensor sums of 1D Laplacians. For more information on
97 % tensor products, see
98 % http://en.wikipedia.org/wiki/Kronecker_sum_of_discrete_Laplacians
99 % For 2D case in MATLAB, see
100 % http://www.mathworks.com/access/helpdesk/help/techdoc/ref/kron.html.
102 % This code is also part of the BLOPEX package:
103 % http://en.wikipedia.org/wiki/BLOPEX or directly
104 % http://code.google.com/p/blopex/
106 % Revision 1.1 changes: rearranged the output variables, always compute
107 % the eigenvalues, compute eigenvectors and/or the matrix on demand only.
109 % $Revision: 1.1 $ $Date: 1-Aug-2011
110 % Tested in MATLAB 7.11.0 (R2010b) and Octave 3.4.0.
112 function [lambda, V, A] = laplacian(varargin)
114 % Input/Output handling.
115 if (nargin < 1 || nargin > 3)
122 error('BLOPEX:laplacian:WrongVectorOfGridPoints',...
123 '%s','Number of grid points must be in a row vector.')
126 error('BLOPEX:laplacian:WrongNumberOfGridPoints',...
127 '%s','Number of grid points must be a 1, 2, or 3')
129 dim=dim2(2); clear dim2;
133 warning('BLOPEX:laplacian:NonIntegerGridSize',...
134 '%s','Grid sizes must be integers. Rounding...');
138 error('BLOPEX:laplacian:NonIntegerGridSize',...
139 '%s','Grid sizes must be positive.');
147 a = whos('regep','f');
148 if sum(a.class(1:4)=='cell') == 4
151 elseif sum(a.class(1:4)=='doub') == 4
157 B = {'DD' 'DD' 'DD'};
161 error('BLOPEX:laplacian:InvalidClass',...
162 '%s','Second input must be either class double or a cell array.');
170 B = {'DD' 'DD' 'DD'};
175 if max(size(m) - [1 1]) ~= 0
176 error('BLOPEX:laplacian:WrongNumberOfEigenvalues',...
177 '%s','The requested number of eigenvalues must be a scalar.');
182 if mint ~= m || mint > maxeigs
183 error('BLOPEX:laplacian:InvalidNumberOfEigs',...
184 '%s','Number of eigenvalues output must be a nonnegative ',...
185 'integer no bigger than number of grid points.');
190 a = whos('regep','B');
191 if sum(a.class(1:4)=='cell') ~= 4 || sum(a.size == [1 dim]) ~= 2
194 BB = zeros(1, 2*dim);
196 if (length(B{i}) == 1)
203 elseif (length(B{i}) == 2)
206 elseif B{i}(1) == 'N'
213 elseif B{i}(2) == 'N'
225 error('BLOPEX:laplacian:InvalidBdryConds',...
226 '%s','Boundary conditions must be a cell of length 3 for 3D, 2',...
227 ' for 2D, 1 for 1D, with values ''DD'', ''DN'', ''ND'', ''NN''',...
231 % Set the component matrices. SPDIAGS converts int8 into double anyway.
232 e1 = ones(u(1),1); %e1 = ones(u(1),1,'int8');
240 % Calculate m smallest exact eigenvalues.
242 if (BB(1) == 1) && (BB(1+dim) == 1)
245 elseif (BB(1) == 2) && (BB(1+dim) == 2)
248 elseif ((BB(1) == 1) && (BB(1+dim) == 2)) || ((BB(1) == 2)...
250 a1 = pi/4/(u(1)+0.5);
254 N = floor((1:u(1))/2)';
257 lambda1 = 4*sin(a1*N).^2;
260 if (BB(2) == 1) && (BB(2+dim) == 1)
263 elseif (BB(2) == 2) && (BB(2+dim) == 2)
266 elseif ((BB(2) == 1) && (BB(2+dim) == 2)) || ((BB(2) == 2)...
268 a2 = pi/4/(u(2)+0.5);
272 N = floor((1:u(2))/2)';
274 lambda2 = 4*sin(a2*N).^2;
280 if (BB(3) == 1) && (BB(6) == 1)
283 elseif (BB(3) == 2) && (BB(6) == 2)
286 elseif ((BB(3) == 1) && (BB(6) == 2)) || ((BB(3) == 2)...
288 a3 = pi/4/(u(3)+0.5);
292 N = floor((1:u(3))/2)';
294 lambda3 = 4*sin(a3*N).^2;
302 lambda = kron(e2,lambda1) + kron(lambda2, e1);
304 lambda = kron(e3,kron(e2,lambda1)) + kron(e3,kron(lambda2,e1))...
305 + kron(lambda3,kron(e2,e1));
307 [lambda, p] = sort(lambda);
313 lambda = lambda(1:m);
320 if nargout > 1 && m > 0 % Calculate eigenvectors if specified in output.
322 p1 = mod(p-1,u(1))+1;
324 if (BB(1) == 1) && (BB(1+dim) == 1)
325 V1 = sin(kron((1:u(1))'*(pi/(u(1)+1)),p1))*(2/(u(1)+1))^0.5;
326 elseif (BB(1) == 2) && (BB(1+dim) == 2)
327 V1 = cos(kron((0.5:1:u(1)-0.5)'*(pi/u(1)),p1-1))*(2/u(1))^0.5;
328 V1(:,p1==1) = 1/u(1)^0.5;
329 elseif ((BB(1) == 1) && (BB(1+dim) == 2))
330 V1 = sin(kron((1:u(1))'*(pi/2/(u(1)+0.5)),2*p1 - 1))...
332 elseif ((BB(1) == 2) && (BB(1+dim) == 1))
333 V1 = cos(kron((0.5:1:u(1)-0.5)'*(pi/2/(u(1)+0.5)),2*p1 - 1))...
337 a = (0.5:1:u(1)-0.5)';
338 V1(:,mod(p1,2)==1) = cos(a*(pi/u(1)*(p1(mod(p1,2)==1)-1)))...
340 pp = p1(mod(p1,2)==0);
342 V1(:,mod(p1,2)==0) = sin(a*(pi/u(1)*p1(mod(p1,2)==0)))...
345 V1(:,p1==1) = 1/u(1)^0.5;
347 V1(:,p1==u(1)) = V1(:,p1==u(1))/2^0.5;
353 p2 = mod(p-p1,u(1)*u(2));
354 p3 = (p - p2 - p1)/(u(1)*u(2)) + 1;
356 if (BB(2) == 1) && (BB(2+dim) == 1)
357 V2 = sin(kron((1:u(2))'*(pi/(u(2)+1)),p2))*(2/(u(2)+1))^0.5;
358 elseif (BB(2) == 2) && (BB(2+dim) == 2)
359 V2 = cos(kron((0.5:1:u(2)-0.5)'*(pi/u(2)),p2-1))*(2/u(2))^0.5;
360 V2(:,p2==1) = 1/u(2)^0.5;
361 elseif ((BB(2) == 1) && (BB(2+dim) == 2))
362 V2 = sin(kron((1:u(2))'*(pi/2/(u(2)+0.5)),2*p2 - 1))...
364 elseif ((BB(2) == 2) && (BB(2+dim) == 1))
365 V2 = cos(kron((0.5:1:u(2)-0.5)'*(pi/2/(u(2)+0.5)),2*p2 - 1))...
369 a = (0.5:1:u(2)-0.5)';
370 V2(:,mod(p2,2)==1) = cos(a*(pi/u(2)*(p2(mod(p2,2)==1)-1)))...
372 pp = p2(mod(p2,2)==0);
374 V2(:,mod(p2,2)==0) = sin(a*(pi/u(2)*p2(mod(p2,2)==0)))...
377 V2(:,p2==1) = 1/u(2)^0.5;
379 V2(:,p2==u(2)) = V2(:,p2==u(2))/2^0.5;
387 if (BB(3) == 1) && (BB(3+dim) == 1)
388 V3 = sin(kron((1:u(3))'*(pi/(u(3)+1)),p3))*(2/(u(3)+1))^0.5;
389 elseif (BB(3) == 2) && (BB(3+dim) == 2)
390 V3 = cos(kron((0.5:1:u(3)-0.5)'*(pi/u(3)),p3-1))*(2/u(3))^0.5;
391 V3(:,p3==1) = 1/u(3)^0.5;
392 elseif ((BB(3) == 1) && (BB(3+dim) == 2))
393 V3 = sin(kron((1:u(3))'*(pi/2/(u(3)+0.5)),2*p3 - 1))...
395 elseif ((BB(3) == 2) && (BB(3+dim) == 1))
396 V3 = cos(kron((0.5:1:u(3)-0.5)'*(pi/2/(u(3)+0.5)),2*p3 - 1))...
400 a = (0.5:1:u(3)-0.5)';
401 V3(:,mod(p3,2)==1) = cos(a*(pi/u(3)*(p3(mod(p3,2)==1)-1)))...
403 pp = p1(mod(p3,2)==0);
405 V3(:,mod(p3,2)==0) = sin(a*(pi/u(3)*p3(mod(p3,2)==0)))...
408 V3(:,p3==1) = 1/u(3)^0.5;
410 V3(:,p3==u(3)) = V3(:,p3==u(3))/2^0.5;
421 V = kron(e2,V1).*kron(V2,e1);
423 V = kron(e3, kron(e2, V1)).*kron(e3, kron(V2, e1))...
424 .*kron(kron(V3,e2),e1);
428 if abs(lambda(m) - w) < maxeigs*eps('double')
429 sprintf('\n%s','Warning: (m+1)th eigenvalue is nearly equal',...
438 if nargout > 2 %also calulate the matrix if specified in the output
440 % Set the component matrices. SPDIAGS converts int8 into double anyway.
441 % e1 = ones(u(1),1); %e1 = ones(u(1),1,'int8');
442 D1x = spdiags([-e1 2*e1 -e1], [-1 0 1], u(1),u(1));
445 D1y = spdiags([-e2 2*e2 -e2], [-1 0 1], u(2),u(2));
449 D1z = spdiags([-e3 2*e3 -e3], [-1 0 1], u(3),u(3));
453 % Set boundary conditions if other than Dirichlet.
456 eval(['D1' char(119 + i) '(1,1) = 1;'])
458 eval(['D1' char(119 + i) '(1,' num2str(u(i)) ') = D1'...
459 char(119 + i) '(1,' num2str(u(i)) ') -1;']);
460 eval(['D1' char(119 + i) '(' num2str(u(i)) ',1) = D1'...
461 char(119 + i) '(' num2str(u(i)) ',1) -1;']);
465 eval(['D1' char(119 + i)...
466 '(',num2str(u(i)),',',num2str(u(i)),') = 1;'])
470 % Form A using tensor products of lower dimensional Laplacians
476 A = kron(Iy,D1x) + kron(D1y,Ix);
480 A = kron(Iz, kron(Iy, D1x)) + kron(Iz, kron(D1y, Ix))...
481 + kron(kron(D1z,Iy),Ix);
487 a = whos('regep','V');
488 disp(['The eigenvectors take ' num2str(a.bytes) ' bytes'])
491 a = whos('regexp','A');
492 disp(['The Laplacian matrix takes ' num2str(a.bytes) ' bytes'])