function varargout = laplacian2(numpts, varargin)
% LAPLACIAN2 Sparse Negative Laplacian in 1D, 2D, or 3D
%
% A = laplacian2(N)
% % Generates a sparse 3D Laplacian matrix for the finite difference method
% % with Dirichlet boundary conditions, from a rectangular regular grid
% % with j x k x l interior grid points where N = [j k l], using the
% % standard 7-point finite-difference scheme, The grid spacing h = 1.
% % If N has length 1 or 2, then 1D and 2D Laplacians are constructed).
%
% options.mass_matrix = true;
% [L, M] = laplacian2(N,options);
% % Generates a sparse 3D Laplacian stiffness and mass matrix for the
% % finite element method using bilinear elements with Dirichlet
% % boundary conditions.
%
% options.mass_matrix = false;
% options.num_eigs = m;
% [A, lambda] = laplacian2(N,options)
% % Also outputs the first m eigenvalues in ascending order if m <=
% % j*k*l.
%
% options.mass_matrix = false;
% options.num_eigs = m;
% [A, lambda, V] = laplacian2(N,options)
% % Also outputs the eigenvectors V.
%
% options.mass_matrix = true;
% options.num_eigs = m;
% [L, M, lambda] = laplacian2(N,options)
% % Outputs the first m eigenvalues for the finite element method
%
% options.mass_matrix = true;
% options.num_eigs = m;
% [L, M, lambda, V] = laplacian2(N,options)
% % Also outputs the eigenvectors V.
%
% options.mass_matrix = true;
% options.bdry_conds = {'DD' 'DN' 'P'}
% options.numeigs = m;
% [L, M, lambda, V] = laplacian2(N,options);
% % Specify the boundary conditions in the x, y, and z directions with a
% % cell array options.bdry_conds. For example,
% % options.bdry_conds = {'DD' 'DN' 'P'} will Dirichlet boundary conditions
% % ('DD') in the x-direction, Dirichlet-Neumann conditions ('DN') in
% % the y-direction and period conditions ('P') in the z-direction.
% % Possible values for the elements of options.bdry_conds are 'DD',
% %'DN', 'ND', 'NN' and 'P'.
%
% A = LAPLACIAN(N,options) produces a 2D or 1D Laplacian matrices
% if the length of N and options.bdry_conds are 2 or 1 respectively.
% It uses the standard 5-point scheme for 2D, and 3-point scheme for 1D.
%
% %%% Examples for the finite difference method: %%%
%
% options.bdry_conds = {'DD' 'NN' 'P'};
% options.numeigs = 30;
% [A,lambda,V] = laplacian2([100,45,55],options);
% % 3D negative Laplacian with mixed boundary conditions.
% options.bdry_conds = {'DD' 'DN'};
% options.numeigs = 30;
% [A,lambda,V] = laplacian2([200 200],options);
% % 2D negative Laplacian with mixed boundary conditions.
% options.bdry_conds = {'DN'};
% options.num_eigs = 30;
% [A,lambda,V] = laplacian2(200,options);
% % 1D negative Laplacian with mixed boundary conditions.
%
% % Example to test if outputs correct eigenvalues and vectors:
% options.bdry_conds = {'DD' 'DN' 'P'};
% options.numeigs = 30;
% [A,lambda,V] = laplacian2([13,10,6],options);
% [Veig D] = eig(full(A)); lambdaeig = diag(D(1:30,1:30));
% max(abs(lambda-lambdaeig)) %checking eigenvalues
% subspace(V,Veig(:,1:30)) %checking the invariant subspace
% subspace(V(:,1),Veig(:,1)) %checking selected eigenvectors
% subspace(V(:,29:30),Veig(:,29:30)) %a multiple eigenvalue
%
% % Example showing equivalence between laplacian.m and built-in MATLAB
% % DELSQ for the 2D case. The output of the last command shall be 0.
% A1 = delsq(numgrid('S',32)); % input 'S' specifies square grid.
% A2 = laplacian2([30,30]);
% norm(A1-A2,inf)
% %%% Examples for the finite element method: %%%
%
% options.mass_matrix = 1;
% options.bdry_conds = {'DD' 'NN' 'P'};
% options.numeigs = 30;
% [L,M,lambda,V] = laplacian2([100,45,55],options);
% % 3D negative Laplacian with mixed boundary conditions.
% options.mass_matrix = 1;
% options.bdry_conds = {'DD' 'DN'};
% options.numeigs = 30;
% [L,M,lambda,V] = laplacian2([200 200],options);
% % 2D negative Laplacian with mixed boundary conditions.
% options.mass_matrix = 1;
% options.bdry_conds = {'DN'};
% options.num_eigs = 30;
% [L,M,lambda,V] = laplacian2(200,options);
% % 1D negative Laplacian with mixed boundary conditions.
%
% % Example to test if outputs correct eigenvalues and vectors:
% options.mass_matrix = 1;
% options.bdry_conds = {'DD' 'DN' 'P'};
% options.numeigs = 31;
% [L,M,lambda,V] = laplacian2([13,10,6],options);
% [Veig D] = eig(full(L),full(M)); lambdaeig = diag(D(1:31,1:31));
% max(abs(lambda-lambdaeig)) %checking eigenvalues
% subspace(V,Veig(:,1:31)) %checking the invariant subspace
% subspace(V(:,1),Veig(:,1)) %checking selected eigenvectors
% subspace(V(:,30:31),Veig(:,30:31)) %a multiple eigenvalue
%
%
% Class support for inputs:
% N - row vector float double
% options.bdryconds - cell array
% options.numeigs - scalar float double
% options.mass_matrix - boolean
%
% Class support for outputs:
% A,L,M - sparse float double,
% lambda and V - full float double.
%
% Note: the actual numerical entries of A fit int8 format, but only
% double data class is currently (2010) supported for sparse matrices.
%
% This program is designed to efficiently compute the (1-3)D negative
% Laplacian on a rectangular grid for Dirichlet, Neumann, and Periodic
% boundary conditions using sums of tensor products of 1D Laplacians.
% For more information on tensor products, see
% http://en.wikipedia.org/wiki/Kronecker_sum_of_discrete_Laplacians
% For 2D case in MATLAB, see
% http://www.mathworks.com/access/helpdesk/help/techdoc/ref/kron.html.
%
% This code is a part of the BLOPEX package:
% http://en.wikipedia.org/wiki/BLOPEX or directly
% http://code.google.com/p/blopex/
% License: BSD
% Copyright 2012 Bryan C. Smith, Andrew V. Knyazev
% $Revision: 2.0 $ $Date: 15-May-2012
% Tested in MATLAB 7.13.0.564 (R2011b)
tic
% Input handling.
if nargin > 2
error('BLOPEX:laplacian:TooManyInputs',...
'%s','Too many input arguments.');
elseif nargin == 0
error('BLOPEX:laplacian:NoInputArguments',...
'%s','Must have at least one input argument.');
end
if nargout > 4
error('BLOPEX:laplacian:TooManyOutputs',...
'%s','Maximum number of outputs is 4.');
end
if nargin == 2
options = varargin{1};
else
options = struct;
end
u = numpts;
dim2 = size(u);
if dim2(1) ~= 1
error('BLOPEX:laplacian:WrongVectorOfGridPoints',...
'%s','Number of grid points must be in a row vector.')
end
if dim2(2) > 3
error('BLOPEX:laplacian:WrongNumberOfGridPoints',...
'%s','Number of grid points must be a 1, 2, or 3')
end
dim=dim2(2); clear dim2;
uint = round(u);
if max(uint~=u)
warning('BLOPEX:laplacian:NonIntegerGridSize',...
'%s','Grid sizes must be integers. Rounding...');
u = uint; clear uint
end
if max(u<=0 )
error('BLOPEX:laplacian:NonIntegerGridSize',...
'%s','Grid sizes must be positive.');
end
if isfield(options,'numeigs')
m = options.numeigs;
else
m = 0;
end
if isfield(options,'bdry_conds')
B = options.bdry_conds;
else
if dim == 1
B = {'DD'};
elseif dim == 2
B = {'DD' 'DD'};
else
B = {'DD' 'DD' 'DD'};
end
end
if isfield(options,'mass_matrix')
mass_matrix = options.mass_matrix;
else
mass_matrix = false;
end
clear options
if max(size(m) - [1 1]) ~= 0
error('BLOPEX:laplacian:WrongNumberOfEigenvalues',...
'%s','The requested number of eigenvalues must be a scalar.');
end
maxeigs = prod(u);
mint = round(m);
if mint ~= m || mint > maxeigs
error('BLOPEX:laplacian:InvalidNumberOfEigs',...
'%s','Number of eigenvalues output must be a nonnegative ',...
'integer no bigger than number of grid points.');
end
m = mint;
bdryerr = 0;
a = whos('regep','B');
if sum(a.class(1:4)=='cell') ~= 4 || sum(a.size == [1 dim]) ~= 2
bdryerr = 1;
else
BB = zeros(1, 2*dim);
for i = 1:dim
if (length(B{i}) == 1)
if B{i} == 'P'
BB(i) = 3;
BB(i + dim) = 3;
else
bdryerr = 1;
end
elseif (length(B{i}) == 2)
if B{i}(1) == 'D'
BB(i) = 1;
elseif B{i}(1) == 'N'
BB(i) = 2;
else
bdryerr = 1;
end
if B{i}(2) == 'D'
BB(i + dim) = 1;
elseif B{i}(2) == 'N'
BB(i + dim) = 2;
else
bdryerr = 1;
end
else
bdryerr = 1;
end
end
end
if bdryerr == 1
error('BLOPEX:laplacian:InvalidBdryConds',...
'%s','Boundary conditions must be a cell of length 3 for 3D, 2',...
' for 2D, 1 for 1D, with values ''DD'', ''DN'', ''ND'', ''NN''',...
', or ''P''.');
end
% Set the component matrices. SPDIAGS converts int8 into double anyway.
w = [-1 2 -1];
if mass_matrix
w2 = [1 4 1]/6;
end
e1 = ones(u(1),1); %e1 = ones(u(1),1,'int8');
horder = 1;
D1x = spdiags(kron(w,e1), -horder:horder, u(1),u(1));
if mass_matrix
M1x = spdiags(kron(w2,e1), -horder:horder, u(1),u(1));
end
if dim > 1
e2 = ones(u(2),1);
D1y = spdiags(kron(w,e2), -horder:horder, u(2),u(2));
if mass_matrix
M1y = spdiags(kron(w2,e2), -horder:horder, u(2),u(2));
end
end
if dim > 2
e3 = ones(u(3),1);
D1z = spdiags(kron(w,e3), -horder:horder, u(3),u(3));
if mass_matrix
M1z = spdiags(kron(w2,e3), -horder:horder, u(3),u(3));
end
end
% Set boundary conditions if other than Dirichlet.
for i = 1:dim
if BB(i) == 2
eval(['D1' char(119 + i) '(1,1) = 1;'])
if mass_matrix
eval(['M1' char(119 + i) '(1,1) = 5/6;'])
end
elseif BB(i) == 3
eval(['D1' char(119 + i) '(1,' num2str(u(i)) ') = D1'...
char(119 + i) '(1,' num2str(u(i)) ') -1;']);
eval(['D1' char(119 + i) '(' num2str(u(i)) ',1) = D1'...
char(119 + i) '(' num2str(u(i)) ',1) -1;']);
if mass_matrix
eval(['M1' char(119 + i) '(1,' num2str(u(i)) ') = M1'...
char(119 + i) '(1,' num2str(u(i)) ') + 1/6;']);
eval(['M1' char(119 + i) '(' num2str(u(i)) ',1) = M1'...
char(119 + i) '(' num2str(u(i)) ',1) + 1/6;']);
end
end
if BB(i+dim) == 2
eval(['D1' char(119 + i)...
'(',num2str(u(i)),',',num2str(u(i)),') = 1;'])
if mass_matrix
eval(['M1' char(119 + i)...
'(',num2str(u(i)),',',num2str(u(i)),') = 5/6;'])
end
end
end
% Form A using tensor products of lower dimensional Laplacians
Ix = speye(u(1));
if dim == 1
A = D1x;
if mass_matrix
M = M1x;
end
elseif dim == 2
Iy = speye(u(2));
if mass_matrix
M = kron(M1y,M1x);
A = kron(Iy,D1x) + kron(D1y,Ix) - (1/3)*kron(D1y,D1x);
else
A = kron(Iy,D1x) + kron(D1y,Ix);
end
elseif dim == 3
Iy = speye(u(2));
Iz = speye(u(3));
if mass_matrix
M = kron(M1z, kron(M1y,M1x));
A = kron(Iz, kron(Iy, D1x)) + kron(Iz, kron(D1y, Ix))...
+ kron(kron(D1z,Iy),Ix) - (1/3)*(kron(Iz, kron(D1y, D1x))...
+ kron(D1z, kron(Iy, D1x)) + kron(kron(D1z,D1y),Ix))...
+ (1/12)*kron(D1z,kron(D1y, D1x));
else
A = kron(Iz, kron(Iy, D1x)) + kron(Iz, kron(D1y, Ix))...
+ kron(kron(D1z,Iy),Ix);
end
end
% Calculate m smallest exact eigenvalues.
if m > 0
if (BB(1) == 1) && (BB(1+dim) == 1)
a1 = pi/2/(u(1)+1);
N = (1:u(1))';
elseif (BB(1) == 2) && (BB(1+dim) == 2)
a1 = pi/2/u(1);
N = (0:(u(1)-1))';
elseif ((BB(1) == 1) && (BB(1+dim) == 2)) || ((BB(1) == 2)...
&& (BB(1+dim) == 1))
a1 = pi/4/(u(1)+0.5);
N = 2*(1:u(1))' - 1;
else
a1 = pi/u(1);
N = floor((1:u(1))/2)';
end
if mass_matrix
sigma1 = 4*sin(a1*N).^2;
mu1 = (1+2*cos(a1*N).^2)/3;
else
lambda1 = 4*sin(a1*N).^2;
end
clear a1
if dim > 1
if (BB(2) == 1) && (BB(2+dim) == 1)
a2 = pi/2/(u(2)+1);
N = (1:u(2))';
elseif (BB(2) == 2) && (BB(2+dim) == 2)
a2 = pi/2/u(2);
N = (0:(u(2)-1))';
elseif ((BB(2) == 1) && (BB(2+dim) == 2)) || ((BB(2) == 2)...
&& (BB(2+dim) == 1))
a2 = pi/4/(u(2)+0.5);
N = 2*(1:u(2))' - 1;
else
a2 = pi/u(2);
N = floor((1:u(2))/2)';
end
if mass_matrix
sigma2 = 4*sin(a2*N).^2;
mu2 = (1+2*cos(a2*N).^2)/3;
else
lambda2 = 4*sin(a2*N).^2;
end
clear a2
else
lambda2 = 0;
sigma2 = 0;
mu2 = 0;
end
if dim > 2
if (BB(3) == 1) && (BB(6) == 1)
a3 = pi/2/(u(3)+1);
N = (1:u(3))';
elseif (BB(3) == 2) && (BB(6) == 2)
a3 = pi/2/u(3);
N = (0:(u(3)-1))';
elseif ((BB(3) == 1) && (BB(6) == 2)) || ((BB(3) == 2)...
&& (BB(6) == 1))
a3 = pi/4/(u(3)+0.5);
N = 2*(1:u(3))' - 1;
else
a3 = pi/u(3);
N = floor((1:u(3))/2)';
end
if mass_matrix
sigma3 = 4*sin(a3*N).^2;
mu3 = (1+2*cos(a3*N).^2)/3;
else
lambda3 = 4*sin(a3*N).^2;
end
clear a3
else
lambda3 = 0;
sigma3 = 0;
mu3 = 0;
end
if dim == 1
if mass_matrix
lambda = sigma1./mu1;
else
lambda = lambda1;
end
elseif dim == 2
if mass_matrix
lambda = (kron(e2,sigma1) + kron(sigma2,e1) - (kron(sigma2,sigma1))/3)./(kron(mu2,mu1));
size(lambda)
else
lambda = kron(e2,lambda1) + kron(lambda2, e1);
end
clear lambda2 mu2 sigma2
else
if mass_matrix
lambda = ( kron(e3,kron(e2,sigma1)) + kron(e3,kron(sigma2,e1))...
+ kron(sigma3,kron(e2,e1)) - (1/3)*(kron(e3,kron(sigma2,sigma1))...
+ kron(sigma3,kron(e2,sigma1)) + kron(sigma3,kron(sigma2,e1)))...
+ (1/12)*kron(sigma3,kron(sigma2,sigma1)) )./kron(mu3,kron(mu2,mu1));
else
lambda = kron(e3,kron(e2,lambda1)) + kron(e3,kron(lambda2,e1))...
+ kron(lambda3,kron(e2,e1));
end
clear lambda2 lambda3 mu2 mu3 sigma2 sigma3
end
clear lambda1 sigma1 mu1
[lambda, p] = sort(lambda);
if m < maxeigs - 0.1
w = lambda(m+1);
else
w = inf;
end
lambda = lambda(1:m);
p = p(1:m)';
else
lambda = [];
end
% Calculate eigenvectors if specified in output.
if nargout == 3 + mass_matrix && m == 0
V = [];
elseif nargout == 3 + mass_matrix
p1 = mod(p-1,u(1))+1;
if (BB(1) == 1) && (BB(1+dim) == 1)
V1 = sin(kron((1:u(1))'*(pi/(u(1)+1)),p1))*(2/(u(1)+1))^0.5;
elseif (BB(1) == 2) && (BB(1+dim) == 2)
V1 = cos(kron((0.5:1:u(1)-0.5)'*(pi/u(1)),p1-1))*(2/u(1))^0.5;
V1(:,p1==1) = 1/u(1)^0.5;
elseif ((BB(1) == 1) && (BB(1+dim) == 2))
V1 = sin(kron((1:u(1))'*(pi/2/(u(1)+0.5)),2*p1 - 1))...
*(2/(u(1)+0.5))^0.5;
elseif ((BB(1) == 2) && (BB(1+dim) == 1))
V1 = cos(kron((0.5:1:u(1)-0.5)'*(pi/2/(u(1)+0.5)),2*p1 - 1))...
*(2/(u(1)+0.5))^0.5;
else
V1 = zeros(u(1),m);
a = (0.5:1:u(1)-0.5)';
V1(:,mod(p1,2)==1) = cos(a*(pi/u(1)*(p1(mod(p1,2)==1)-1)))...
*(2/u(1))^0.5;
pp = p1(mod(p1,2)==0);
if ~isempty(pp)
V1(:,mod(p1,2)==0) = sin(a*(pi/u(1)*p1(mod(p1,2)==0)))...
*(2/u(1))^0.5;
end
V1(:,p1==1) = 1/u(1)^0.5;
if mod(u(1),2) == 0
V1(:,p1==u(1)) = V1(:,p1==u(1))/2^0.5;
end
end
if dim > 1
p2 = mod(p-p1,u(1)*u(2));
p3 = (p - p2 - p1)/(u(1)*u(2)) + 1;
p2 = p2/u(1) + 1;
if (BB(2) == 1) && (BB(2+dim) == 1)
V2 = sin(kron((1:u(2))'*(pi/(u(2)+1)),p2))*(2/(u(2)+1))^0.5;
elseif (BB(2) == 2) && (BB(2+dim) == 2)
V2 = cos(kron((0.5:1:u(2)-0.5)'*(pi/u(2)),p2-1))*(2/u(2))^0.5;
V2(:,p2==1) = 1/u(2)^0.5;
elseif ((BB(2) == 1) && (BB(2+dim) == 2))
V2 = sin(kron((1:u(2))'*(pi/2/(u(2)+0.5)),2*p2 - 1))...
*(2/(u(2)+0.5))^0.5;
elseif ((BB(2) == 2) && (BB(2+dim) == 1))
V2 = cos(kron((0.5:1:u(2)-0.5)'*(pi/2/(u(2)+0.5)),2*p2 - 1))...
*(2/(u(2)+0.5))^0.5;
else
V2 = zeros(u(2),m);
a = (0.5:1:u(2)-0.5)';
V2(:,mod(p2,2)==1) = cos(a*(pi/u(2)*(p2(mod(p2,2)==1)-1)))...
*(2/u(2))^0.5;
pp = p2(mod(p2,2)==0);
if ~isempty(pp)
V2(:,mod(p2,2)==0) = sin(a*(pi/u(2)*p2(mod(p2,2)==0)))...
*(2/u(2))^0.5;
end
V2(:,p2==1) = 1/u(2)^0.5;
if mod(u(2),2) == 0
V2(:,p2==u(2)) = V2(:,p2==u(2))/2^0.5;
end
end
else
V2 = ones(1,m);
end
if dim > 2
if (BB(3) == 1) && (BB(3+dim) == 1)
V3 = sin(kron((1:u(3))'*(pi/(u(3)+1)),p3))*(2/(u(3)+1))^0.5;
elseif (BB(3) == 2) && (BB(3+dim) == 2)
V3 = cos(kron((0.5:1:u(3)-0.5)'*(pi/u(3)),p3-1))*(2/u(3))^0.5;
V3(:,p3==1) = 1/u(3)^0.5;
elseif ((BB(3) == 1) && (BB(3+dim) == 2))
V3 = sin(kron((1:u(3))'*(pi/2/(u(3)+0.5)),2*p3 - 1))...
*(2/(u(3)+0.5))^0.5;
elseif ((BB(3) == 2) && (BB(3+dim) == 1))
V3 = cos(kron((0.5:1:u(3)-0.5)'*(pi/2/(u(3)+0.5)),2*p3 - 1))...
*(2/(u(3)+0.5))^0.5;
else
V3 = zeros(u(3),m);
a = (0.5:1:u(3)-0.5)';
V3(:,mod(p3,2)==1) = cos(a*(pi/u(3)*(p3(mod(p3,2)==1)-1)))...
*(2/u(3))^0.5;
pp = p1(mod(p3,2)==0);
if ~isempty(pp)
V3(:,mod(p3,2)==0) = sin(a*(pi/u(3)*p3(mod(p3,2)==0)))...
*(2/u(3))^0.5;
end
V3(:,p3==1) = 1/u(3)^0.5;
if mod(u(3),2) == 0
V3(:,p3==u(3)) = V3(:,p3==u(3))/2^0.5;
end
end
else
V3 = ones(1,m);
end
if dim == 1
V = V1;
elseif dim == 2
V = kron(e2,V1).*kron(V2,e1);
else
V = kron(e3, kron(e2, V1)).*kron(e3, kron(V2, e1))...
.*kron(kron(V3,e2),e1);
end
end
if m ~= 0
if abs(lambda(m) - w) < maxeigs*eps('double')
sprintf('\n%s','Warning: (m+1)th eigenvalue is nearly equal',...
' to mth.')
end
end
disp(' ')
toc
a = whos('regexp','A');
if mass_matrix
disp(['The stiffness matrix takes ' num2str(a.bytes) ' bytes'])
a = whos('regep','M');
disp(['The mass matrix takes ' num2str(a.bytes) ' bytes'])
else
disp(['The Laplacian matrix takes ' num2str(a.bytes) ' bytes'])
end
if nargout == 3 + mass_matrix
a = whos('regep','V');
disp(['The eigenvectors take ' num2str(a.bytes) ' bytes'])
end
disp(' ')
k = 1;
varargout{1} = A;
if mass_matrix
k = k + 1;
varargout{k} = M;
end
k = k + 1;
varargout{k} = lambda;
if nargout == 3 + mass_matrix
k = k + 1;
varargout{k} = V;
end