File Exchange

image thumbnail

Laplacian in 1D, 2D, or 3D

version 1.2.0.0 (5.27 KB) by Andrew Knyazev
Sparse (1-3)D Laplacian on a rectangular grid with exact eigenpairs.

10 Downloads

Updated 15 May 2015

View License

The code computes the exact eigenpairs of (1-3)D negative Laplacian on a rectangular finite-difference 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

Comments and Ratings (10)

Hi Varun Shankar, I am not familiar with the "ghost point based implementation on a vertex-centered 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

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 vertex-centered grid? It certainly looks like that in 1D, but I couldn't tell for certain in 2D.

Thanks!

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

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

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?

Thanks for the reply Andrew, you are right.

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

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

Thomas Clark

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!

Updates

1.2.0.0

added a conversion to a toolbox

1.2.0.0

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

1.1.0.0

updated description

MATLAB Release Compatibility
Created with R2010b
Compatible with any release
Platform Compatibility
Windows macOS Linux

Discover Live Editor

Create scripts with code, output, and formatted text in a single executable document.


Learn About Live Editor