LAPLACIAN
The Discrete Laplacian Operator

LAPLACIAN is a FORTRAN77 library which carries out computations related to the discrete Laplacian operator, including full or sparse evaluation, evaluation for unequally spaced data sampling points, application to a set of data samples, solution of associated linear systems, eigenvalues and eigenvectors, and extension to 2D and 3D geometry.

For a (twice-continuously differentiable) function U(X) of a one-dimensional variable X, the (continuous) Laplacian operator L(U)(X) is simply the second derivative:

        L(U)(X) = d^2 U(X)/dX^2

for a function of two dimensions, we have

        L(U)(X,Y) = d^2 U(X,Y)/dX^2 + d^2 U(X,Y)/dY^2

and so on. The Laplacian is important mathematically because it arises naturally in the mathematical model of many important physical systems, such as the variation of temperature in a heated object.

The discrete Laplacian is an approximation to the continuous Laplacian that is appropriate when data is known or sampled only at finitely many points. It is often the case that these points are evenly spaced in a line or grid.

The one-dimensional discrete Laplacian Suppose that a function U(X) is known at three points X-H, X and X+H. Then the discrete Laplacian operator, applied to this data, is simply the standard approximation to the second derivative:

        L(U)(X) = ( + 2 U(X) 
                    -   U(X-H)
                    -   U(X+H) ) / H^2

The 2-dimensional analog, assuming that (X,Y) data is available at the uniform spacing H, is:

        L(U)(X,Y) = ( + 4 U(X,Y)
                      -   U(X-H,Y-H) - U(X-H,Y+H)
                      -   U(X+H,Y-H) - U(X+H,Y+H) ) / H^2

and similar results apply for higher dimensions.

If the data is not available at equally spaced points, then the computation becomes somewhat more involved. However, simply using Taylor series, we have:

        U(X+H1) = U(X) + U' * H1 + U'' * H1^2 / 2 + U''' H1^3 / 6 + ...
        U(X+H2) = U(X) + U' * H2 + U'' * H2^2 / 2 + U''' H2^3 / 6 + ...

from which we can determine that:

        H2 * U(X+H1) - H1 * U(X+H2) - ( H2 - H1 ) * U(X) =
          ( H2 * H1^2 / 2 - H1 * H2^2 / 2 ) * U''
        + ( H2 * H1^3 / 6 - H1 * H3^3 / 6 ) * U'''
        + higher order terms

and hence, we have the approximation:

        L(U)(X) = ( H2 * U(X+H1) - H1 * U(X+H2) - ( H2 - H1 ) * U(X) )
                / ( H2 * H1^2 / 2 - H1 * H2^2 / 2 )

Again, corresponding results can be determined for higher dimensions, in cases where the data is sampled along coordinate lines, but with nonuniform spacing.

Now, let us consider the 1-dimensional case, where U(X) is available at N+2 points equally spaced by H, and indexed from 0 to N+1. Let U now be the vector of values U(X(0)), U(X(2)), ..., U(X(N+1)), and let U(i) indicate the function value U(X(i)). For all but the first and last indices, it is easy to estimate the second derivative, using the formula:

        L(Ui) = U''(i) = ( - U(i-1) + 2 U(i) - U(i+1) ) / h^2

This is a linear operation. It takes N+2 values U and produces N values L(Ui), for i = 1 to N. It has the matrix form:

        L(Ui) = L * U

where the N by N+2 matrix L has the form (if N = 4, and H = 1 ):

         -1  2 -1  0  0  0
          0 -1  2 -1  0  0
          0  0 -1  2 -1  0
          0  0  0 -1  2 -1

Square matrices are much more convenient for analysis, so let us assume either that U(1) and U(N) are zero, or that we can somehow neglect or defer the analysis of the first and last columns. In that case, we have our first example of a discrete Laplacian matrix, in this case for N = 4 and H = 1:

           2 -1  0  0 
          -1  2 -1  0 
           0 -1  2 -1 
           0  0 -1  2

The purpose of this library is to study issues related to matrices of this form, including:

producing a copy of the matrix in full storage mode;
producing a copy of the matrix in sparse storage mode;
handling cases where U(0) or U(N+1) are not given, but instead derivative information or periodic boundary conditions are applied;
determining the result of applying the discrete Laplacian to data; that is, multiplying the matrix times a vector of data;
solving a linear system involving the discrete Laplacian;
determining the eigenvalues and eigenvectors;
handling cases where the spacing is nonuniform;
handling cases in 2D or 3D.

Licensing:

The computer code and data files made available on this web page are distributed under the GNU LGPL license.

Languages:

LAPLACIAN is available in a C version and a C++ version and a FORTRAN77 version and a FORTRAN90 version and a MATLAB version.

Related Data and Programs:

TEST_MAT, a FORTRAN77 library which defines test matrices for which some of the determinant, eigenvalues, inverse, null vectors, P*L*U factorization or linear system solution are already known.

Source Code:

laplacian.f, the source code.
laplacian.sh, BASH commands to compile the source code.

Examples and Tests:

laplacian_prb.f, a sample calling program.
laplacian_prb.sh, BASH commands to compile and run the sample program.
laplacian_prb_output.txt, the output file.

List of Routines:

CHOLESKY_UPPER_ERROR determines the error in an upper Cholesky factor.
EIGEN_ERROR determines the error in a (right) eigensystem.
INVERSE_ERROR determines the error in an inverse matrix.
L1DD_APPLY applies the 1D DD Laplacian to a vector.
L1DD_CHOLESKY computes the Cholesky factor of the 1D DD Laplacian.
L1DD_EIGEN returns eigeninformation for the 1D DD Laplacian.
L1DD_FULL stores the 1D DD Laplacian as a full matrix.
L1DD_FULL_INVERSE stores the inverse of the 1D DD Laplacian.
L1DN_APPLY applies the 1D DN Laplacian to a vector.
L1DN_CHOLESKY computes the Cholesky factor of the 1D DN Laplacian.
L1DN_EIGEN returns eigeninformation for the 1D DN Laplacian.
L1DN_FULL stores the 1D DN Laplacian as a full matrix.
L1DN_FULL_INVERSE stores the inverse of the 1D DN Laplacian.
L1ND_APPLY applies the 1D ND Laplacian to a vector.
L1ND_CHOLESKY computes the Cholesky factor of the 1D ND Laplacian.
L1ND_EIGEN returns eigeninformation for the 1D ND Laplacian.
L1ND_FULL stores the 1D ND Laplacian as a full matrix.
L1ND_FULL_INVERSE stores the inverse of the 1D ND Laplacian.
L1NN_APPLY applies the 1D NN Laplacian to a vector.
L1NN_CHOLESKY computes the Cholesky factor of the 1D NN Laplacian.
L1NN_EIGEN returns eigeninformation for the 1D NN Laplacian.
L1NN_FULL stores the 1D NN Laplacian as a full matrix.
L1PP_APPLY applies the 1D PP Laplacian to a vector.
L1PP_CHOLESKY computes the Cholesky factor of the 1D PP Laplacian.
L1PP_EIGEN returns eigeninformation for the 1D PP Laplacian.
L1PP_FULL stores the 1D PP Laplacian as a full matrix.
LAPLACIAN_1D_UNEVEN_APPLY applies the 1D Discrete Laplacian to a vector.
R8MAT_PRINT prints an R8MAT.
R8MAT_PRINT_SOME prints some of an R8MAT.
R8VEC_PRINT prints an R8VEC.
TIMESTAMP prints the current YMDHMS date as a time stamp.

You can go up one level to the FORTRAN77 source codes.

Last revised on 29 October 2013.