Numerical Solutions of the spectral problem for arbitrary self-adjoint extensions of the 1D Schrödinger equationThis work has been partially supported by the Spanish MICINN grant MTM2010-21186-C02-02, QUITEMAD programme P2009 ESP-1594. AI has been partially supported by Fundación Caja Madrid. JMP has been supported by the UCIIIM University through the PhD Program Grant M02-0910

Numerical Solutions of the spectral problem for arbitrary self-adjoint extensions of the 1D Schrödinger equation††thanks: This work has been partially supported by the Spanish MICINN grant MTM2010-21186-C02-02, QUITEMAD programme P2009 ESP-1594. AI has been partially supported by Fundación Caja Madrid. JMP has been supported by the UCIIIM University through the PhD Program Grant M02-0910

A. Ibort 222Department of Mathematics, Univ. of California at Berkeley, Berkeley CA 94720, USA. On leave of absence from Depto. de Matemáticas, Univ. Carlos III de Madrid, Avda. de la Universidad 30, 28911 Leganés, Madrid, Spain.    J. M. Pérez-Pardo 333Depto. de Matemáticas, Univ. Carlos III de Madrid, Avda. de la Universidad 30, 28911 Leganés, Madrid, Spain.
Abstract

A numerical algorithm to solve the spectral problem for arbitrary self-adjoint extensions of 1D regular Schrödinger operators is presented. It is shown that the set of all self-adjoint extensions of 1D regular Schrödinger operators is in one-to-one correspondence with the group of unitary operators on the finite dimensional Hilbert space of boundary data, and they are characterized by a generalized class of boundary conditions that include the well-known Dirichlet, Neumann, Robin, (quasi-)periodic boundary conditions, etc. The numerical algorithm is based on a nonlocal boundary extension of the finite element method and their convergence is proved. An appropriate basis of boundary functions must be introduced to deal with arbitrary boundary conditions and the conditioning of its computation is analyzed. Some significant numerical experiments are also discussed as well as the comparison with some standard algorithms. In particular it is shown that appropriate perturbations of standard boundary conditions for the free particle leads to the theoretically predicted result of very large absolute values of the negative groundlevels of the system as well as the localization of the corresponding eigenvectors at the boundary (edge states).

1 Introduction

The study of the self-adjointness of Schrödinger operators has been a fundamental mathematical problem since the beginning of Quantum Mechanics and there is a vast literature on the subject (see for instance [Re75], the review [Si00] and references therein). In spite of this, there is a continuous flow of new results and even surprises (see for instance the recent papers where some apparent paradoxical aspects of the spectrum of certain self-adjoint extensions of the Schrödinger operator in 2D are analyzed [Be08], [Ma09], [Be09]).

Consider the evolution of a quantum system on a -dimensional Riemannian manifold with boundary under the influence of a potential which is given by the Schrödinger equation with the Hamiltonian operator of the system given by

 H=−ℏ22mΔη+V(x)=−ℏ22m1√|η|∂∂xj√|η|ηjk∂∂xk+V(x). \hb@xt@.01(1.1)

Here stands for the Laplace-Beltrami operator on defined by the metric tensor given by , and . The second order differential operator is formally self-adjoint, however in order to define a unitary evolution of the quantum state , a self-adjoint extension of it must be specified. If denotes one of such self-adjoint extensions, because of Stone’s theorem, a one-parameter group of unitary operators exists such that and the quantum evolution of a given initial state will be uniquely determined as . Self-adjoint extensions of the operator are usually chosen by fixing the boundary values of the functions where the operator acts, typically Dirichlet or Neumann boundary conditions. However they are by no means the most general choice of boundary conditions determining self-adjoint extensions of the operator and a variety of other possibilities exist. The development of quantum information technologies makes relevant the discussion of more general classes of boundary conditions, i.e., of general self-adjoint extensions for the operator , as well as the numerical computation of their spectrum in order to integrate Schrödinger’s equation. For instance it has been recently proposed a physical implementation of a universal quantum computer by using scattering states on a quantum graph [Ci09], i.e., of a Schrödinger operator on a graph. It is also relevant to mention here the celebrated exponential algorithmic speedup by a quantum walk by Childs et al [Ci03]. The self-adjoint extensions defining a quantum system on a graph, otherwise called quantum Kirchoff’s rules [Ko00], [K03] are just one more instance of these possibilities. Another one is provided by the use of absorbing boundary conditions [Ba04] or by diverse optical cavities implementations [Kn03] of quantum walks.

