LAPLACIAN Sparse Negative Laplacian in 1D, 2D, or 3D
[~,~,A]=LAPLACIAN(N) generates a sparse negative 3D Laplacian matrix
with Dirichlet boundary conditions, from a rectangular cuboid regular
grid with j x k x l interior grid points if N = [j k l], using the
standard 7-point finite-difference scheme, The grid size is always
one in all directions.
[~,~,A]=LAPLACIAN(N,B) specifies boundary conditions with a cell array
B. For example, B = {'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 B are 'DD', 'DN', 'ND', 'NN' and 'P'.
LAMBDA = LAPLACIAN(N,B,M) or LAPLACIAN(N,M) outputs the m smallest
eigenvalues of the matrix, computed by an exact known formula, see
http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_of_the_second_derivative
It will produce a warning if the mth eigenvalue is equal to the
(m+1)th eigenvalue. If m is absebt or zero, lambda will be empty.
[LAMBDA,V] = LAPLACIAN(N,B,M) also outputs orthonormal eigenvectors
associated with the corresponding m smallest eigenvalues.
[LAMBDA,V,A] = LAPLACIAN(N,B,M) produces a 2D or 1D negative
Laplacian matrix if the length of N and B are 2 or 1 respectively.
It uses the standard 5-point scheme for 2D, and 3-point scheme for 1D.
% Examples:
[lambda,V,A] = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
% Everything for 3D negative Laplacian with mixed boundary conditions.
laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
% or
lambda = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
% computes the eigenvalues only
[~,V,~] = laplacian([200 200],{'DD' 'DN'},30);
% Eigenvectors of 2D negative Laplacian with mixed boundary conditions.
[~,~,A] = laplacian(200,{'DN'},30);
% 1D negative Laplacian matrix A with mixed boundary conditions.
% Example to test if outputs correct eigenvalues and vectors:
[lambda,V,A] = laplacian([13,10,6],{'DD' 'DN' 'P'},30);
[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] = laplacian([30,30]);
norm(A1-A2,inf)
Class support for inputs:
N - row vector float double
B - cell array
M - scalar float double
Class support for outputs:
lambda and V - full float double, A - sparse 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 eigenvalues,
eigenvectors, and the sparse matrix of the (1-3)D negative Laplacian
on a rectangular grid for Dirichlet, Neumann, and Periodic boundary
conditions using tensor sums 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 also part of the BLOPEX package:
http://en.wikipedia.org/wiki/BLOPEX or directly
http://code.google.com/p/blopex/
Package: specfun
