text:
# SPECTRAL POLLUTION

## 1 Introduction

It is well known that computing the eigenvalues of a self-adjoint bounded or differential operator $`H`$ acting in an infinite dimensional space $``$ is not straightforward. The standard truncation (or Ritz) method is to take a large finite-dimensional test function space $``$ with an orthonormal basis $`\{\varphi _r\}_{r=1}^n`$ and to compute the eigenvalues of the $`n\times n`$ matrix $`A`$ given by

$$A_{r,s}=H\varphi _r,\varphi _s.$$

Unfortunately one often encounters the problem of spectral pollution: the matrix may possess eigenvalues which are unrelated to any spectral properties of the original operator. This is an important problem in several areas of applied mathematics, as is illustrated by recent papers, . There are methods of detecting and resolving such problems, but they are not easy to understand and for that reason are not well appreciated. Our main result in this paper is the proof that, at least in some situations, the geometrically motivated method of Davies,  is mathematically equivalent to the method of Zimmermann and Mertins,  which is based on Lehmann’s method  and its extension by Goerisch . We also provide various other ideas and examples to motivate the method and show how well it works.

## 2 Formulation of the Problem

We describe the context of the paper at a rather abstract level. The assumptions are designed to apply to a range of problems, including those involving second order elliptic operators, as will become clear in our applications.

Let $`A`$ be a bounded, self-adjoint, linear operator acting on a Hilbert space $``$, and let $`P`$ be the orthogonal projection onto a closed linear subspace $``$ of $``$. The truncation of $`A`$ to $``$, denoted $`A_{}`$, is defined to be the operator $`PAP`$ restricted to $``$. (Similar definitions may be made for unbounded $`A`$ and for quadratic forms.) If $``$ is sufficiently large, one might hope that the spectrum of $`A_{}`$ is close to that of $`A`$. The Rayleigh-Ritz, or variational, theorems establish precisely this for the part of the spectrum of $`A`$ outside $`[\alpha ,\beta ]`$, where $`\alpha `$ (resp. $`\beta `$) is the smallest (resp. biggest) point in the essential spectrum of $`A`$.

Within $`[\alpha ,\beta ]`$ the situation is far less pleasant. One says that $`\lambda `$ is a point of spectral pollution of $`A`$ for the increasing sequence of subspaces $`_n`$ with union dense in $``$ if there exist $`\lambda _n\mathrm{Spec}(A__n)`$ such that $`\lambda _n\lambda `$ as $`n\mathrm{}`$ but $`\lambda \mathrm{Spec}(A)`$. We will see that spectral pollution is commonplace, but the following two elementary lemmas sometimes enable one to determine $`\mathrm{Spec}(A)`$ from its truncations.

###### Lemma 1

If $`(IP)AP<\epsilon `$ and $`\lambda \mathrm{Spec}(A_{})`$ then

$$(\lambda \epsilon ,\lambda +\epsilon )\mathrm{Spec}(A)\mathrm{}.$$

Proof One considers $`A`$ as a perturbation of $`PAP+(IP)A(IP)`$.

###### Lemma 2

If $`f=1`$ and $`Af\lambda f<\epsilon `$ then

$$(\lambda \epsilon ,\lambda +\epsilon )\mathrm{Spec}(A)\mathrm{}.$$

Proof This uses the resolvent norm estimate

$$(A\lambda I)^1=\mathrm{dist}\{\lambda ,\mathrm{Spec}(A)\}^1.$$

Lemma 2 is originally due to D. Weinstein. Methods of estimating $`\mathrm{Spec}(A)`$ starting from a subspace $``$ may be divided into two general categories. Prior methods involve proving that if $`_n`$ is a particular increasing sequence of subspaces then the spectrum of $`A__n`$ converges to that of $`A`$. There should also be an explicit error bound, or at least some information about the rate of convergence. When applied to differential operators such methods often involve finite element subspaces and Sobolev space embedding theorems. They have two disadvantages. The first is that the convergence proofs are often very difficult and assume conditions which are specific to the particular applications. The other is that the error bounds are often far larger than the true errors. Nevertheless valuable and deep theoretical results can sometimes be obtained using such methods, .

Posterior methods, on the other hand, involve choosing a subspace $``$ and carrying out computations which provide rigorous bounds on the location of $`\mathrm{Spec}(A)`$. If these are good enough one stops, but if they are not one chooses a different or larger subspace $``$. The problem here is that there is no prior guarantee that useful bounds will ever be obtained. Nevertheless in practice the method is often a very good one. Besides the result by D. Weinstein (Lemma 2), the more refined spectral bounds by T. Kato , Aronszajn-A.Weinstein , Lehmann  and Goerisch  fall into this category. Two particular posterior methods have recently been the focus of some attention. One uses the second order spectrum, introduced in , and is the subject of the companion paper by Levitin and Shargorodsky, . This paper studies the other, which depends upon a systematic study of how certain residuals depend upon a real parameter, . We describe a close relationship between this method and that of Zimmermann and Mertins, .

We next set the technical context of the paper. Let $`H`$ be a non-negative unbounded self-adjoint operator acting on $``$, and suppose that $`0\mathrm{Spec}(H)`$. Now suppose that $`K`$ is a second non-negative self-adjoint operator, possibly with a different operator domain, but satisfying a quadratic form relative bound. In other words if $`Q_H`$ and $`Q_K`$ are their two quadratic forms then $`\mathrm{Dom}(Q_H)\mathrm{Dom}(Q_K)`$ and

$$0Q_K(f)aQ_H(f)$$

for all $`f\mathrm{Dom}(Q_H)=\mathrm{Dom}(H^{1/2})`$. (In fact $`Q_K`$ need not be closed, and hence need not be determined by an operator $`K`$, in the analysis below.)

The eigenvalue problem $`\lambda Hf=Kf`$ may be recast in the weak, or variational, form. In other words one seeks $`f`$, $`\lambda `$ such that

$`\lambda Q_H(f,g)=Q_K(f,g)`$ (1)
for all $`g\mathrm{Dom}(Q_H)`$.

In many cases the operator $`K`$ is of the following form. There exists a closed operator $`D:`$ with operator domain $`\mathrm{Dom}(Q_H)`$ such that $`H=D^{}D`$ and $`K=D^{}MD`$, where $`M`$ is a non-negative bounded operator acting on $``$. Under our hypotheses $`\mathrm{Ker}(D)=0`$ and $`\mathrm{Ran}(D)`$ is closed. There is no loss of generality in assuming that $`\mathrm{Ran}(D)=`$, by simply reducing the size of $``$ and restricting $`M`$. In this case the issue is to find $`f`$ such that

$$\lambda Df,Dg=MDf,Dg$$

for all $`g`$, or equivalently to find the eigenvalues of the bounded operator self-adjoint $`M`$.

Many non-negative differential operators $`H`$ of order $`2n`$ can be written in the form $`H=D^{}D`$, where $`D`$ is a differential operator of order $`n`$. Indeed, in a formal sense the assumption is always satisfied with $`D=H^{1/2}`$ and $`M=H^{1/2}KH^{1/2}`$. This observation is rarely useful, because $`H^{\pm 1/2}`$ are usually not easy to describe in explicit terms. If $`H`$ and $`K`$ are both differential operators, and $`K`$ is of lower order than $`H`$, then one often finds that $`M`$ is compact. Its essential spectrum is then simply $`\{0\}`$, and the pollution problem does not arise.

Before continuing, we mention that for Sturm-Liouville equations one can often bypass the problems in this paper by solving the initial value problem and computing the Prufer angles of the solution. This can enable one to locate eigenvalues in a gap of the essential spectrum. The method is not available in higher dimensions.

