## Abstract

Procedures for constructing boundary integral equations equivalent to linear boundary-value problems governed by partial differential equations are well established. In this paper, it is demonstrated how these procedures can be extended to linear boundary-value problems defined on lattices and governed by algebraic (‘difference’) equations. The boundary equations that arise are then themselves algebraic equations. Such ‘boundary algebraic equations’ (BAEs) are derived for fundamental boundary-value problems defined on both perfect lattices and lattices with defects. It is demonstrated that key advantages of representing a continuum boundary-value problem as an equation on the boundary, such as favourable spectral properties and minimal problem size, are preserved in the lattice environment. Certain spectral properties of BAEs are established rigorously, whereas others are supported by numerical experiments.

## 1. Introduction

In this paper, we present solution methods for potential (conduction) problems defined on periodic lattices or structures made of such lattices. Conceptually, these methods are similar to boundary integral equation (BIE) methods of classical potential theory, as the governing difference equation, defined on the entire domain, is replaced by an equivalent algebraic equation defined on the boundary. For continuum potential problems, the boundary equation is an integral equation whose kernel is derived from the fundamental solution of the Laplace operator. For a discrete potential problem, the boundary equation is an algebraic equation whose coefficient matrix is derived from the fundamental solution of the discrete potential problem.

The idea of replacing difference equations defined on a regular lattice with boundary algebraic equations (BAEs) goes back at least to Saltzer (1958). In the current terminology (e.g. Bonnet 1999), Saltzer developed an ‘indirect’ BAE corresponding to the potential problem defined on a square lattice. Saltzer motivated his work by supposing that a small, but dense, linear algebraic problem associated with a BAE would be preferable to a large, but sparse, linear algebraic problem associated with the original difference equation. However, a numerical example presented by Saltzer lent no support to this supposition and therefore the idea was abandoned.

The premise of this paper is that BAEs are worth revisiting in light of recent developments that led to the emergence of fast methods for solving BIEs. These include fast matrix–vector multiplication methods (such as those described by Rokhlin 1985; Barnes & Hut 1986; Greengard & Rokhlin 1987, 1997) that exploit certain properties of the harmonic kernel to compute the matrix–vector product approximately, using arithmetic operations and storage. Here, *N*_{b} is the number of unknowns associated with the boundary discretization. In addition, it has been realized that BIEs give rise to much better conditioned linear algebraic problems than finite-element or other widely used methods based on domain discretization. In particular, for matrices associated with the second-kind BIE involving the double-layer kernel, the condition number is independent of the problem size. Synergetic combinations of fast matrix–vector multiplication methods and favourable spectral properties give rise to fast iterative methods, capable of solving linear algebraic problems associated with BIEs using arithmetic operations and storage. Such methods significantly improve the competitiveness of boundary-element and other numerical methods based on BIEs.

We believe that BAEs are also amenable to fast iterative methods. It is clear that methods for applying the relevant system matrix to a vector can be developed, simply because the continuum and discrete fundamental solutions become indistinguishable at large distances. Moreover, the BAEs derived in this paper appear to be at least as well-conditioned as the corresponding BIEs, making them highly amenable to iterative solvers. We support this claim with extensive numerical experiments and provide rigorous proofs of key spectral properties.

Preliminary versions of the material presented in this paper have appeared in non-archival publications (Martinsson 2002; Martinsson & Rodin 2002*b*). BAEs have also been applied to modelling of defects and interfaces in lattices by Haq *et al*. (2006, 2007). In this regard, let us mention that matching atomistic models to lattices is much easier than to continuum, and that BAEs provide a natural framework for matching.

For demonstration purposes, we focus on heat conduction problems defined on two-dimensional lattices. Extensions to more complex two- and three-dimensional lattices should be transparent. The rest of the paper is organized as follows. In §2, we introduce the notation and define the discrete Dirichlet boundary-value problem. In §3, we introduce the discrete single- and double-layer kernels. In §4, we construct direct and indirect BAEs for the homogeneous Dirichlet boundary-value problem. In §5, we generalize the BAEs to lattices with defects. In §6, we state and prove two theorems that reveal the superior spectral properties of matrices associated with BAEs. In §7, we present numerical examples that provide a better understanding of theorems considered in §6. In §8, we summarize the main results.