In this paper we present an algorithm based on the finite element method to solve numerically the spectral problem for all  possible self-adjoint extensions of 1D regular Schrödinger operators. In this particular instance, the Schrödinger operator takes the form of a Sturm-Liouville operator. There exist many methods for approaching this problem (see [Ac09], [Ch99], [Ch09], [Fa57], [Pr73] just to mention a few of them), however they cannot be used to solve the problem of general self-adjoint extensions of it, as explained below. Most of the previous algorithms solve the wide class of self-adjoint extensions given by the equations of the form: , , where and are the endpoints of the interval. For instance, (quasi)-periodic boundary conditions can not be implemented with the above setting.

In this work we will take the approach developed in [As05] to describe the set of self-adjoint extensions of Laplace-Beltrami operators. We will show in section LABEL:Selfadjoint_extensions_of_Schrodinger_operators that in 1D the set of self-adjoint extensions of the Schrödinger operator is in one-to-one correspondence with the group of unitary operators on the space of boundary data for the problem and we will provide explicit constructions of such correspondence so that it can be actually used in the implementation of the corresponding numerical algorithms. Similar results are valid in higher dimensions, however there are additional technical difficulties related to the functional spaces where the domains of the operators are defined, that make the presentation and the discussion of the results more difficult. Because of this and in order to make the presentation of the results as simple and as short as possible we will consider in this paper just the 1D situation leaving the discussion in higher dimensions for later works. Thus the manifold will be assumed to be one-dimensional and compact and the potential function will be assumed to be regular and bounded below.

The numerical algorithm which is developed to solve the eigenvalue problem for all self-adjoint extensions of the Hamiltonian operator described in terms of boundary data is based on the finite element method, even though the treatment of the boundary basis functions is novel. The determination of the boundary basis function is analyzed as well as its conditioning. The stability and convergence of the method is proved by describing the corresponding self-adjoint extensions in weak form. A number of experiments are presented that are consistent with the obtained results and shows the efficiency of the corresponding computer algorithms. These results are compared with the standard methods in some significant examples showing again the good behavior of the proposed algorithm. In particular it is shown that for a family of self-adjoint extensions the groundlevel of the system must go to , a behaviour that is described nicely by the algorithm presented here as well as the presence of edge states, i.e., eigenfuctions which are localized at the boundary.

The paper is organized as follows. Section LABEL:Selfadjoint_extensions_of_Schrodinger_operators is devoted to prove the general theorems that will allow for a convenient description of the self-adjoint extensions of the Schrödinger operator in 1D in terms of boundary data. The main results of the finite element method for the eigenvalue problem for completely general self-adjoint extensions are discussed in section LABEL:FEM. An appropriate basis of nonlocal boundary functions is introduced there in order to implement the boundary conditions as described by the general theory. Moreover, in this section, results concerning the convergence and the stability of the numerical scheme are also proved. In section LABEL:Numerical_experiments_and_conclusions we present some numerical experiments displaying relevant features of the algorithm, as well as a comparison with other methods.

2 Self-adjoint extensions of Schrödinger operators

As it was stated in the introduction, we will restrict our attention to the case of Schrödinger operators on 1D compact manifolds and regular potentials bounded below. The Schrödinger operator for a particle moving on a smooth manifold with boundary and Riemannian metric is given by the Hamiltonian operator defined in Eq. (LABEL:hamiltonian). Notice first that a compact 1D manifold consists of a finite number of closed intervals , . Each interval will have the form and the boundary of the manifold is given by the family of points . Functions on are determined by vectors of complex valued functions . A Riemannian metric on is given by specifying a Riemannian metric on each interval , this is, by a positive smooth function on the interval , i.e., . Then the inner product on takes the form and the Hilbert space of square integrable functions on is given by .