## 3 Reality of the Phenomenon

Experience shows that spectral pollution is not just a theoretical problem. It occurs on a regular basis for operators of physical importance, such as truncations of periodic Schrödinger operators to finite regions; see for example . In this section we demonstrate how bad the phenomenon can be. Later in the paper we produce more realistic examples.

The following well-known argument (see e.g. ) shows that spectral pollution can even occur in an extreme manner for basis functions which are finite linear combinations of exact eigenfunctions, and which are therefore in the operator domain.

###### Lemma 3

Let $`(\lambda _n^{(0)})_{n𝐍}`$ and $`(\lambda _n^{(1)})_{n𝐍}`$ be two (possibly constant) sequences of eigenvalues of problem (1), converging to limits $`\lambda ^{(0)}`$ and $`\lambda ^{(1)}`$, respectively, satisfying $`\lambda ^{(0)}<\lambda ^{(1)}`$. Moreover, let $`\mu `$ satisfy $`\lambda ^{(0)}<\mu <\lambda ^{(1)}`$. Then, there exists a sequence $`(\psi _n)`$ of trial functions which all are finite linear combinations of exact eigenfunctions of (1), such that

$`Q_H(\psi _n,\psi _m)=\delta _{nm},Q_K(\psi _n,\psi _m)=\mu \delta _{nm},`$ (2)
i.e., $`\mu `$ occurs as the only eigenvalue of the Ritz method, for every dimension.

Proof Without loss of generality, let $`\lambda _n^{(0)}<\mu <\lambda _n^{(1)}`$ for all $`n𝐍`$, and let $`(f_n^{(0)})`$ and $`(f_n^{(1)})`$ denote sequences of eigenfunctions of (1) corresponding to $`(\lambda _n^{(0)})`$ and $`(\lambda _n^{(1)})`$, respectively, such that

$`Q_H(f_n^{(i)},f_m^{(i)})=\delta _{nm},Q_K(f_n^{(i)},f_m^{(i)})=\lambda _n^{(i)}\delta _{nm}.`$
Finally, choose $`\alpha _n[0,2\pi )`$ such that

$`(\mathrm{cos}^2\alpha _n)\lambda _n^{(0)}+(\mathrm{sin}^2\alpha _n)\lambda _n^{(1)}=\mu ,`$
and

$`\psi _n:=(\mathrm{cos}\alpha _n)f_n^{(0)}+(\mathrm{sin}\alpha _n)f_n^{(1)}.`$
(2) is then quite obvious.

In the same way, also arbitrary sequences $`(\mu _n)`$ between $`(\lambda _n^{(0)})`$ and $`(\lambda _n^{(1)})`$ can be shown to occur as spurious eigenvalues.

If we make the weaker assumption that $`\lambda ^{(0)}`$ and $`\lambda ^{(1)}`$ are in the essential spectrum of $`M`$, then a similar argument can be followed. One can ensure that $`D\psi _n`$ lie in the domain of all powers of the operator $`M`$.

One may criticize the above example on the grounds that the sequence $`(\psi _n)`$ of basis functions constructed above may be far from complete. This gap can be closed as follows, at least if problem (1) has a $`Q_H`$-orthonormal and complete system of eigenfunctions:

In this case, let $`(f_n^{(2)})_{n𝐍}`$ denote a third sequence of eigenfunctions such that $`(f_n^{(0)}),(f_n^{(1)})`$, and $`(f_n^{(2)})`$ together are $`Q_H`$-orthonormal and complete. Then, with $`(\alpha _n)`$ as above, define

$`\phi _1:=(\mathrm{cos}\alpha _1)f_1^{(0)}+(\mathrm{sin}\alpha _1)f_1^{(1)},`$
and for $`n1`$:

$`\phi _{4n2}`$ $`:=`$ $`(\mathrm{sin}\alpha _n)f_n^{(0)}+(\mathrm{cos}\alpha _n)f_n^{(1)},`$
$`\phi _{4n1}`$ $`:=`$ $`(\mathrm{cos}\alpha _{2n})f_{2n}^{(0)}+(\mathrm{sin}\alpha _{2n})f_{2n}^{(1)},`$
$`\phi _{4n}`$ $`:=`$ $`(\mathrm{cos}\alpha _{2n+1})f_{2n+1}^{(0)}+(\mathrm{sin}\alpha _{2n+1})f_{2n+1}^{(1)},`$
$`\phi _{4n+1}`$ $`:=`$ $`f_n^{(2)}.`$
Then, $`(\phi _n)_{n𝐍}`$ is $`Q_H`$-orthonormal and complete, and the Ritz method (with basis functions $`\phi _n`$) of dimension $`N`$ provides $`\mu `$ as a $`\frac{N}{4}`$-fold eigenvalue.

## 4 Methods of Avoiding Pollution

There are various ways of identifying and avoiding spurious eigenvalues. Suppose that $`(\mathrm{Dom}(Q_H),Q_H)`$ is a Hilbert space, with continuous embedding into $``$. Let $`(_n)_{n𝐍}`$ denote an increasing sequence of subspaces, with $`_{n𝐍}_n`$ dense in $`\mathrm{Dom}(Q_H)`$. Let $`(\mu _n,\chi _n)_{n𝐍}`$ be a sequence of approximate eigenpairs to the problem (1) obtained by the Ritz method:

$`\chi _n_n,\mu _nQ_H(\chi _n,g)=Q_K(\chi _n,g)\text{ for all }g_n.`$ (3)
###### Lemma 4

Suppose that $`\mu _n\mu 𝐑`$ (or that a subsequence converges to $`\mu `$), with $`\mu `$ not an eigenvalue of (1); $`\mu `$ is therefore either in the essential spectrum or a point of spectral pollution.

Then, if we normalize $`\chi _n`$ by

$`\text{ }\chi _n\text{ }_H:=\sqrt{Q_H(\chi _n,\chi _n)}=1,`$ (4)
we have

$`\chi _n0\text{ weakly in }(\mathrm{Dom}(Q_H),Q_H).`$ (5)
Note If in particular the embedding $`\mathrm{Dom}(Q_H)`$ is compact, we therefore obtain $`\chi _n0`$ in $`(,.)`$, i.e.,

$`{\displaystyle \frac{\chi _n}{\text{ }\chi _n\text{ }_H}}0.`$ (6)
Since in practice one often performs the Ritz method for various different dimensions $`n`$, one has the chance of getting evidence for (6) by evaluating the quotient for various $`n`$.

Proof We show that each subsequence of $`(\chi _n)`$ has a subsequence converging weakly to zero. Indeed, (4) implies that a subsequence $`(\chi _{n_k})`$ of a given subsequence converges weakly in $`\mathrm{Dom}(Q_H)`$ to some $`\chi \mathrm{Dom}(Q_H)`$. Since $`Q_K0`$ is bounded in $`(\mathrm{Dom}(Q_H),Q_H)`$, this implies $`Q_K(\chi _{n_k},g)Q_K(\chi ,g)`$ for each $`g\mathrm{Dom}(Q_H)`$. Consequently,

$`\mu _{n_k}Q_H(\chi _{n_k},g)Q_K(\chi _{n_k},g)\mu Q_H(\chi ,g)Q_K(\chi ,g)`$ (7)
for each $`g\mathrm{Dom}(Q_H)`$. By (3), the left-hand side of (7) is zero for $`g_n`$ and $`n_kn`$. Therefore, (7) implies

$`\mu Q_H(\chi ,g)=Q_K(\chi ,g)\text{ for all }g{\displaystyle \underset{n𝐍}{}}_n`$
and thus, by density, for all $`g\mathrm{Dom}(Q_H)`$. Since $`\mu `$ is no eigenvalue of (1), the assertion $`\chi =0`$ follows.

