Sparse (13)D Laplacian on a rectangular grid with exact eigenpairs.
The code computes the exact eigenpairs of (13)D negative Laplacian on a rectangular finitedifference grid for combinations of Dirichlet, Neumann, and Periodic boundary conditions using explicit formulas from
http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_of_the_second_derivative
The code can also compute the sparse matrix itself, using Kronecker sums of 1D Laplacians. For more information on tensor sums, see
http://en.wikipedia.org/wiki/Kronecker_sum_of_discrete_Laplacians
Example, compute everything for 3D negative Laplacian with mixed boundary conditions:
[lambda,V,A] = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
Compute only the eigenvalues:
lambda = laplacian([100,45,55],{'DD' 'NN' 'P'}, 20);
Compute the matrix only:
[~,~,A] = laplacian([100,45,55],{'DD' 'NN' 'P'});
GNU OCTAVE compatible.
This code is a part of the BLOPEX eigensolver package, see
http://en.wikipedia.org/wiki/BLOPEX
or go directly to
http://code.google.com/p/blopex/
Copyright owners: Bryan C. Smith and Andrew V. Knyazev
1.2  added a conversion to a toolbox 

1.2  Revision 1.1 changes: rearranged the output variables, always compute the eigenvalues, compute eigenvectors and/or the matrix on demand only. 

1.1  updated description 
Andrew Knyazev (view profile)
Hi Varun Shankar, I am not familiar with the "ghost point based implementation on a vertexcentered grid". Our implementation of the Neumann BCs in 1D gives the first raw [1 1 0 ... 0] to my recollection. 2D and 3D cases are computed just as products of 1D.
Varun Shankar (view profile)
Hi Andrew,
I have a quick question. Does the implementation of Neumann BCs herein end up being equivalent to the standard ghost point based implementation on a vertexcentered grid? It certainly looks like that in 1D, but I couldn't tell for certain in 2D.
Thanks!
Andrew Knyazev (view profile)
Dear Ales,
Not hearing back from you, I have performed a few tests myslef, in 1D, 2D and 3D. Please the results quoted below. I see no problems with the code, and close this bug case.
Thanks, Andrew
+++++++++++++++++++++++++++++++++++
octave:52> [lambda,V,A] = laplacian([2,4,5],{'P' 'P' 'P'}, 20);
ans =
Warning: (m+1)th eigenvalue is nearly equal
to mth.
The eigenvectors take 6400 bytes
The Laplacian matrix takes 3524 bytes
octave:53> whos
Variables in the current scope:
Attr Name Size Bytes Class
==== ==== ==== ===== =====
A 40x40 3524 double
V 40x20 6400 double
lambda 20x1 160 double
Total is 1100 elements using 10084 bytes
octave:54> [lambda,V,A] = laplacian([2,4],{'P' 'P' }, 2);
ans =
Warning: (m+1)th eigenvalue is nearly equal
to mth.
The eigenvectors take 128 bytes
The Laplacian matrix takes 516 bytes
octave:55> whos
Variables in the current scope:
Attr Name Size Bytes Class
==== ==== ==== ===== =====
A 8x8 516 double
V 8x2 128 double
lambda 2x1 16 double
Total is 58 elements using 660 bytes
octave:56> [lambda,V,A] = laplacian([2],{'P'}, 2);
The eigenvectors take 32 bytes
The Laplacian matrix takes 60 bytes
octave:57> whos
Variables in the current scope:
Attr Name Size Bytes Class
==== ==== ==== ===== =====
A 2x2 60 double
V 2x2 32 double
lambda 2x1 16 double
Total is 10 elements using 108 bytes
Andrew Knyazev (view profile)
Re: no solution when periodic BC are applied.
Dear Ales,
Could you please provide a specific complete example, quoting the function call, the full output, and MATLAB version? I cannot fix a problem without being able to reproduce the problem myslef.
Thanks, Andrew
Ales (view profile)
Dear Andrew:
Thank you for a nice software. I ran into a problem  there is no solution when periodic BC are applied. I think it is endemic to these types of problems. Do you know why?
Marios Karaoulis (view profile)
Thanks for the reply Andrew, you are right.
Andrew Knyazev (view profile)
Re: Marios Karaoulis
It depends on boundary conditions. Your example apparently uses Neumann boundary condition as is http://en.wikipedia.org/wiki/Laplacian_matrix#As_an_approximation_to_the_negative_continuous_Laplacian
while the code default us the Dirichlet boundary conditions. To get your matrix, please use
[~,~,A]=laplacian([3 3],{'NN' 'NN'})
full(A) shows what you want:
2 1 0 1 0 0 0 0 0
1 3 1 0 1 0 0 0 0
0 1 2 0 0 1 0 0 0
1 0 0 3 1 0 1 0 0
0 1 0 1 4 1 0 1 0
0 0 1 0 1 3 0 0 1
0 0 0 1 0 0 2 1 0
0 0 0 0 1 0 1 3 1
0 0 0 0 0 1 0 1 2
Marios Karaoulis (view profile)
Hi, useful add on, but one question.
Let's say we have a 3x3 grid, then why all the element on the diagonal are 4? Should be something like this
[2 1 0 1 0 0 0 0 0
1 3 1 0 1 0 0 0
0 1 2 0 0 1 0 0 0 ....
...]
Only line 5 should have 4 on it's diagonal.
Mr Smart (view profile)
Thomas Clark (view profile)
This works extremely well, extremely flexibly and with no fuss.
I'm SO upset I didn't find this earlier. I've basically coded an alternative, already, but it isn't close to being as good.
Thank you very much!