laplacian, a FORTRAN90 code 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^2for a function of two dimensions, we have
L(U)(X,Y) = d^2 U(X,Y)/dX^2 + d^2 U(X,Y)/dY^2and 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^2The 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^2and 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 termsand 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^2This 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 * Uwhere 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 -1Square 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:
The computer code and data files made available on this web page are distributed under the GNU LGPL license.
laplacian is available in a C version and a C++ version and a FORTRAN90 version and a MATLAB version.
TEST_MAT, a FORTRAN90 code 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.