One may find the eigenvalues of a bounded operator $`M`$ in a gap $`(\alpha ,\beta )`$ of its essential spectrum by considering instead the operator $`N=(M\rho I)^1`$, for some $`\rho (\alpha ,\beta )`$. The eigenvalues of $`M`$ in $`(\alpha ,\beta )`$ are associated to eigenvalues of $`N`$ which lie outside its essential spectrum; the latter may therefore be found by using the Rayleigh-Ritz procedure. An apparent disadvantage of this approach is the need to find an expression for the operator $`N`$, but this may be avoided, .

In this paper we explore a geometrically motivated method of finding eigenvalues, which traces its origins back to Kato, , and Weidmann, \[6, Cor. 6.20\]. One of our main goals is to describe its close relationship to ideas of Goerisch and Lehmann, and to developments by various other authors, . Suppose that the self-adjoint operator $`M`$ acting on the Hilbert space $``$ has a single eigenvalue $`\mu `$ which lies within a gap $`(\alpha ,\beta )`$ in its spectrum. If

$$|\lambda \mu |<\mathrm{min}\{|\lambda \alpha |,|\lambda \beta |\}$$

then one may apply the Rayleigh-Ritz method to the operator $`B=(M\lambda I)^2`$ for any $`\lambda `$ close enough to $`\mu `$. The smallest eigenvalue of $`B`$ will then be $`(\lambda \mu )^2`$, from which $`\mu `$ can be calculated (up to a square root ambiguity). The idea is to choose the value of $`\lambda `$ to maximize the numerical efficiency of the method. One’s intuition is to take $`\lambda `$ as close as possible to the eigenvalue $`\mu `$, but we shall see that this is not the correct strategy.

As in , we investigate the $`\lambda `$-dependence of the procedure above by introducing the function

$$F(\lambda )=\mathrm{min}\{Mf\lambda f:f\text{ and }f=1\}.$$

If one could evaluate $`F`$ numerically then the formula

$$F(\lambda )=\mathrm{dist}\{\lambda ,\mathrm{Spec}(M)\}$$
(8)
would enable one to locate the spectrum of $`M`$. To find the eigenvalues of $`M`$ one needs to combine this with some theoretical input which distinguishes between the point spectrum and essential spectrum of $`M`$.

Following , we introduce the approximate functions

$$F_n(\lambda )=\mathrm{min}\{Mf\lambda f:f_n\text{ and }f=1\}.$$

These may be used to obtain useful spectral information, as follows.

###### Lemma 5

If $`M`$ is bounded then for any $`\epsilon >0`$ there exists $`N_\epsilon `$ such that $`nN_\epsilon `$ implies

$$F(s)F_n(s)F(s)+\epsilon $$
(9)
for all $`s𝐑`$.

Proof It follows from their definition that $`F_n(s)`$ decrease monotonically to $`F(s)`$ as $`n\mathrm{}`$ for each $`s𝐑`$. It also follows from its definition that

$$|F_n(\lambda )F_n(\mu )||\lambda \mu |$$

for all $`\lambda ,\mu 𝐑`$. Equivalently $`|F^{}(\lambda )|1`$ for almost all $`\lambda 𝐑`$ (in the sense of Lebesgue). Since the functions $`F_n`$ are equicontinuous they must converge uniformly on compact subsets of $`𝐑`$. Suppose that $`nN_\epsilon `$ implies (9) for $`MsM`$. If $`s>M`$ then

$`F(s)`$ $``$ $`F_n(s)`$
$``$ $`F_n(M)+sM`$
$``$ $`F(M)+\epsilon +sM`$
$`=`$ $`F(s)+\epsilon .`$
A similar argument applies if $`s<M`$.

###### Lemma 6

Given (9), let $`\lambda `$ be an eigenvalue of $`M`$ and let there be no other point of $`\mathrm{Spec}(M)`$ in $`[\lambda 4\epsilon ,\lambda +4\epsilon ]`$. Then $`F_n`$ has a local minimum at some point $`\sigma [\lambda \epsilon ,\lambda +\epsilon ]`$. If $`\mu =\sigma F_n(\sigma )`$ and $`\nu =\sigma +F_n(\sigma )`$ then

$$\lambda 2\epsilon \mu \lambda \nu \lambda +2\epsilon .$$

Proof Since $`F(\lambda )=0`$ we see that $`0F_n(\lambda )\epsilon `$. However if $`\epsilon <|s\lambda |2\epsilon `$ we have

$$\epsilon <|s\lambda |=F(s)F_n(s).$$

Hence the minimum of $`F_n`$ within the interval $`[\lambda 2\epsilon ,\lambda +2\epsilon ]`$ must occur for $`\sigma `$ within $`[\lambda \epsilon ,\lambda +\epsilon ]`$. Moreover $`0F_n(\sigma )\epsilon `$. The second statement of the lemma follows immediately.

The functions $`F_n()`$ may be computed as follows. One introduces

$$F_{m,n}(\lambda )=\mathrm{min}\{P_m(Mf\lambda f):f_n\text{ and }f=1\}$$

where $`mn`$. If $`m=n`$ then this achieves no more than finding the eigenvalues of the truncation of $`M`$ to $`_n`$ (i.e. applying the Ritz method), but if $`mn`$ it may be far superior. One may also compute $`F_n(\lambda )`$ exactly by replacing $`_m`$ above by the finite-dimensional subspace

$$\stackrel{~}{}_n=M(_n)+_n.$$

The projection $`P_m`$ in the argument below is then replaced by the projection $`\stackrel{~}{P}_n`$ onto $`\stackrel{~}{}_n`$. Whether or not this Lanczos process is easily implemented depends upon the example.

The evaluation of $`F_{m,n}(\lambda )`$ is straightforward. One has the formula

$$F_{m,n}(\lambda )=G(m,n,\lambda )^{1/2}$$

where $`G(m,n,\lambda )`$ is the smallest eigenvalue of the $`n\times n`$ matrix

$$B(m,n,\lambda )=P_nMP_mMP_n2\lambda P_nMP_n+\lambda ^2P_n.$$

In practice one finds potential eigenvalues $`\lambda `$ of $`M`$ by using the Ritz method. One then examines $`F_{m,n}(\lambda )`$ for those values. If it is sufficiently small then one has good evidence that there is indeed an eigenvalue near the candidate value, while if it is not small then one has no such evidence. Using only a finite number of matrix entries of the operator $`M`$ no more could ever be established, but this method is far better than simply using the Ritz method without any checks.

The above provides evidence for the existence of an eigenvalue near a specified number. One can also obtain evidence for its multiplicity, or more precisely for the number of eigenvalues in a very short interval, by an elaboration of the above ideas which is described in .

## 5 Bounding the eigenvalues

In  we discussed how to obtain rigorous bounds on the spectrum of $`M`$ (or of problem (1), respectively) from the functions $`F_n()`$, thus avoiding the problem of spectral pollution. In this paper we develop the ideas further. If the subspace $`_n`$ is large enough then $`F_n(\lambda )`$ has local minima near the eigenvalues of $`H`$, and its derivatives are close to $`\pm 1`$ on either side of the eigenvalues. The following exact formula for the derivative may be useful, for example if applying Newton’s method. We have

$$F_n^{}(\lambda )=\frac{G_n^{}(\lambda )}{2\sqrt{G_n(\lambda )}}$$

where $`G_n(\lambda )`$ is the smallest eigenvalue of the $`n\times n`$ matrix