## 2. Problem statement

Consider an infinite square lattice in whose nodes occupy locations with integer coordinates . Each node is connected to its four nearest neighbours with bonds of unit conductivity (figure 1). A finite lattice structure is formed from the infinite lattice by cutting certain bonds. The algorithm for forming a finite structure is as follows (figure 2*a*).

Identify a finite set of nodes that define the structure.

Remove all bonds with one node belonging to and the other to . Thus, the structure is formed by bonds with both nodes belonging to .

Identify the set formed by nodes that used to be connected to at least one removed bond. We refer to Γ as the boundary nodes and to as internal nodes.

Note that this algorithm is generic, as it does not rely on the lattice details. Furthermore, it can be easily extended to infinite structures.

The model problem is formulated as a homogeneous Dirichlet problem. That is, the structure contains no heat sources in its interior and it is subjected to certain given nodal temperatures on Γ. The corresponding mathematical statement is
2.1
where
2.2
In equation (2.2), *e*_{1}=[1,0] and *e*_{2}=[0,1]. The stated problem has a unique solution for any *g*.

## 3. Single- and double-layer kernels

The fundamental solution of *A* is well known (e.g. Maradudin *et al.* 1960; Joyce 1994; Glasser & Boersma 2000; Martinsson 2002; Martinsson & Rodin 2002*a*),
3.1
This function satisfies
3.2
where , and approaches the fundamental solution of the Laplace operator for large |*m*|,
3.3
Unlike the fundamental solution of the Laplace operator, **G**(*m*) is finite at the origin; in fact, **G**(0)=0.

The single-layer kernel is defined by
3.4
Following the classical definition from the continuum case, we define the double-layer kernel **D**(*m*,*n*) as the incoming flux at a node *n*∈Γ due to a unit source placed at a node . In the classical definition, the term ‘incoming’ necessitates the definition of an oriented elementary surface at *n*. For lattices, this surface is associated with cross sections of the bonds removed from the node *n*, while forming the structure. Thus, the incoming flux can be computed by summing up the fluxes through the removed bonds,
3.5
where is the set of nodes that used to be connected to the node *n* (figure 2*b,c*).

## 4. Boundary algebraic equations

### a Indirect boundary algebraic equations

An indirect BAE corresponding to equation (2.1) can be constructed using the ansatz
4.1
where ϕ is a yet unspecified function on Γ. As constructed, the function *u* automatically satisfies [*A**u*](*m*)=0 for any *m*∈Ω. Therefore, the governing equation for ϕ is obtained by enforcing the boundary conditions
4.2
The function ϕ has a clear physical interpretation—it represents a set of sources applied to the infinite lattice at Γ, so that the nodal temperatures satisfy equation (2.1) in .

The continuum equivalent of equation (4.2) is a first-kind equation involving a compact integral operator that is positive, but not positive definite (in other words, points in its spectrum get arbitrarily close to zero, but do not reach zero). Upon discretization, a compact positive operator is approximated by a positive definite matrix whose smallest eigenvalue is finite, but which depends on the mesh size. As the mesh size grows, the smallest eigenvalue tends to zero, leading to ill-conditioned equations. Results of numerical experiments presented in §6 indicate that this pattern of behaviour is also exhibited by the matrix associated with the discrete single-layer kernel. In particular, this matrix is positive definite, but its condition number increases with the problem size (for further details, refer to Atkinson 1997).

Alternatively, one can formulate an indirect BAE using the ansatz
4.3
The function *u* defined by the ansatz (4.3) automatically satisfies [*A**u*](*m*)=0 for *m*∈Ω (as did the single-layer potential (4.1)), so the governing equation for ψ has the form
4.4
The function ψ can be given a physical interpretation, but it is not necessary for our purposes.

