octave_packages/specfun-1.1.0/laplacian.m

   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
 112 function [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('  ')
 495 endfunction