$`B_n(\lambda )`$ $`=`$ $`P_nM^{}MP_n2\lambda P_nMP_n+\lambda ^2P_n`$ (10)
$`=`$ $`D_n2\lambda M_n+\lambda ^2I`$ (11)
where the second line is understood as acting in $`_n`$ and $`M_n=M__n`$. If $`f_n(\lambda )`$ is the normalized eigenvector corresponding to the eigenvalue $`G_n(\lambda )`$ then it is well-known that

$`G_n^{}(\lambda )`$ $`=`$ $`B_n^{}(\lambda )f_n(\lambda ),f_n(\lambda )`$ (12)
$`=`$ $`2\lambda 2M_nf_n(\lambda ),f_n(\lambda ).`$ (13)
A computation of the graph of the function $`F_n()`$ is bound to be slow if $`n`$ is large. It is best therefore to attempt this only for reasonably small $`n`$, in order to understand approximately where the spectrum is located. One then passes to larger values of $`n`$ for which $`F_n(\lambda )`$ is only computed for carefully selected values of $`\lambda `$.

We next discuss the computation of the graph of $`F_n()`$ for reasonably small $`n`$. The procedure below is guided by the fact that it is optimal if $`D_n=M_n^2`$, or equivalently if $`(IP_n)MP_n=0`$. It should therefore be close to optimal if $`D_nM_n^2=P_nM(IP_n)MP_n`$ is small. The starting point is to find all of the eigenvalues $`\lambda _r`$ of $`M_n`$ and then to compute $`G(\lambda _r)`$ for each $`r`$. If $`G(\lambda _r)`$ is not sufficiently small for some $`r`$ then $`\lambda _r`$ is rejected as spurious. Let the remaining eigenvalues be listed in increasing order. The corresponding normalized eigenvectors $`g_r`$ of $`M_n`$ are used as starting vectors for an inverse power iteration to find $`f_r=f(\lambda _r)`$. As a check it should be verified that $`f_r`$ and $`g_r`$ are approximately equal, and that the set $`f_r`$ is approximately orthonormal as $`r`$ varies.

Now assume that one wants a list of values of $`F_n()`$ over a set of points covering the interval $`(\lambda _r,\lambda _{r+1})`$, in order to plot a graph. This is done separately on either side of the midpoint $`\gamma _r`$ of the interval, and we consider the left hand half. One chooses some subdivision $`\lambda _r=\sigma _1<\sigma _2<\mathrm{}<\sigma _m<\gamma _r`$ and then computes $`G(\sigma _r)`$ in order. At each stage the starting vector for the inverse iteration should be the eigenvector for the previous stage. In the right hand half interval one deals with the points in decreasing order.

The reason for using different procedures on either side of the midpoint becomes clear if one considers the exact case, in which $`D_n=M_n^2`$; the eigenvector is then constant in each subinterval and changes discontinuously as $`\lambda `$ increases through $`\gamma _r`$. In the inexact case one should expect $`f(\lambda )`$ to change rapidly near $`\lambda =\gamma _r`$, and it is therefore advisable to start the iteration with the vector $`f(\lambda ^{})+f(\lambda ^{\prime \prime })`$, where $`\lambda ^{}`$ and $`\lambda ^{\prime \prime }`$ are values already computed on either side of $`\gamma _r`$.

One also needs to evaluate $`F_n(\lambda )`$ for large values of $`n`$ at individual points $`\lambda `$ satisfying $`\lambda _r<\lambda <\lambda _{r+1}`$. Once again this may be done using inverse power iteration, but using $`f_r+f_{r+1}`$ as the initial vector.

We finally discuss how the computed values of $`F_n(\lambda )`$ are used to locate the spectrum of $`M`$. We use

$$\mathrm{Spec}(M)[\lambda F_n(\lambda ),\lambda +F_n(\lambda )]\mathrm{}$$

for all $`\lambda 𝐑`$, which follows from the spectral theorem, and the obvious

$$F_n(\lambda )F(\lambda )$$

for all $`\lambda 𝐑`$. Unless $`F_n(\lambda )`$ and $`F(\lambda )`$ are approximately equal, the above provide little information. One has two choices:

(i) One can supplement the computation of $`F_n(\lambda )`$ by prior information about the location of the spectrum of $`M`$.

(ii) One can examine the graph of $`F_n(\lambda )`$ to see whether its form is consistent with the assumption that it is approximately equal to $`F(\lambda )`$, using (8). If not then one can increase the value of $`n`$, or revert to (i).

Of course the former procedure is to be preferred, but the latter is likely to be easier. We illustrate both in the example discussed in the next section.

For the rest of the paper we make the following hypothesis.

(S) We have $`F_n(s)>F(s)`$ and $`|F_n^{}(s)|<1`$ for all $`s𝐑`$.

This can be ensured by replacing $`F_n`$ by

$$\widehat{F}_n(s)=F_n(s)+\epsilon \mathrm{e}^{F_n(s)}$$

for any $`\epsilon >0`$. Increasing $`F_n(s)`$ slightly makes the spectral bounds which we prove slightly weaker, but does not involve an essential loss of generality. We will use the following consequence of the hypothesis (S).

$$\mathrm{Spec}(M)(\lambda F_n(\lambda ),\lambda +F_n(\lambda ))\mathrm{}.$$

Our goal is to determine those eigenvalues of $`M`$ which lie in a given interval $`(\alpha ,\beta )`$; $`\alpha `$ and $`\beta `$ might sometimes be known points of the essential spectrum of $`M`$. Let $`\sigma _r`$, $`1rR`$ be an increasing sequence of numbers lying in the interval $`(\alpha ,\beta )`$; they are our first approximations to the eigenvalues of $`M`$ within the interval. They might be the non-spurious eigenvalues $`\lambda _r`$ of $`M_n`$ already described, or they might be the local minima of the function $`F_n(\lambda )`$, obtained by solving $`F_n^{}(\lambda )=0`$ using the bisection method. In the simplest cases these two methods of choosing $`\sigma _r`$ yield approximately the same values. If this fails one can still proceed, using the second choice of $`\sigma _r`$ just mentioned. Putting $`\mu _r=\sigma _rF_n(\sigma _r)`$ and $`\nu _r=\sigma _r+F_n(\sigma _r)`$ for $`1rR`$, we have

$$\mathrm{Spec}(M)(\mu _r,\nu _r)\mathrm{}$$

for such $`r`$. We also put $`\nu _0=\alpha `$ and $`\mu _{R+1}=\beta `$. Further progress depends upon the following hypothesis. Following (i) above, this might have been proved, possibly using a homotopy method. Or following (ii), it might simply be motivated by the numerical evidence, and tested by the consistency of the consequences drawn using it.

(H) We have

$$\nu _0<\mu _1<\nu _1<\mu _2<\mathrm{}<\mu _R<\nu _R<\mu _{R+1}.$$

The intervals $`(\nu _{r1},\mu _r)`$ do not meet $`\mathrm{Spec}(M)`$ for $`1rR+1`$, while $`[\mu _r,\nu _r]`$ each contain a single eigenvalue $`m_r`$ of $`M`$ for $`1rR`$.

Situations in which the operator $`M`$ has eigenvalues which are too close to be resolved using the subspace $`_n`$ can often be handled using a modified procedure, or more simply by increasing $`n`$.

We describe an iterative procedure for obtaining steadily more accurate enclosures of the eigenvalues. The following lemma provides the starting point for the iteration. The constructions are most easily understood by plotting the graph of a typical function $`F_n`$. See Figure 1 for an example in which $`\alpha =0`$, $`\beta =1`$ are eigenvalues of infinite multiplicity, and $`R=1`$.

###### Lemma 7

Assuming (S) and (H), we have

$$\mu _r<\mu _r^{}<m_r<\nu _r^{}<\nu _r$$