On the surface, equation (4.4) appears to be a discrete version of a first-kind equation. This is somewhat unexpected because, in the continuum setting, the ansatz corresponding to equation (4.3) leads to an equation of the second kind. Nevertheless, results of numerical experiments indicate that the spectral properties of equation (4.4) are similar to those exhibited by equations of the second kind. In particular, for many geometries, the condition number associated with equation (4.4) is independent of the problem size.

### b Direct boundary algebraic equations

A direct BAE equivalent to equation (2.1) can be derived using the reciprocity theorem. To this end, we consider two thermal states on . The first state involves the temperature distribution *u* that solves equation (2.1) and the source distribution *f* along Γ that represents the incoming heat flux; at this stage, *f* is unknown. The second thermal state is defined as the restriction to of the thermal state induced in the infinite lattice by a unit source applied at a node *m*∈Γ. Accordingly, for *n*∈Γ, the temperature and incoming flux of the second state are *u*(*n*)=**S**(*n*,*m*) and *f*(*n*)=**D**(*m*,*n*), respectively. The reciprocity theorem applied to the two states yields
4.5
The above equation parallels the direct BIE of potential theory. Note that, for prescribed Dirichlet data, equation (4.5) is similar to equation (4.2). The difference is that equation (4.2) involves sources applied to the infinite lattice, whereas equation (4.5) involves sources applied to the finite structure .

## 5. Defects

The BAE formulations derived in §4 can be extended to include lattices with inclusions and defects such as missing bars. The resulting boundary formulations frequently have as good spectral properties as the ones for the original lattice. Moreover, when the inclusions affect only a small part of the lattice, the advantage stemming from a reduction of problem size will remain.

To illustrate the idea, let us consider the lattice equation
5.1
where *A*′ represents a lattice equilibrium operator in which a small number of entries have been changed compared to the operator *A* representing the perfect lattice. The idea is to rewrite equation (5.1) as an equation involving the unperturbed operator *A*,
5.2
and couple it with an equation for the ‘fictitious’ body load μ,
5.3
The net effect is to replace the original equation (5.1), whose only unknown is *u*, with a system of equations formed by equations (5.2) and (5.3), which now has both *u* and μ as unknowns. The advantage is that the system formed by equations (5.2) and (5.3) is readily handled using the methods of §4.

To demonstrate the procedure, let us consider a lattice in which there are *J* ‘defective’ bars. For simplicity, we assume that there are no defective bars sharing a node, and all defective bars are connecting internal nodes only. Let us denote the conductivity of the *j*th defective bar by *r*_{j} and the bounding nodes by and . Then,
We now set
5.4
and define the function μ via
To construct the BAE, we represent *u* in the form
5.5
We seek to determine {σ(*n*)}_{n∈Γ} and so that equations (5.2) and (5.3) are satisfied. That *u* satisfies the equation *A* *u*=μ follows directly from the definition (5.5). The condition that *u*=*g* on Γ takes the form
5.6
and finally, the equivalent of equation (5.3) is obtained by inserting equation (5.5) into equation (5.4),
5.7
Together, equations (5.6) and (5.7) determine σ and μ.

## 6. Spectral properties of the boundary algebraic equations

In this section, we present two theorems related to the spectral properties of the matrices associated with BAEs. For the single-layer operator (4.2), we establish bounds for both the largest and the smallest singular values. For the double-layer operator (4.4), we establish only an upper bound for the largest singular value. In §7, we discuss how the theoretical results given in this section are compared with those of numerical experiments. It appears that our bounds for the single-layer potential are sharp and consistent with correct asymptotic behaviour, whereas those for the double-layer potential could probably be improved.