On each subinterval the differential operator takes the form of a Sturm-Liouville operator

 Hα=−1Wαddxpα(x)ddx+Vα(x), \hb@xt@.01(2.1)

with smooth coefficients (in what follows we are taking the physical constants and such that ), , and is formally self-adjoint in the sense that

 ∫M¯¯¯¯ΨΔηΦ√ηdx=−∫Md¯¯¯¯ΨdxdΦdx√ηdx=∫M(¯¯¯¯¯¯¯¯¯¯¯ΔηΨ)Φ√ηdx, \hb@xt@.01(2.2)

for , complex-valued smooth functions with compact support in the interior of . In fact, the differential expression (LABEL:laplace) defines a symmetric operator on with dense domain . The operator is closable on and its minimal closed extension, that we will keep denoting by , has domain , the Sobolev space of order 2 on with functions and first derivatives vanishing at the boundary. The adjoint operator has dense domain on . In the particular instance of 1D the domain of the adjoint operator is exactly given by , the Sobolev space of order 2 on , i.e., the space of functions on possessing first and second weak derivatives that are square integrable (notice that in higher dimensions this is not true because the domain of the adjoint could contain functions which are not in , see for instance [Gr68]). Clearly and , thus the operator is symmetric but not self-adjoint on . Notice that the maximal extension of the operator , this is the operator defined in is not self-adjoint either because .

Hence we would like to determine whether there exist self-adjoint extensions of the symmetric operator , i.e., operators with domain such that , and . As it will become clear from the discussion below, there exist self-adjoint extensions of and all of them are determined by appropriate boundary conditions for the functions in .

Von Neumann’s theorem establishes (see for instance [We80], Thm. 8.12) that there is a one-to-one correspondence between self-adjoint extensions of the Laplace-Beltrami operator and unitary operators , where the deficiency spaces are defined as:

 N±={Ψ∈H2(M)∣Δ†ηΨ=±iΨ}. \hb@xt@.01(2.3)

Thus, given a unitary operator , the domain of the operator is given by

 D=D0⊕(I+K)N+, \hb@xt@.01(2.4)

and the extended operator takes the explicit form:

 ΔD(Ψ0⊕(I+K)ξ+)=ΔηΨ0⊕i(I−K)ξ+, \hb@xt@.01(2.5)

for all and .

Unfortunately von Neumann’s theorem is not always suitable for the numerical computation of general self-adjoint extensions of the Laplace-Beltrami operator because one needs to determine first the deficiency spaces and, moreover, the domain given by Eq. (LABEL:Neumanndomain) is algebraically not well behaved. We can take however a different route inspired in the classical treatment of formally self-adjoint differential operators. If we rewrite Eq. (LABEL:laplace) for functions , instead of , a simple computation shows:

 ∫M¯¯¯¯ΨΔηΦ√ηdx=∫M(¯¯¯¯¯¯¯¯¯¯¯ΔηΨ)Φ√ηdx+∫∂M(¯¯¯¯ψ˙φ−¯¯¯¯˙ψφ)vol∂η, \hb@xt@.01(2.6)

where stands for the induced Riemannian metric at the boundary, , , where is the exterior normal vector to , and respectively for and .

We thus obtain the Lagrange boundary form for the Laplace-Beltrami operator:

 Σ((ψ,˙ψ),(φ,˙φ))=∫∂M¯¯¯¯ψ˙φvol∂η−∫∂M¯¯¯¯˙ψφvol∂η=⟨ψ,˙φ⟩L2(∂M)−⟨˙ψ,φ⟩L2(∂M). \hb@xt@.01(2.7)

In what follows, if there is no risk of confusion, we will omit the subscript that denotes the inner product on the boundary manifold with respect to the measure defined by the volume form and we will simply write . The Lagrange boundary bilinear form defines a continuous bilinear form on the Hilbert space ,

 Σ((ψ1,ψ2),(φ1,φ2))=⟨ψ1,φ2⟩−⟨ψ2,φ1⟩,∀(ψ1,ψ2),(φ1,φ2)∈L2(∂M)⊕L2(∂M). \hb@xt@.01(2.8)

