Code covered by the BSD License  

Highlights from
Laplacian in 1D, 2D, or 3D

5.0

5.0 | 4 ratings Rate this file 84 Downloads (last 30 days) File Size: 5.27 KB File ID: #27279
image thumbnail

Laplacian in 1D, 2D, or 3D

by

 

17 Apr 2010 (Updated )

Sparse (1-3)D Laplacian on a rectangular grid with exact eigenpairs.

| Watch this File

File Information
Description

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

MATLAB release MATLAB 7.11 (R2010b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (8)
20 May 2014 Andrew Knyazev

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

17 May 2014 Andrew Knyazev

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

16 May 2014 Ales

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?

14 Mar 2012 Marios Karaoulis

Thanks for the reply Andrew, you are right.

12 Mar 2012 Andrew Knyazev

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

11 Mar 2012 Marios Karaoulis

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.

12 Sep 2011 Mr Smart  
07 Nov 2010 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
30 Apr 2010

updated description

01 Aug 2011

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

Contact us