for $`1rR`$, where $`\nu _r^{}`$ and $`\mu _r^{}`$ are defined by solving the following equations.

$`F_n(s_r)`$ $`=`$ $`s_r\nu _{r1}`$
$`\nu _r^{}`$ $`=`$ $`s_r+F_n(s_r)`$
$`F_n(t_r)`$ $`=`$ $`\mu _{r+1}t_r`$
$`\mu _r^{}`$ $`=`$ $`t_rF_n(t_r)`$
Proof Let $`1rR`$. The function

$$\varphi (s)=F_n(s)s+\nu _{r1}$$

satisfies $`\varphi ^{}(s)=F_n^{}(s)1<0`$ for all $`s`$, $`\varphi (\nu _{r1})=F_n(\nu _{r1})>0`$ and

$$\varphi (\sigma _r)=F_n(\sigma _r)\sigma _r+\nu _{r1}=\nu _{r1}\mu _r<0.$$

Hence $`\nu _{r1}<s_r<\sigma _r`$. We next observe that

$$F_n(s_r)F_n(\sigma _r)<\sigma _rs_r$$

may be rewritten in the form $`\nu _r^{}<\nu _r`$. Also

$$\mathrm{Spec}(M)(s_rF_n(s_r),s_r+F_n(s_r))\mathrm{}$$

may be rewritten as

$$\mathrm{Spec}(M)(\nu _{r1},\nu _r^{})\mathrm{}.$$

Using (H) we now deduce that $`m_r<\nu _r^{}`$. The proof that $`\mu _r<\mu _r^{}<m_r`$ is similar.

Improving the above result iteratively depends upon setting up the inductive hypothesis properly. Our hypothesis below allows the possibility of increasing the value of $`n`$ from one stage of the iteration to the next, in other words of increasing the size of the test function space as the iteration proceeds. In practice one would choose a fairly small initial value of $`n`$ and apply the lemma repeatedly until the iteration stops yielding improvements, and then increase $`n`$ significantly and repeat the iteration with that value of $`n`$ until improvements cease.

###### Theorem 8

Let $`\mu _r`$, $`\mu _r^{}`$, $`\nu _r`$, $`\nu _r^{}`$, $`s_r`$ and $`t_r`$ satisfy

$$\mu _r\mu _r^{}<m_r<\nu _r^{}\nu _r$$

and

$`F_n(s_r)`$ $``$ $`s_r\nu _{r1}`$
$`\nu _r^{}`$ $``$ $`s_r+F_n(s_r)`$
$`F_n(t_r)`$ $``$ $`\mu _{r+1}t_r`$
$`\mu _r^{}`$ $``$ $`t_rF_n(t_r)`$
for $`1rR`$. Put $`\nu _0^{}=\sigma _0`$ and $`\mu _{R+1}^{}=\sigma _{R+1}`$. Given (S) and (H), define $`\mu _r^{\prime \prime }`$, $`\nu _r^{\prime \prime }`$, $`s_r^{}`$ and $`t_r^{}`$ for $`1rR`$ by

$`F_n(s_r^{})`$ $`=`$ $`s_r^{}\nu _{r1}^{}`$
$`\nu _r^{\prime \prime }`$ $`=`$ $`s_r^{}+F_n(s_r^{})`$
$`F_n(t_r^{})`$ $`=`$ $`\mu _{r+1}^{}t_r^{}`$
$`\mu _r^{\prime \prime }`$ $`=`$ $`t_r^{}F_n(t_r^{})`$
Then

$$\mu _r\mu _r^{}\mu _r^{\prime \prime }<m_r<\nu _r^{\prime \prime }\nu _r^{}\nu _r$$

for $`1rR`$.

Proof Let $`1rR`$. The function

$$\varphi (s)=F_n(s)s+\nu _{r1}^{}$$

satisfies $`\varphi ^{}(s)<0`$ for all $`s𝐑`$ and

$$\varphi (\nu _{r1}^{})=F_n(\nu _{r1}^{})>0.$$

Also

$$\varphi (s_r)=F_n(s_r)s_r+\nu _{r1}^{}\nu _{r1}^{}\nu _{r1}0.$$

Hence $`\nu _{r1}^{}<s_r^{}s_r`$. We next observe that

$$F_n(s_r^{})F_n(s_r)s_rs_r^{}$$

implies

$$\nu _r^{\prime \prime }=F_n(s_r^{})+s_r^{}F_n(s_r)+s_r\nu _r^{}.$$

Also

$$\mathrm{Spec}(M)(s_r^{}F_n(s_r^{}),s_r^{}+F_n(s_r^{}))\mathrm{}$$

is equivalent to

$$\mathrm{Spec}(M)(\nu _{r1}^{},\nu _r^{\prime \prime })\mathrm{}$$

and implies that $`m_r<\nu _r^{\prime \prime }`$. The proof that $`\mu _r^{}\mu _r^{\prime \prime }<m_r`$ is similar.

At each stage in the implementation of the iterative process one has further consistency checks, namely the numerical procedure need not yield the inequality $`\mu _r^{\prime \prime }<\nu _r^{\prime \prime }`$ for all $`r`$. If this does not happen then (H) must be false; that is the operator $`M`$ must have more spectrum than assumed.

Plotting the graph of a typical function $`F_n`$ suggests certain results about the location of $`s_r`$ and $`t_r`$; these are proved below.

###### Lemma 9

Assume (S) and (9). Given $`\lambda 𝐑`$, let $`\lambda ^{}`$ denote the smallest point in $`\mathrm{Spec}(M)`$ greater than $`\lambda `$. Assume $`\lambda ^{}>\lambda +4\epsilon `$ and let $`s`$ be the solution of $`F_n(s)=s\lambda `$. Then

$$\frac{\lambda +\lambda ^{}}{2}<s<\frac{\lambda +\lambda ^{}}{2}+2\epsilon .$$

Proof We combine the fact that the function $`\varphi (s)=F_n(s)s+\lambda `$ is strictly monotone decreasing and continuous with the inequalities

$$\varphi ((\lambda +\lambda ^{})/2)>F((\lambda +\lambda ^{})/2)(\lambda +\lambda ^{})/2+\lambda =0$$

and

$$\varphi ((\lambda +\lambda ^{})/2+2\epsilon )<F((\lambda +\lambda ^{})/2+2\epsilon )+2\epsilon (\lambda +\lambda ^{})/22\epsilon +\lambda =2\epsilon <0.$$

###### Theorem 10

Suppose that $`M`$ has a single eigenvalue $`m`$ in $`(\alpha ,\beta )`$. Define $`s_n`$ and $`t_n`$ by

$`F_n(s_n)`$ $`=`$ $`s_n\alpha `$
$`\nu _n`$ $`=`$ $`s_n+F_n(s_n)(=2s_n\alpha )`$
$`F_n(t_n)`$ $`=`$ $`\beta t_n`$
$`\mu _n`$ $`=`$ $`t_nF_n(t_n)(=2t_n\beta ).`$
Then

$$\mu _n<m<\nu _n$$

for all $`n1`$. Moreover

$$\underset{n\mathrm{}}{lim}s_n=(\alpha +m)/2,\underset{n\mathrm{}}{lim}t_n=(m+\beta )/2$$
(14)
and

$$\underset{n\mathrm{}}{lim}\mu _n=\underset{n\mathrm{}}{lim}\nu _n=m.$$

Proof The first statement follows from Lemma 7 with $`R=1`$, $`\nu _0=\alpha `$ and $`\mu _2=\beta `$. The other statements are then consequences of Lemma 9.

## 6 The Zimmermann-Mertins Paper