The trace map given by , is continuous and induces an isomorphism from onto (see for instance [Ad75], Thm. 7.20), where and are Sobolev spaces of fractional order. The previous observations provide a simple characterization of self-adjoint extensions of the operator in 1D.

Theorem 2.1

There is a one-to-one correspondence between self-adjoint extensions of the Laplace-Beltrami operator and non-trivial maximal closed isotropic subspaces of contained in . The correspondence being explicitly given by .

Proof. Let be the domain of a self-adjoint extension of the operator . Consider the subspace consisting on the set of pairs of functions that are respectively the restriction to of a function and its normal derivative. Notice that the subspace is closed in . Hence, because is an homeomorphism from to , is a closed subspace of . Because of Eq. (LABEL:boundary), it is clear that and finally the subspace is maximally isotropic in , because if this were not the case, there will be a closed isotropic subspace containing . Then defines a domain containing such that the operator would be symmetric on it, in contradiction with the self-adjointness assumption.

Let be a maximal closed -isotropic subspace in . Then consider the subspace of functions such that . It is clear that for any pair of functions , on , because is isotropic with respect to , Eq. (LABEL:boundary) gives and the operator is therefore symmetric in . Moreover, because of the maximality of in it is easy to see that , hence it is self-adjoint.

The key observation now is that, because , , where is the number of subintervals. Hence, the self-adjoint extensions of the Laplace-Beltrami operator are in one-to-one correspondence with the maximal isotropic subspaces of the Lagrange boundary form . It is a well known result (see [As05], [Br08] and references therein) that the maximal isotropic subspaces of the bilinear form (LABEL:lagrange2) are in one-to-one correspondence with the graphs of unitary operators through the relationship

 φ1−iφ2=U(φ1+iφ2).

This result is easily derived by realizing that the unitary transformation , defined by , transforms the Lagrange bilinear form into the bilinear form:

 ~Σ((φ+,φ−),(ψ+,ψ−))=−i[⟨φ+,ψ+⟩−⟨φ−,ψ−⟩], \hb@xt@.01(2.9)

where the notation has been used. Then it follows immediately that maximal isotropic closed subspaces of are mapped by into graphs of unitary operators of .

The last step in discussing all self-adjoint extensions of the Hamiltonian consists in realizing that because is compact the operator multiplication by a regular function is essentially self-adjoint and its unique self-adjoint extension has domain . Hence, the self-adjoint extensions of coincide with the self-adjoint extensions of .

We can summarize the preceding analysis by stating that, under the conditions above, the domain of a self-adjoint extension of the Schrödinger operator in 1D is defined as a closed subspace of functions on satisfying:

 ψ−i˙ψ=U(ψ+i˙ψ) \hb@xt@.01(2.10)

for a given unitary operator . Equation (LABEL:asorey) represent general quantum Kirchoff’s rules, as they define the most general self-adjoint extension of the quantum free motion on a family of intervals. Notice for instance that corresponds to Neumann’s boundary conditions and determines Dirichlet’s boundary conditions. Other unitary matrices will define different graphs, but many other will not. The formula above, Eq. (LABEL:asorey), provides a powerful and effective computational tool to deal with general self-adjoint extensions of Schrödinger operators, (even in dimensions higher than 1), and will be used extensively in the rest of the paper.

Once we have determined a self-adjoint extension of the Schrödinger operator , we can obtain its unitary evolution as indicated in the Introduction by computing the flow . It is well-known that the Dirichlet extension of the Laplace-Beltrami operator has a pure discrete spectrum because of the compactness of the manifold and the ellipticity of the operator, hence all self-adjoint extensions have a pure discrete spectrum (see [We80], Thm. 8.18). Then the spectral theorem for the self-adjoint operator states:

 HU=∞∑k=1λkPk,

where is the orthogonal projector onto the finite-dimensional eigenvector space corresponding to the eigenvalue . The unitary flow is given by:

 Ut=∞∑k=1e−itλk/ℏPk.

