1function [lambda, V, A] = laplacian(varargin) 2 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/ 81 82% Revision 1.1 changes: rearranged the output variables, always compute 83% the eigenvalues, compute eigenvectors and/or the matrix on demand only. 84 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. 89 90tic 91 92% Input/Output handling. 93if nargin > 3 94 error('BLOPEX:laplacian:TooManyInputs',... 95 '%s','Too many input arguments.'); 96elseif nargin == 0 97 error('BLOPEX:laplacian:NoInputArguments',... 98 '%s','Must have at least one input argument.'); 99end 100 101if nargout > 3 102 error('BLOPEX:laplacian:TooManyOutputs',... 103 '%s','Maximum number of outputs is 3.'); 104end 105 106u = varargin{1}; 107dim2 = size(u); 108if dim2(1) ~= 1 109 error('BLOPEX:laplacian:WrongVectorOfGridPoints',... 110 '%s','Number of grid points must be in a row vector.') 111end 112if dim2(2) > 3 113 error('BLOPEX:laplacian:WrongNumberOfGridPoints',... 114 '%s','Number of grid points must be a 1, 2, or 3') 115end 116dim=dim2(2); clear dim2; 117 118uint = round(u); 119if max(uint~=u) 120 warning('BLOPEX:laplacian:NonIntegerGridSize',... 121 '%s','Grid sizes must be integers. Rounding...'); 122 u = uint; clear uint 123end 124if max(u<=0 ) 125 error('BLOPEX:laplacian:NonIntegerGridSize',... 126 '%s','Grid sizes must be positive.'); 127end 128 129if nargin == 3 130 m = varargin{3}; 131 B = varargin{2}; 132elseif 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 151else 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; 160end 161 162if max(size(m) - [1 1]) ~= 0 163 error('BLOPEX:laplacian:WrongNumberOfEigenvalues',... 164 '%s','The requested number of eigenvalues must be a scalar.'); 165end 166 167maxeigs = prod(u); 168mint = round(m); 169if 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.'); 173end 174m = mint; 175 176bdryerr = 0; 177a = whos('regep','B'); 178if sum(a.class(1:4)=='cell') ~= 4 || sum(a.size == [1 dim]) ~= 2 179 bdryerr = 1; 180else 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 209end 210 211if 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''.'); 216end 217 218% Set the component matrices. SPDIAGS converts int8 into double anyway. 219e1 = ones(u(1),1); %e1 = ones(u(1),1,'int8'); 220if dim > 1 221 e2 = ones(u(2),1); 222end 223if dim > 2 224 e3 = ones(u(3),1); 225end 226 227% Calculate m smallest exact eigenvalues. 228if 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)'; 302else 303 lambda = []; 304end 305 306V = []; 307if 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 422end 423 424A = []; 425if nargout > 2 %also calculate 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 470end 471 472disp(' ') 473toc 474if ~isempty(V) 475 a = whos('regep','V'); 476 disp(['The eigenvectors take ' num2str(a.bytes) ' bytes']) 477end 478if ~isempty(A) 479 a = whos('regexp','A'); 480 disp(['The Laplacian matrix takes ' num2str(a.bytes) ' bytes']) 481end 482disp(' ') 483 484