We have repeatedly referred to solving the equations $`F_n(s)=s\nu `$ and $`F_n(t)=\mu t`$ for various $`\mu ,\nu `$. This may be carried out by bisection or by Newton’s method. By examining the paper of Zimmermann and Mertins, , one discovers a method which is probably better in most circumstances; see (21) below.

We assume that $`(\nu ,\mu )`$ is an interval which contains a single eigenvalue $`m`$ of $`M`$, but no other point of Spec $`(M)`$. E.g. in the situation described in hypothesis (H), we have in mind tho choose $`\nu =\nu _{r1},\mu =\mu _{r+1},m=m_r`$, for any $`r\{1,\mathrm{},R\}`$. Suppose for simplicity that $`\nu `$ and $`\mu `$ are not eigenvalues of $`M`$.

Let $`_n`$ denote an $`n`$-dimensional trial function space. To ensure that the exact eigenfunction corresponding to $`m`$ is ‘not too far’ from $`_n`$, we make the additional assumption

$`\nu <\underset{v_n\{0\}}{\mathrm{max}}{\displaystyle \frac{Mv,v}{v,v}},\underset{v_n\{0\}}{\mathrm{min}}{\displaystyle \frac{Mv,v}{v,v}}<\mu .`$ (15)
(Without (15), the Lehmann method gives nothing, and the $`F_n`$-method gives a lower bound for $`m`$ less than $`\nu `$, resp. an upper bound larger than $`\mu `$, which makes it more or less useless.)

The right-definite Lehmann method may be applied to the situation in which only one eigenvalue is to be bounded. We follow \[20, Theorem 4.2\], where the spectral parameter $`\rho `$ is once chosen to be $`\mu `$, and once $`\nu `$. Let $`\{v_1,\mathrm{},v_n\}`$ denote a basis of $`_n`$, and let $`\tau ^{}`$ be the smallest eigenvalue of the $`n\times n`$ matrix eigenvalue problem

$`(Mv_i\mu v_i,v_k)_{i,k}x=\tau (Mv_i\mu v_i,Mv_k\mu v_k)_{i,k}x`$ (16)
and $`\tau ^+`$ the largest eigenvalue of

$`(Mv_i\nu v_i,v_k)_{i,k}x=\tau (Mv_i\nu v_i,Mv_k\nu v_k)_{i,k}x.`$ (17)
Then

$`\mu +{\displaystyle \frac{1}{\tau ^{}}}m\nu +{\displaystyle \frac{1}{\tau ^+}}.`$ (18)
(Note that $`\tau ^{}<0<\tau ^+`$, according to (15)).

The following theorem states that the $`F_n`$-method gives exactly the right-definite Temple-Lehmann-Maehly bounds, in the Zimmermann-Mertins setting. See also \[7, Theorem 6\].

###### Theorem 11

Let $`\tau ^{}`$ and $`\tau ^+`$ as above, and let $`s`$ and $`t`$ denote the solutions of (cf. Lemma 7):

$`F_n(s)=s\nu ,F_n(t)=\mu t.`$ (19)
Then,

$`s+F_n(s)=\nu +{\displaystyle \frac{1}{\tau ^+}},tF_n(t)=\mu +{\displaystyle \frac{1}{\tau ^{}}}.`$ (20)
Hence

$`s=\nu +{\displaystyle \frac{1}{2\tau ^+}},t=\mu +{\displaystyle \frac{1}{2\tau ^{}}}.`$ (21)
Proof For $`\rho 𝐑`$, define

$`A_2(\rho )`$ $`:=`$ $`(Mv_i\rho v_i,Mv_k\rho v_k)_{i,k=1,\mathrm{},n},`$
$`A_1(\rho )`$ $`:=`$ $`(Mv_i\rho v_i,v_k)_{i,k=1,\mathrm{},n},A_0:=(v_i,v_k)_{i,k=1,\mathrm{},n}.`$
Straightforward calculations and use of (19) give, for all $`x𝐑^n`$,

$`x^TA_2(s)xF_n(s)^2x^TA_0x`$ $`=`$ $`x^TA_2(\nu )x+2(\nu s)x^TA_1(\nu )x`$ (22)
$`+[(\nu s)^2F_n(s)^2]x^TA_0x`$
$`=`$ $`x^TA_2(\nu )x2F_n(s)x^TA_1(\nu )x.`$
By the definition of $`F_n(s)`$, the left-hand side of (22) is non-negative for all $`x`$, and zero for some specific $`x`$. This is therefore true also for the right-hand side, i.e.,

$`{\displaystyle \frac{x^TA_1(\nu )x}{x^TA_2(\nu )x}}{\displaystyle \frac{1}{2F_n(s)}}\mathrm{for}\mathrm{all}x𝐑^n\{0\},`$
with equality for some specific $`x`$. Thus, $`1/(2F_n(s))`$ is the maximum of the Rayleigh quotient of problem (17), i.e.,

$`{\displaystyle \frac{1}{2F_n(s)}}=\tau ^+.`$ (23)
Using (19) again, we obtain

$`s+F_n(s)=\nu +2F_n(s)=\nu +{\displaystyle \frac{1}{\tau ^+}}.`$
A similar method yields

$`x^TA_2(t)xF_n(t)^2x^TA_0x=x^TA_2(\mu )x+2F_n(t)x^TA_1(\mu )x,`$
and the same arguments as before yield

$`{\displaystyle \frac{x^TA_1(\mu )x}{x^TA_2(\mu )x}}{\displaystyle \frac{1}{2F_n(t)}}\mathrm{for}\mathrm{all}x𝐑^n\{0\},`$
with equality for some specific $`x`$, which implies

$`{\displaystyle \frac{1}{2F_n(t)}}=\tau ^{}`$ (24)
and thus

$`tF_n(t)=\mu 2F_n(t)=\mu +{\displaystyle \frac{1}{\tau ^{}}}.`$
The final statement of the theorem is obtaining by combining (19), (23) and (24).

## 7 An Example in One Dimension

Suppose that the operator $`H`$ has a complete orthonormal set of eigenfunctions $`\{\varphi _r\}_{r=1}^{\mathrm{}}`$, whose corresponding eigenvalues $`\lambda _r`$ are all positive. Suppose also that the eigenvalues are listed in increasing order and tend to infinity. If we put

$$\psi _r=\lambda _r^{1/2}D\varphi _r$$
(25)
then $`\{\psi _r\}_{r=1}^{\mathrm{}}`$ forms an orthonormal basis in $``$ and the matrix elements of $`M`$ with respect to this basis are

$$M_{r,s}=\lambda _r^{1/2}\lambda _s^{1/2}Q_K(\varphi _r,\varphi _s)$$

The truncation method of computing the eigenvalues consists of choosing $`n`$ large enough and then finding the eigenvalues of the matrix $`M_n`$ obtained from $`M`$ by considering only those coefficients for which $`r,sn`$. This is equivalent to finding the spectrum of $`M_n=P_nMP_n`$, where $`P_n`$ is the orthogonal projection onto $`_n=\mathrm{lin}\{\psi _1,\mathrm{},\psi _n\}`$. It is well known that every point of $`\mathrm{Spec}(M)`$ is the limit of eigenvalues of $`M_n`$ as $`n\mathrm{}`$. Spectral pollution arises from the fact that the converse is not true.

As an example we consider the operator $`H`$ acting in $`L^2(0,\pi )`$ according to the formula

$$Hf(x)=\frac{\mathrm{d}^2f}{\mathrm{d}x^2}$$

and subject to Dirichlet boundary conditions. We then have

$$\varphi _n(x)=\sqrt{2/\pi }\mathrm{sin}(nx)$$

and $`\lambda _n=n^2`$ for $`n=1,2,\mathrm{}`$. We take $`K`$ to be the differential operator