Hence all that remains to be done is to solve the eigenvalue problem:

 HUΨ=λΨ, \hb@xt@.01(2.11)

for the Schrödinger operator . We devote the rest of this paper to provide an efficient numerical algorithm to solve Eq. (LABEL:eigenvalue_problem).

3 Finite element method for the eigenvalue problem in the interval

Our aim now is to find a numerical algorithm to solve the spectral problem for the self-adjoint extension of the 1D Schrödinger operator determined by the unitary matrix in (LABEL:eigenvalue_problem) and given explicitly by the differential system:

 ⎧⎪ ⎪⎨⎪ ⎪⎩−ΔηΦ+VΦ=λΦ,Φ∈H2(M),φ−i˙φ=U(φ+i˙φ),φ=Φ∣∣∂M,˙φ=dΦdν∣∣∂M,U∈U(2n).. \hb@xt@.01(3.1)

From now on we will keep the dot notation for the outward normal derivative at the boundary and use primes to denote the standard derivative . Because of the particularly simple form that the Lagrange boundary form takes in one dimension, it will be convenient to introduce the following notation (check with the right-hand side of (LABEL:lagrange)):

 [¯¯¯¯ΨΦ′]∣∣∂M:=⟨ψ,˙φ⟩L2(∂M)=n∑α=1¯¯¯¯Ψ(bα)Φ′(bα)−¯¯¯¯Ψ(aα)Φ′(aα).

It is convenient to look for weak solutions of (LABEL:schrodinger). Taking the inner product of (LABEL:schrodinger) with a vector on the dense domain of smooth functions on and integrating by parts we obtain easily the corresponding weak problem:

 ⎧⎪ ⎪⎨⎪ ⎪⎩⟨Ψ′,Φ′⟩−[¯¯¯¯ΨΦ′]∂M+⟨Ψ,VΦ⟩=λ⟨Ψ,Φ⟩∀Ψ∈C∞(M),φ−i˙φ=U(φ+i˙φ),φ=Φ∣∣∂M,˙φ=dΦdν∣∣∂M,U∈U(2n). \hb@xt@.01(3.2)

It is easy to show that is a solution of (LABEL:schrodinger) iff it solves (LABEL:weak). In the next subsections we develop an algorithm that approximates such solutions and, moreover, we prove its stability and convergence.

3.1 Finite elements for general self-adjoint boundary conditions

To solve (LABEL:weak) we construct a family of finite-dimensional subspaces of functions of satisfying the boundary conditions (LABEL:asorey). Such finite-dimensional subspaces will be constructed using finite elements. The finite element model that we use is given by the unit interval, the space of linear polynomials on , and the vertex set .

The domain of our problem is the manifold which consists of the disjoint union of the intervals , . For each we will construct a nondegenerate subdivision as follows. Let be the integer defined as , where denotes the integer part of , , and . We will assume that each , and . Let us denote by the multi index . Then satisfies

 N≤|r|≤N+n.

Now we will subdivide each interval into subintervals of length:

 hα=Lαrα+1.

The nondegeneracy condition supposes that it exists such that for all and for all

 diamBIα,a≥ρdiamIα,a,

where is the largest ball contained in such that is star shaped with respect to . In our particular case this is satisfied trivially since

 diamBIα,a=diamIα,a.

(It would be possible to use a set of independent steps , one for each interval; however, this could create some technical difficulties later on that we prevent in this way.) Each subinterval contains nodes that will be denoted as

 x(α)k=aα+hαk,k=0,…,rα+1.

3.1.1 Bulk functions

Consider now the family of piecewise linear functions , , that are zero at all nodes except at the th node of the interval , where it has value , or more explicitly,

 f(α)k(x)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩s,x=x(α)k−1+shα,0≤s≤1,1−s,x=x(α)k+shα,0≤s≤1,0,otherwise.

Notice that these functions are differentiable on each subinterval. All these functions satisfy trivially the boundary conditions (LABEL:asorey) because they and their normal derivatives vanish at the endpoints of each interval. They are localized around the inner nodes of the intervals. We will call these functions bulk functions.

