File Exchange

## Laplacian in 1D, 2D, or 3D

version 1.2 (5.27 KB) by

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

Updated

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

Copyright owners: Bryan C. Smith and Andrew V. Knyazev

Andrew Knyazev

### Andrew Knyazev (view profile)

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

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

Thanks!

Andrew Knyazev

### 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

### 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

### 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

### Marios Karaoulis (view profile)

Thanks for the reply Andrew, you are right.

Andrew Knyazev

### 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

### 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

Thomas Clark

### 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!