$$Kf(x)=\frac{\mathrm{d}}{\mathrm{d}x}\left\{b(x)\frac{\mathrm{d}f}{\mathrm{d}x}\right\}$$

where $`b()`$ is any non-negative bounded function on $`(0,\pi )`$. If $`D`$ is the differentiation operator with domain $`W_0^{1,2}(0,\pi )`$ then $`H=D^{}D`$ and $`K=D^{}BD`$, where $`B`$ is the operator of multiplication by $`b()`$ acting in $`L^2(0,\pi )`$. The range of $`D`$ is not $`L^2(0,\pi )`$ but the subspace $``$ of all functions $`f`$ such that

$$_0^\pi f(x)dx=0.$$

This corresponds to the fact that in the present context (25) translates to

$$\psi _n(x)=\sqrt{2/\pi }\mathrm{cos}(nx).$$

Since $`n=1,2,\mathrm{}`$ these do not form a complete orthonormal set in $`L^2(0,\pi )`$.

The bounded operator $`M`$ is given for $`f`$ by

$$Mf=bfbf,\psi _0\psi _0$$

where $`\psi _0(x)=\pi ^{1/2}`$. It follows by trace class perturbation techniques that $`M`$ and $`B`$ have the same essential spectrum. This equals the closure of the set of essential values of $`b`$, or the closure of its actual set of values if it is piecewise continuous. However $`M`$ may also have a further eigenvalue in any gap of its essential spectrum.

###### Lemma 12

The isolated eigenvalues of $`M`$ are of multiplicity $`1`$ and are the nonzero solutions of the equation

$$\mathrm{\Xi }(\lambda ):=_0^\pi \frac{b(x)}{b(x)\lambda }dx=\pi .$$

In particular if $`b`$ is the characteristic function of the measurable set $`E(0,\pi )`$ then $`M`$ has only one isolated eigenvalue, namely $`\lambda =(\pi |E|)/\pi `$.

Proof The task is to find all $`fL^2(0,\pi )`$ which satisfy $`f,\psi _0=0`$ and

$$bfbf,\psi _0\psi _0=\lambda f$$
(26)
where $`\lambda `$ is not in the essential range of $`b`$. We may normalize $`f`$ by $`bf,\psi _0=1`$ since if the LHS equals $`0`$, we see from (26) that $`\lambda `$ lies in the essential range of $`b`$. We then insert

$$f=(b\lambda )^1\psi _0$$

into the normalization equation. The function $`\mathrm{\Xi }`$ is analytic and has a positive derivative in each gap in the spectrum of $`M`$, so there is at most one eigenvalue in each gap.

The entries of the matrix $`A_n=P_nMP_n`$ of (11) are given for $`1r,sn`$ by

$$A_{n,r,s}=\frac{2}{\pi }_0^\pi b(x)\mathrm{cos}(rx)\mathrm{cos}(sx)dx.$$

The matrix entries of $`D_n=P_nM^{}MP_n`$ are calculated using the formula

$$M^{}Mf,g=_0^\pi b(x)^2f(x)\overline{g(x)}dxbf,\psi _0\psi _0,bg,$$

valid for all $`f,g`$.

We consider the exactly soluble example