3.1.2 Boundary functions

We will add to the set of bulk functions a family of functions that implement nontrivially the boundary conditions determining the self-adjoint extension. These functions will be called boundary functions and the collection of all of them will be denoted by . Contrary to bulk functions, boundary functions need to be “delocalized” so that they can fulfill any possible self-adjoint extension’s boundary condition.

Because the endpoints , of the intervals and the adjacent nodes, and , are going to play a prominent role in what follows, we introduce some notation that takes care of them. We will consider an index that labels the endpoints of the intervals. Now for each vector consider the following functions (see Figure LABEL:function_beta):

 β(w)(x)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩v2α−1+s(w2α−1−v2α−1),x=x(α)0+shα,w2α−1(1−s),x=x(α)1+shα,sw2α,x=x(α)rα−1+shα,w2α+s(v2α−w2α),x=x(α)rα+shα,,0,x(α)2≤x≤x(α)rα−1,0≤s≤1,1≤α≤n.

Each function of the previous family is determined (apart from the vector ) by the vector that collects the values of at the endpoints of the subintervals. If we denote by the vectors such that , , the vectors are just the standard basis for . The corresponding functions will now be denoted simply by . Notice that each boundary function is completely characterized by the unique nonextremal node where it does not vanish111If , it is the node , and if , it corresponds to the node . and the values at the endpoints. We denote by , , the boundary values of the functions above.

3.1.3 The boundary matrix

The extremal values of the boundary functions are undefined, but we are going to show that the conditions (LABEL:asorey) imposed on the boundary functions constitute a determinate system of linear equations for them. Because the boundary functions are constructed to be piecewise linear, the normal derivatives of these functions at the boundary can be obtained easily. For the left boundaries of the intervals, i.e., at the points , we have

 dβ(k)dν∣∣∣x=aα=−1hα(w(k)2α−1−v(k)2α−1), \hb@xt@.01(3.3)

and respectively, for the right boundaries,

 dβ(k)dν∣∣∣x=bα=1hα(v(k)2α−w(k)2α)=−1hα(w(k)2α−v(k)2α). \hb@xt@.01(3.4)

Thus the vector containing the normal derivatives of the function , consistently denoted by , is given by

 ˙β(k)l=−1hl(w(k)l−v(k)l)=−1hl(δlk−v(k)l), \hb@xt@.01(3.5)

where we use again the consistent notation , if , or . For each boundary function , the boundary conditions (LABEL:asorey) read simply as the system of equations on the components of the vector ,

 v(k)−i˙β(k)=U(v(k)+i˙β(k)),

or, componentwise,

 v(k)l(1−ihl)+ihlw(k)l=2n∑j=1Ulj[v(k)j(1+ihj)−ihjw(k)j],l=1,…,2n.

Collecting coefficients and substituting the expressions we get

 2n∑j=1[(1−ihj)δlj−Ulj(1+ihj)]v(k)j=[−ihkδlk−Ulkihk]. \hb@xt@.01(3.6)

This last equation can be written as the matrix linear system:

 FV=C \hb@xt@.01(3.7)

with a matrix whose entries are given by , . The th column of contains the boundary values of the boundary function . The matrix with entries

 Flj=(1−ihj)δlj−Ulj(1+ihj),

will be called the boundary matrix of the subdivision of the domain determined by the integer , and

 Clk=−ihk(δlk+Ulk)

defines the inhomogeneous term of the linear system (LABEL:boundary_equation). Using a compact notation we get

 F=diag(1−i/h)−Udiag(1+i/h),C=−i(I+U)diag(1/h)

where denotes the vector whose components are . Notice that depends just on and the integer defining the size of the discretization leading to the approximate weak spectral problem defined below.

3.2 Conditioning of the boundary matrix

Before addressing the construction of the approximate spectral problem we will study the behavior of system (LABEL:boundary_equation) under perturbations; in other words, we will compute the condition number of the boundary matrix and show that it is small enough to ensure the accuracy of the numerical determination of our family of boundary functions . The relative condition number we want to compute is

 κ(F)=∥F∥∥F−1∥.

