laplacian_matrix, a Fortran77 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. The approximation error is O(h^2).
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:
The information on this web page is distributed under the MIT license.
laplacian_matrix is available in a C version and a C++ version and a Fortran77 version and a Fortran90 version and a MATLAB version and an Octave version.
test_matrix, a Fortran77 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.