For our purposes, it is advantageous to introduce the forward ∂_{j} and the backward difference operators defined as
where, as before, *e*_{1}=[1,0] and *e*_{2}=[0,1]. With this notation, we have
Moreover, if *u* and *v* are two lattice potentials that decay sufficiently fast at infinity, by carrying out straightforward algebraic manipulations, one can establish the result
6.1

The first theorem concerns the matrix associated with the single-layer potential (4.2). It provides an upper bound for its largest singular value and a lower bound for its smallest singular value.

## Theorem 6.1

*Let Ω be a connected domain in whose boundary Γ contains N nodes. Let S be the N×N matrix corresponding to the linear system (4.2). Then, any singular value λ of S satisfies*

## Proof.

The upper bound follows directly from the inequality
which holds for all . As Ω is connected, |*m*−*n*|≤*N* for any *m*,*n*∈Γ; therefore, the magnitude of any entry of the *N*×*N* matrix *S* is bounded by and the magnitude of any singular value is bounded by .

For the lower bound, we begin with noting that *S* is symmetric and therefore it is sufficient to prove that, for any boundary data σ,
6.2
We prove equation (6.2) using an energy argument. Given a boundary flux data σ, let *u* denote the temperature on defined by the formula (4.1). Then,
6.3
Here, the identity (a) is justified by the fact that [*A* *u*](*m*)=0 when *m*∉Γ and (b) follows from equation (6.1). For each *m*∈Γ, the flux equilibrium equation
implies the inequality
6.4
Now, equation (6.2) follows from equation (6.3) with equation (6.4). ▪

## Remark 6.2

Theorem 6.1 implies that the condition number of the single-layer operator cannot grow faster than as problem size *N* grows. Numerical experiments indicate that the condition number in fact does grow as , but that the logarithmic factor appears to be because of a lone outlying singular value and that the ‘effective’ condition number grows as *O*(*N*), see §7 and, in particular, remark 7.1.

The following theorem states that the norm of a matrix associated with the double-layer kernel is typically very moderate. Our proof of this statement is restricted to contours Γ whose nodes can be numbered in such a way that the distance between the *indices* of any two nodes is a lower bound for the Euclidean distance between those nodes in space. This restriction rules out contours characterized by large aspect ratios.

## Theorem 6.3

*Let c be a real number such that 0<c≤1. Suppose that Γ is a boundary of a lattice domain Ω with the following property: there exists an ordering of the nodes in Γ such that*
6.5
*where*
*Then, the matrix D of the double-layer potential (4.3) associated with Ω satisfies*
6.6

## Proof.