In our case, the boundary matrix can be expressed as

 F=¯D−UD=(I−UD¯D−1)¯D

with . Notice that the product is a unitary matrix which we will denote as or simply if we do not want to emphasize the dependence of . Thus, and

 ∥F∥=∥(I−U0)¯D∥≤∥I−U0∥∥¯D∥≤2∥D∥.

On the other hand,

 ∥F−1∥=∥¯D−1(I−U0)−1∥≤∥¯D−1∥∥(I−U0)−1∥=∥D−1∥minλ∈spec(U0){|1−λ|}

and thus we obtain

 κ(F)≤κ(D)2minλ∈spec(U0){|1−λ|}. \hb@xt@.01(3.8)

As is a diagonal matrix its condition number is given by

 κ(D)=√1h2min+1√1h2max+1≤hmaxhmin

with () the biggest (smallest) step of the discretization determined by . We get finally,

 κ(F)≤hmaxhmin2|1−λ| \hb@xt@.01(3.9)

with the closest element of the spectrum of to . Of course, because is unitary, it may happen that is in its spectrum, so that the condition number is not bounded. Because the matrix depends on , its eigenvalues will depend on too. We want to study the dependence of the closest eigenvalue to 1, or 1 for that matter, with respect to perturbations of the vector .

Lemma 3.1

Suppose that is an eigenvector with eigenvalue of and that the perturbed matrix , for small enough, is such that . Then to first order in .

Proof. Clearly, if and is small enough, there exist a vector , with , such that . Then we have

 U0δX+δUX0+δUδX=δX.

Because and is unitary, , and then multiplying on the left by and keeping only first order terms, we get the desired condition:

Because of the previous lemma, if is a unitary perturbation of such that , for any eigenvector with eigenvalue 1, then . Now if we consider a unitary perturbation of such that we want to estimate how far 1 is from the spectrum of . Consider the eigenvalue equation for the perturbed matrix. The perturbed eigenvalue will satisfy

 (U0+δU)(X0+δX)=(1+δλ)(X0+δX). \hb@xt@.01(3.10)

Multiplying on the left by and solving for it follows that

 |δλ|≥|XH0δUX0+XH0δUδX||1+XH0δX|≥|XH0δUX0|−|XH0δUδX|1+∥δX∥ \hb@xt@.01(3.11)

for small enough. Taking into account the particular form of the matrix we have that and therefore . Moreover, because is a diagonal matrix, and its singular values are the modulus of its diagonal entries. Hence we have the following proposition.

Proposition 3.2

Given the matrix with eigenvalue 1, then 1 is not an eigenvalue of any unitarily perturbed matrix with small enough and . Moreover there exists a constant such that the perturbation of such eigenvalue satisfies the lower bound,

 |δλ|≥σmin(δ(D¯D−1))−Cσ2max(δ(D¯D−1))1+Cσmax(δ(D¯D−1))>0.

Proof. Because of Lemma LABEL:perturbation_1 it is sufficient to show that . But this is an easy consequence of the fact that . Furthermore there exists a constant such that ; hence taking small enough we get and the bound follows from (LABEL:bounddelta).

Now we can apply Proposition LABEL:lower_bound to (LABEL:KF) and if we neglect terms and we finally get the desired bound for the condition number

 κ(F)≤hmaxhmin1δh. \hb@xt@.01(3.12)

Then, if for a given we obtain a boundary matrix which is bad conditioned, it suffices to change the size of the discretization, i.e., to increase , to improve the condition number. Of course, if is already quite large, then the bound (LABEL:kappa_h) could be useless. For typical values , it can be taken as to provide condition numbers .

3.3 The spectral pencil

For any we define the finite-dimensional approximation space as the linear span of the bulk and boundary functions, i.e., . All functions and are linearly independent; thus the dimension of will be . It is convenient to rearrange the elements of the basis above as follows:

 β(1),f(1)1,…,f(1)r1,β(2),β(3),f(2)1,…,f(2)r2,β(4),…,β(2n−1),f(n