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. 4## 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. 15## 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. 26 27% LAPLACIAN Sparse Negative Laplacian in 1D, 2D, or 3D 28% 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. 34% 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'. 40% 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. 46% 47% [LAMBDA,V] = LAPLACIAN(N,B,M) also outputs orthonormal eigenvectors 48% associated with the corresponding m smallest eigenvalues. 49% 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. 53% 54% % Examples: 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); 58% % or 59% lambda = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20); 60% % computes the eigenvalues only 61% 62% [~,V,~] = laplacian([200 200],{'DD' 'DN'},30); 63% % Eigenvectors of 2D negative Laplacian with mixed boundary conditions. 64% 65% [~,~,A] = laplacian(200,{'DN'},30); 66% % 1D negative Laplacian matrix A with mixed boundary conditions. 67% 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 75% 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]); 80% norm(A1-A2,inf) 81% 82% Class support for inputs: 83% N - row vector float double 84% B - cell array 85% M - scalar float double 86% 87% Class support for outputs: 88% lambda and V - full float double, A - sparse float double. 89% 90% Note: the actual numerical entries of A fit int8 format, but only 91% double data class is currently (2010) supported for sparse matrices. 92% 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. 101% 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/ 105 106% Revision 1.1 changes: rearranged the output variables, always compute 107% the eigenvalues, compute eigenvectors and/or the matrix on demand only. 108 109% $Revision: 1.1 $ $Date: 1-Aug-2011 110% Tested in MATLAB 7.11.0 (R2010b) and Octave 3.4.0. 111 112function [lambda, V, A] = laplacian(varargin) 113 114 % Input/Output handling. 115 if (nargin < 1 || nargin > 3) 116 print_usage; 117 endif 118 119 u = varargin{1}; 120 dim2 = size(u); 121 if dim2(1) ~= 1 122 error('BLOPEX:laplacian:WrongVectorOfGridPoints',... 123 '%s','Number of grid points must be in a row vector.') 124 end 125 if dim2(2) > 3 126 error('BLOPEX:laplacian:WrongNumberOfGridPoints',... 127 '%s','Number of grid points must be a 1, 2, or 3') 128 end 129 dim=dim2(2); clear dim2; 130 131 uint = round(u); 132 if max(uint~=u) 133 warning('BLOPEX:laplacian:NonIntegerGridSize',... 134 '%s','Grid sizes must be integers. Rounding...'); 135 u = uint; clear uint 136 end 137 if max(u<=0 ) 138 error('BLOPEX:laplacian:NonIntegerGridSize',... 139 '%s','Grid sizes must be positive.'); 140 end 141 142 if nargin == 3 143 m = varargin{3}; 144 B = varargin{2}; 145 elseif nargin == 2 146 f = varargin{2}; 147 a = whos('regep','f'); 148 if sum(a.class(1:4)=='cell') == 4 149 B = f; 150 m = 0; 151 elseif sum(a.class(1:4)=='doub') == 4 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 = f; 160 else 161 error('BLOPEX:laplacian:InvalidClass',... 162 '%s','Second input must be either class double or a cell array.'); 163 end 164 else 165 if dim == 1 166 B = {'DD'}; 167 elseif dim == 2 168 B = {'DD' 'DD'}; 169 else 170 B = {'DD' 'DD' 'DD'}; 171 end 172 m = 0; 173 end 174 175 if max(size(m) - [1 1]) ~= 0 176 error('BLOPEX:laplacian:WrongNumberOfEigenvalues',... 177 '%s','The requested number of eigenvalues must be a scalar.'); 178 end 179 180 maxeigs = prod(u); 181 mint = round(m); 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.'); 186 end 187 m = mint; 188 189 bdryerr = 0; 190 a = whos('regep','B'); 191 if sum(a.class(1:4)=='cell') ~= 4 || sum(a.size == [1 dim]) ~= 2 192 bdryerr = 1; 193 else 194 BB = zeros(1, 2*dim); 195 for i = 1:dim 196 if (length(B{i}) == 1) 197 if B{i} == 'P' 198 BB(i) = 3; 199 BB(i + dim) = 3; 200 else 201 bdryerr = 1; 202 end 203 elseif (length(B{i}) == 2) 204 if B{i}(1) == 'D' 205 BB(i) = 1; 206 elseif B{i}(1) == 'N' 207 BB(i) = 2; 208 else 209 bdryerr = 1; 210 end 211 if B{i}(2) == 'D' 212 BB(i + dim) = 1; 213 elseif B{i}(2) == 'N' 214 BB(i + dim) = 2; 215 else 216 bdryerr = 1; 217 end 218 else 219 bdryerr = 1; 220 end 221 end 222 end 223 224 if bdryerr == 1 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''',... 228 ', or ''P''.'); 229 end 230 231 % Set the component matrices. SPDIAGS converts int8 into double anyway. 232 e1 = ones(u(1),1); %e1 = ones(u(1),1,'int8'); 233 if dim > 1 234 e2 = ones(u(2),1); 235 end 236 if dim > 2 237 e3 = ones(u(3),1); 238 end 239 240 % Calculate m smallest exact eigenvalues. 241 if m > 0 242 if (BB(1) == 1) && (BB(1+dim) == 1) 243 a1 = pi/2/(u(1)+1); 244 N = (1:u(1))'; 245 elseif (BB(1) == 2) && (BB(1+dim) == 2) 246 a1 = pi/2/u(1); 247 N = (0:(u(1)-1))'; 248 elseif ((BB(1) == 1) && (BB(1+dim) == 2)) || ((BB(1) == 2)... 249 && (BB(1+dim) == 1)) 250 a1 = pi/4/(u(1)+0.5); 251 N = 2*(1:u(1))' - 1; 252 else 253 a1 = pi/u(1); 254 N = floor((1:u(1))/2)'; 255 end 256 257 lambda1 = 4*sin(a1*N).^2; 258 259 if dim > 1 260 if (BB(2) == 1) && (BB(2+dim) == 1) 261 a2 = pi/2/(u(2)+1); 262 N = (1:u(2))'; 263 elseif (BB(2) == 2) && (BB(2+dim) == 2) 264 a2 = pi/2/u(2); 265 N = (0:(u(2)-1))'; 266 elseif ((BB(2) == 1) && (BB(2+dim) == 2)) || ((BB(2) == 2)... 267 && (BB(2+dim) == 1)) 268 a2 = pi/4/(u(2)+0.5); 269 N = 2*(1:u(2))' - 1; 270 else 271 a2 = pi/u(2); 272 N = floor((1:u(2))/2)'; 273 end 274 lambda2 = 4*sin(a2*N).^2; 275 else 276 lambda2 = 0; 277 end 278 279 if dim > 2 280 if (BB(3) == 1) && (BB(6) == 1) 281 a3 = pi/2/(u(3)+1); 282 N = (1:u(3))'; 283 elseif (BB(3) == 2) && (BB(6) == 2) 284 a3 = pi/2/u(3); 285 N = (0:(u(3)-1))'; 286 elseif ((BB(3) == 1) && (BB(6) == 2)) || ((BB(3) == 2)... 287 && (BB(6) == 1)) 288 a3 = pi/4/(u(3)+0.5); 289 N = 2*(1:u(3))' - 1; 290 else 291 a3 = pi/u(3); 292 N = floor((1:u(3))/2)'; 293 end 294 lambda3 = 4*sin(a3*N).^2; 295 else 296 lambda3 = 0; 297 end 298 299 if dim == 1 300 lambda = lambda1; 301 elseif dim == 2 302 lambda = kron(e2,lambda1) + kron(lambda2, e1); 303 else 304 lambda = kron(e3,kron(e2,lambda1)) + kron(e3,kron(lambda2,e1))... 305 + kron(lambda3,kron(e2,e1)); 306 end 307 [lambda, p] = sort(lambda); 308 if m < maxeigs - 0.1 309 w = lambda(m+1); 310 else 311 w = inf; 312 end 313 lambda = lambda(1:m); 314 p = p(1:m)'; 315 else 316 lambda = []; 317 end 318 319 V = []; 320 if nargout > 1 && m > 0 % Calculate eigenvectors if specified in output. 321 322 p1 = mod(p-1,u(1))+1; 323 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))... 331 *(2/(u(1)+0.5))^0.5; 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))... 334 *(2/(u(1)+0.5))^0.5; 335 else 336 V1 = zeros(u(1),m); 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)))... 339 *(2/u(1))^0.5; 340 pp = p1(mod(p1,2)==0); 341 if ~isempty(pp) 342 V1(:,mod(p1,2)==0) = sin(a*(pi/u(1)*p1(mod(p1,2)==0)))... 343 *(2/u(1))^0.5; 344 end 345 V1(:,p1==1) = 1/u(1)^0.5; 346 if mod(u(1),2) == 0 347 V1(:,p1==u(1)) = V1(:,p1==u(1))/2^0.5; 348 end 349 end 350 351 352 if dim > 1 353 p2 = mod(p-p1,u(1)*u(2)); 354 p3 = (p - p2 - p1)/(u(1)*u(2)) + 1; 355 p2 = p2/u(1) + 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))... 363 *(2/(u(2)+0.5))^0.5; 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))... 366 *(2/(u(2)+0.5))^0.5; 367 else 368 V2 = zeros(u(2),m); 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)))... 371 *(2/u(2))^0.5; 372 pp = p2(mod(p2,2)==0); 373 if ~isempty(pp) 374 V2(:,mod(p2,2)==0) = sin(a*(pi/u(2)*p2(mod(p2,2)==0)))... 375 *(2/u(2))^0.5; 376 end 377 V2(:,p2==1) = 1/u(2)^0.5; 378 if mod(u(2),2) == 0 379 V2(:,p2==u(2)) = V2(:,p2==u(2))/2^0.5; 380 end 381 end 382 else 383 V2 = ones(1,m); 384 end 385 386 if dim > 2 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))... 394 *(2/(u(3)+0.5))^0.5; 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))... 397 *(2/(u(3)+0.5))^0.5; 398 else 399 V3 = zeros(u(3),m); 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)))... 402 *(2/u(3))^0.5; 403 pp = p1(mod(p3,2)==0); 404 if ~isempty(pp) 405 V3(:,mod(p3,2)==0) = sin(a*(pi/u(3)*p3(mod(p3,2)==0)))... 406 *(2/u(3))^0.5; 407 end 408 V3(:,p3==1) = 1/u(3)^0.5; 409 if mod(u(3),2) == 0 410 V3(:,p3==u(3)) = V3(:,p3==u(3))/2^0.5; 411 end 412 413 end 414 else 415 V3 = ones(1,m); 416 end 417 418 if dim == 1 419 V = V1; 420 elseif dim == 2 421 V = kron(e2,V1).*kron(V2,e1); 422 else 423 V = kron(e3, kron(e2, V1)).*kron(e3, kron(V2, e1))... 424 .*kron(kron(V3,e2),e1); 425 end 426 427 if m ~= 0 428 if abs(lambda(m) - w) < maxeigs*eps('double') 429 sprintf('\n%s','Warning: (m+1)th eigenvalue is nearly equal',... 430 ' to mth.') 431 432 end 433 end 434 435 end 436 437 A = []; 438 if nargout > 2 %also calulate the matrix if specified in the output 439 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)); 443 if dim > 1 444 % e2 = ones(u(2),1); 445 D1y = spdiags([-e2 2*e2 -e2], [-1 0 1], u(2),u(2)); 446 end 447 if dim > 2 448 % e3 = ones(u(3),1); 449 D1z = spdiags([-e3 2*e3 -e3], [-1 0 1], u(3),u(3)); 450 end 451 452 453 % Set boundary conditions if other than Dirichlet. 454 for i = 1:dim 455 if BB(i) == 2 456 eval(['D1' char(119 + i) '(1,1) = 1;']) 457 elseif BB(i) == 3 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;']); 462 end 463 464 if BB(i+dim) == 2 465 eval(['D1' char(119 + i)... 466 '(',num2str(u(i)),',',num2str(u(i)),') = 1;']) 467 end 468 end 469 470 % Form A using tensor products of lower dimensional Laplacians 471 Ix = speye(u(1)); 472 if dim == 1 473 A = D1x; 474 elseif dim == 2 475 Iy = speye(u(2)); 476 A = kron(Iy,D1x) + kron(D1y,Ix); 477 elseif dim == 3 478 Iy = speye(u(2)); 479 Iz = speye(u(3)); 480 A = kron(Iz, kron(Iy, D1x)) + kron(Iz, kron(D1y, Ix))... 481 + kron(kron(D1z,Iy),Ix); 482 end 483 end 484 485 disp(' ') 486 if ~isempty(V) 487 a = whos('regep','V'); 488 disp(['The eigenvectors take ' num2str(a.bytes) ' bytes']) 489 end 490 if ~isempty(A) 491 a = whos('regexp','A'); 492 disp(['The Laplacian matrix takes ' num2str(a.bytes) ' bytes']) 493 end 494 disp(' ') 495endfunction 496