$$b(x)=\{\begin{array}{cc}1\hfill & \text{if }0x\alpha \hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$

where $`0<\alpha <\pi `$ is arbitrary. The operator $`M`$ then has two-point essential spectrum $`\{0,1\}`$ together with an eigenvalue $`m=(\pi \alpha )/\pi `$ of multiplicity $`1`$. If one puts $`\alpha =\pi /2`$ then $`m=1/2`$; numerical calculations show that as $`n`$ increases $`A_n`$ has an increasing number of eigenvalues moving inwards from $`0`$ and $`1`$ and apparently converging slowly to $`0.5`$. The following table lists the four eigenvalues of $`A_n`$ closest to $`0.5`$.

$$\begin{array}{ccccc}n& \lambda _1& \lambda _2& \lambda _3& \lambda _4\\ 50& 0.0526& 0.4632& 0.5368& 0.9474\\ 100& 0.0691& 0.4750& 0.5250& 0.9309\\ 200& 0.0856& 0.4830& 0.5170& 0.9144\\ 400& 0.1019& 0.4884& 0.5116& 0.8981\end{array}$$

The use of the techniques of this paper not only avoids the problem of spectral pollution, but also provides accurate results for much smaller values of $`n`$. We start by putting $`\alpha =\pi /2`$ and $`n=10`$. The function $`F_{10}(\lambda )`$ has three local minima. They are given in the following table, in which the numbers are all accurate to eight decimal places.

$$\begin{array}{cc}\lambda & F_{10}(\lambda )\\ 0.00& 0.00000679\\ 0.50& 0.10049280\\ 1.00& 0.00000679\end{array}$$

The central minimum is exactly equal to $`1/2`$ as a result of the unusual symmetry for $`\alpha =\pi /2`$. The hypothesis (H) is therefore valid with $`\nu _0=0`$, $`\mu _1=0.39950720`$, $`\nu _1=0.60049280`$ and $`\mu _2=1`$. An application of Lemma 7 now yields $`s_1=0.25`$ and $`t_1=0.75`$ to six decimal places, which confirms that the eigenvalue equals $`0.5`$ to the same accuracy.

We repeat the calculation for $`\alpha =1.0`$, in order to see whether the unexpected accuracy of the results are caused by the symmetry. Two representative graphs are shown in Figure 1. The three local minima of $`F_8(\lambda )`$ are now as follows, accurate to six decimal places.

$$\begin{array}{cc}\lambda & F_8(\lambda )\\ 0& 0\\ 0.697669& 0.106699\\ 0.999313& 0.026051\end{array}$$

The hypothesis (H) is therefore valid with $`\nu _0=0`$, $`\mu _1=0.590970`$, $`\nu _1=0.804368`$ and $`\mu _2=1`$. An application of Lemma 7 now yields $`s_1=0.340845`$ and $`t_1=0.840844`$ to six decimal places. This implies that

$$0.681688\mu 0.681690.$$

The accuracy of the upper and lower bounds on the eigenvalue $`\mu `$ depends on the rate of convergence in (14), which is closely associated to the rate of convergence in

$$\underset{n\mathrm{}}{lim}F_n(s)=F(s),\underset{n\mathrm{}}{lim}F_n(t)=F(t)$$

where $`s=m/2`$ and $`t=(m+1)/2`$. For reasons which we do not understand the rate of convergence at these two points is much faster than at other points in $`(0,1)`$. This phenomenon is illustrated in the following table, in which $`u`$ and $`v`$ are ‘general’ points.

Putting $`\alpha =1.0`$ we have $`m=0.681690`$, $`s=0.340845`$ and $`t=0.840845`$. We also put $`u=0.5`$ and $`v=0.75`$.

$$\begin{array}{ccccc}n& F_n(s)& F_n(t)& F_n(u)& F_n(v)\\ 2& 0.340873& 0.269294& 0.254218& 0.433012\\ 4& 0.340845& 0.162575& 0.238589& 0.146875\\ 6& 0.340845& 0.159219& 0.221431& 0.137537\\ 8& 0.340845& 0.159156& 0.210113& 0.116463\\ 10& 0.340845& 0.159155& 0.207660& 0.108933\\ 100& 0.340845& 0.159155& 0.187173& 0.077840\\ \mathrm{}& 0.340845& 0.159155& 0.181690& 0.068310\end{array}$$

The final row lists the values of $`F()`$.

The surprisingly rapid convergence of $`F_n(s)`$ and $`F_n(t)`$ does not occur only for the two values of $`\alpha `$ presented. A further computation indicates that if we put $`s=m/2`$ then

$$s=F(s)F_{10}(s)F(s)+10^6$$

for all $`0<\alpha <2.3`$. For larger $`\alpha `$ the difference $`F_{10}(s)F(s)`$ increases steadily. We also have

$$s=F(s)F_{50}(s)F(s)+10^8$$

for all $`0<\alpha <2.9`$.

## 8 An Example with Variable Coefficients

It might be thought that the cause of the spectral pollution in the above examples was the fact that the basis chosen was badly adapted to the operator.

In this section we consider the same type of example as in the last, but using a different type of basis, constructed using the Krylov subspace method. This example does not display spectral pollution, but this is not because the associated matrix is tridiagonal. A recent example of Denisov and Simon shows that spectral pollution can occur in the following sense even for tridiagonal (Jacobi) matrices, . They construct a bounded, self-adjoint, tridiagonal operator $`A`$ acting on $`l^2(𝐙^+)`$ with spectrum $`[5,1][1,5]`$. The existence of a gap in its spectrum is not a priori obvious. For certain values of $`n`$, the truncation $`A_n`$ of $`A`$ to the space $`_n`$ of sequences with support in $`\{1,n\}`$ has a single eigenvalue in $`(1,1)`$. The set of limit points of these eigenvalues is equal to the whole interval $`[1,1]`$.

We replace the interval $`(0,\pi )`$ by $`(2,2)`$ and put

$$b(x)=\{\begin{array}{cc}\alpha ^{}+\beta ^{}x\hfill & \text{ if }x<0\hfill \\ \alpha ^++\beta ^+x\hfill & \text{ if }x>0.\hfill \end{array}$$

The essential spectrum of $`b`$, and of $`M`$, is thus $`[\alpha ^{}2\beta ^{},\alpha ^{}][\alpha ^+,\alpha ^++2\beta ^+]`$. Assuming that $`\beta ^\pm >0`$ and $`\alpha ^{}<\alpha ^+`$, there is an eigenvalue $`\lambda `$ in the spectral gap $`(\alpha ^{},\alpha ^+)`$. This is evaluated by solving

$$\mathrm{\Xi }(\lambda ):=_2^2\frac{b(x)}{b(x)\lambda }dx=4.$$
(27)
as in Lemma 12.

We next write down the basis which we will use. Let $`\{𝒫_n\}_{n=1}^{\mathrm{}}`$ denote the Legendre polynomials, and define

$$Q_n(x)=\sqrt{\frac{2n+1}{2}}𝒫_n(x)$$

so that

$$xQ_n(x)=\gamma _{n+1}Q_{n+1}(x)+\gamma _nQ_{n1}(x)$$

for all $`n0`$, where

$$\gamma _n=\frac{n}{\sqrt{4n^21}}.$$

Then define $`Q_n^\pm `$ on $`(2,2)`$ by

$`Q_n^+(x)=\{\begin{array}{cc}Q_n(x1)\hfill & \text{ if }x>0\hfill \\ 0\hfill & \text{ if }x<0,\hfill \end{array}`$
$`Q_n^{}(x)=\{\begin{array}{cc}0\hfill & \text{ if }x>0\hfill \\ Q_n(x+1)\hfill & \text{ if }x<0.\hfill \end{array}`$
Finally define $`\psi _n`$ for all $`n𝐙`$ by

$$\psi _n=\{\begin{array}{cc}Q_n^+\hfill & \text{ if }n1\hfill \\ Q_{|n|}^{}\hfill & \text{ if }n1\hfill \\ (Q_0^++Q_0^{})/\sqrt{2}\hfill & \text{ if }n=0.\hfill \end{array}$$

It is easy to check that $`\{\psi _n\}_{n𝐙}`$ is a complete orthonormal set in

$$=\{fL^2(2,2):_2^2f(x)dx=0\}.$$

Let $`P_n`$ be the orthogonal projection of $``$ onto $`_n=\mathrm{lin}\{\psi _n,\mathrm{},\psi _n\}`$. Then it follows from standard properties of the Legendre polynomials that $`M(_n)_{n+1}`$. Hence $`M`$ has a tridiagonal matrix. Routine calculations yield

$$M_{r,r}=\{\begin{array}{cc}\alpha ^++\beta ^+\hfill & \text{ if }r1\hfill \\ \alpha ^{}\beta ^{}\hfill & \text{ if }r1\hfill \\ (\alpha ^++\beta ^++\alpha ^{}\beta ^{})/2\hfill & \text{ if }r=0.\hfill \end{array}$$

Also

$$M_{r,r+1}=M_{r+1,r}=\{\begin{array}{cc}\beta ^+\gamma _{r+1}\hfill & \text{ if }r1\hfill \\ \beta ^+\gamma _1/\sqrt{2}\hfill & \text{ if }r=0\hfill \\ \beta ^{}\gamma _1/\sqrt{2}\hfill & \text{ if }r=1\hfill \\ \beta ^{}\gamma _{|r|}\hfill & \text{ if }r2.\hfill \end{array}$$

All other coefficients of $`M`$ vanish. We observe that $`M_{r,r}`$ and $`M_{r,r+1}`$ have different limits as $`r\pm \mathrm{}`$. Operators whose coefficients have different asymptotics on the left and right have an interesting spectral and scattering theory, studied, for example, in .

We finally compute $`B_n(\lambda )`$ from (11), obtaining the formula

$$B_n(\lambda )=(A_n\lambda I)^2+E_n$$

where $`E_n=P_nM(IP_n)MP_n`$. Explicitly we have

$$E_{n,r,s}=\{\begin{array}{cc}(\beta ^{}\gamma _{n+1})^2\hfill & \text{ if }r=s=n\hfill \\ (\beta ^+\gamma _{n+1})^2\hfill & \text{ if }r=s=n\hfill \\ 0\hfill & \text{ otherwise.}\hfill \end{array}$$

We carried out numerical calculations for the case $`\alpha ^+=1`$, $`\alpha ^{}=1`$, $`\beta ^+=40`$, $`\beta ^{}=30`$, $`n=50`$, so that $`dim(_n)=2n+1`$. The essential spectrum of $`M`$ is $`[61,1][1,81]`$ and by solving (27) we find that $`M`$ also has a single eigenvalue $`\mu 0.453261434`$. Plotting the function $`F(\lambda )`$ on $`[0,1]`$ shows that it has a well-defined local minimum near this point and nowhere else in the interval. From the values

$`F_{50}(0.23)`$ $`=`$ $`0.22326196`$
$`F_{50}(0.72)`$ $`=`$ $`0.26673868`$
we obtain the enclosures

$$0.72F_{50}(0.72)\mu 0.23+F_{50}(0.23)$$

or, more explicitly $`\mu =0.45326_{13}^{20}`$. Of course all these computations could be carried out using interval arithmetic, but we did not do so. These upper and lower bounds on $`\lambda `$ compare with the completely uncontrolled values for $`\lambda `$ obtained by the truncation method.

$$\begin{array}{cc}n& \mu \\ 10& 0.50275311\\ 20& 0.45504606\\ 50& 0.45326153\\ 75& 0.45326143\\ 100& 0.45326143\end{array}$$

In this example we have found that by evaluating $`F_{50}(\lambda )`$ at only two points, we obtained upper and lower bounds on $`\mu `$ which were of comparable accuracy to the uncontrolled value obtained by truncation. The full plot of $`F_n(\lambda )`$ gives valuable geometrical insights, but may be carried out for a smaller value of $`n`$ than is used for the final bounds; this is, of course, much faster.

Acknowledgments We would like to thank M Levitin and E Shargorodsky for valuable discussions.