One can verify that the double-layer kernel **D** associated with any contour satisfies
6.7
By combining equations (6.5) and (6.7) and the condition 0<*c*≤1, we obtain
6.8
Equation (6.8) implies that the norm of *D* is bounded from above by 1/*c* times the norm of the *N*×*N* matrix *D*^{+} defined by
As *D*^{+} is a convolution matrix, bounds on its operator norm can easily be obtained in Fourier space. For instance,
▪

## Remark 6.4

In the numerical experiments presented in §7, the norm of the double-layer operator appears to be uniformly bounded and not grow with problem size *N*, as indicated by equation (6.6). For the specific geometries considered, this can be explained by replacing equation (6.7) by a sharper estimate on the decay of **D**(*m*,*n*) as |*m*−*n*| grows. To be precise, equations (3.3) and (3.5) imply that
where the vector is defined as a lattice equivalent of the outward pointing normal to Γ at the boundary point *n*,
For lattice domains in which, loosely speaking, the ‘radius of curvature’ of the boundary is much larger than the lattice spacing, the term (*m*−*n*)⋅ν is small in the critical part of the sum (6.8) and no logarithmic factor arises.

## Remark 6.5

In complete analogy with the continuum case, one can establish that the double-layer kernel has an eigenvalue equal to unity and that the corresponding eigenvector is a constant vector.

## 7. Numerical examples

In this section, we investigate the spectral properties of the BAEs constructed in §§4 and 5. We consider six different lattice geometries. For each geometry, we analysed a sequence of domains of increasing size, assembled the matrices associated with the relevant BAEs and computed the largest and the smallest singular values. The six geometries are the following:

*L-shape*. A sequence of self-similar L-shaped domains shown in figure 3*a*.*Ellipse*. A sequence of quasi-elliptical domains shown in figure 3*b*. The nodes forming this domain are inside ellipses with the aspect ratio 0.75.*Rectangle*. A sequence of rectangular domains whose short side was kept fixed at seven units (figure 3*c*). This sequence of domains violates restrictions of theorem 6.3 as the aspect ratio and the problem size increase.*Slit*. A sequence of square domains with a row of vertical bars removed. The length of the row was kept fixed at half of the size of the square domain (figure 3*d*).*Soft slit*. This sequence of domains is the same as the slit sequence, but rather than removing the bars in the slit, they are replaced with bars with conductivity equal to one-half.*Random*. A sequence of square domains in which a number of bars are removed at random (figure 3*e*). The bars are removed so that no node connects to more than one removed bar, and no bar is connected to a boundary node. The number of bars removed equals the number of nodes on one side of the square.

For each geometry, a sequence of domains with the number of boundary nodes *N* between 100 and 800 was analysed. The inclusions in the slit, soft slit and random lattices were handled using the formulation described in §5.

The largest and the smallest singular values of the matrices associated with the single-layer potential are shown in figure 4. The corresponding plot of the smallest and the largest singular values of the double-layer potential are plotted in figure 5.

The numerical results for the single-layer potential indicate that the theoretical results in theorem 6.1 are sharp in the sense that they have the correct asymptotic behaviour. The constants are not quite optimal, however. For instance, it appears that the correct lower bound for the spectrum of the single-layer potential is , rather than 1/8. Note that the theorem applies only to the first three examples considered (the L-shape, ellipse and rectangle) because the geometries with inclusions were modelled using the augmented equations (5.6) and (5.7). In particular, the decay of the smallest singular value for the slit problem does not contradict theorem 6.1.

## Remark 7.1

Although the single-layer potential leads to fairly ill-conditioned equations, the situation is significantly better than figure 4 indicates. It turns out that the matrix *S* has a single outlying singular value that is much larger than the others. This singular value corresponds to a singular vector that is almost, but not exactly, constant. Recalling that (cf. equation (6.3))
it seems likely that this term arises from the first term in the asymptotic expansion of the potential at infinity,
This hypothesis suggests that when an iterative solver is used, it may be better to work with the operator obtained by projecting away from the constant functions
7.1
where *e* is a vector of unit length whose entries are all identical and handles the constant vector separately. As expected, the spectral properties of *S*′ turn out to be much better than those of *S*, as illustrated in figure 6. In fact, numerical experiments indicate that the quotient between the largest and the smallest non-zero singular value of *S*′ scales as *O*(*N*), rather than the scaling of *S*.

Considering next the numerical results for the double-layer operator, the experiments indicate that, for the geometries considered, the top singular value is uniformly bounded and does not grow logarithmically with problem size, as indicated by theorem 6.3 (cf. remark 6.4). Moreover, it appears that the smallest singular value is uniformly bounded away from zero for all geometries except the rectangle, whose aspect ratio deteriorates, and the slit problem. Note that, for the remaining geometries, it is not only the case that the condition number is *O*(1) as , but also its actual numerical value is less than 10.

## 8. Conclusions

In this paper, we established that difference equations defined on a finite lattice can be replaced with algebraic equations defined on the boundary only. The resulting BAEs are much smaller in size than the original problem and appear to give rise to well-conditioned linear algebraic problems. In particular, there are indications that the condition number for the BAE based on the discrete double-layer kernel is independent of the problem size. Furthermore, it should be possible to develop fast methods for computing matrix–vector products associated with BAEs and consequently fast iterative methods.

## Footnotes

- Received January 15, 2009.
- Accepted April 28, 2009.

- © 2009 The Royal Society