A finite element computation of eigenvalues of elliptic operators on compact manifolds

Stefan Bechtluft-Sachs

September 27, 2004

Abstract

A procedure to numerically compute the small eigenvalues of a first order self adjoint elliptic operator acting on sections of a Hermitian vector bundle over a compact Riemannian manifild. Explicit error bounds for the piecewise linear finite elements are derived.

1 Introduction

On a compact manifold M carrying a complex vector bundle E we consider an elliptic first order partial differential operator P with smooth coefficients. With the objective to compute its small eigenvalues, we approximate P by its restriction P^{□} to a suitable finite dimensional subspace V ⊂ L^{2}(E). We will describe the case where the finite elements v ∈ V are piecewise linear with respect to a given triangulation |K| = M and a bundle embedding E ⊂ M × ℂ^{L} of the coefficient bundle E into the trivial bundle of rank L.

This difficulty in defining the finite elements does not arise in the case of elliptic operators on domains in ℝ^{n}. The error estimates for the Laplace operator on domains in ℝ^{n} have been derived in [L] for spline approximations of arbitrary order. Using the approximation results of [N] an extension of our estimates to higher order splines should be straightforward, but will not be given here.

The main objective of this paper is the discretisation error, always assuming that the eigenvalues of the finite dimensional approximation P^{□} are computed exactly. On the other hand the procedure works for any selfadjoint elliptic differential operator. In particular we do not assume that the spectrum of P be bounded from below or even positive. Therefore the procedure can be applied to any geometric operator such as (twisted) Dirac operators on a Riemannian manifold.

The outline of the paper is as follows. In section 2 we formulate our main Theorem 1 and explain how it can be used to approximately compute the eigenvalues of P in a given interval [−Λ, Λ]. In section 3 we define the piecewise linear finite elements for a vector bundle over a compact manifold. In section 4 we derive the explicit formulae for the error estimates in Theorem 1. These depend on pointwise estimates for a (local) parametrix for P and its remainder which we derive in section 5 in a form adapted to our purpose. The existence of such estimates follows from the Sobolev inequality and elliptic regularity but we need the explicit expressions for the constants.

2 Computation of Small Eigenvalues

We choose a Riemannian metric on M and endow E ⊂ M × ℂ^{L} with a Hermitian metric by restricting the standard Hermitian metric of M × ℂ^{L}. Usually we denote by | · | the pointwise norm, by ∥ · ∥ the L^{2}-norm and by ∥ · ∥_{∞} the supremum of | · |. The computation of the small eigenvalues of P hinges on the following

Theorem 1 Assume that f ∈ C^{∞}(E) is a unit eigenvector of P with eigenvalue λ, i.e. P f = λf, ∥f∥ = 1. Then there is v ∈ V satisfying pointwise estimates

∥f − v∥_{∞} ≤ δ

and

∥P v − λv∥_{∞} ≤ ϵ + |λ|δ .

The values of δ = δ(λ) and ϵ = ϵ(λ) are explicitely given by (4.2) and (4.3) in section 4 and increase monotoneously with λ.

In particular one finds an almost eigenvector of P^{□}, i.e. v ∈ V such that ∥v∥ ≥ 1 − δ(vol(M))^{1/2} and

∥P^{□}v − λv∥^{2} ≤ (ϵ + |λ|δ)^{2}vol(M) .  (2.1)

Note that one finds eigenvectors v of P^{□} if P is self-adjoint because P^{□} then is symmetric. This follows by expanding an almost eigenvector v in a basis of eigenvectors of P^{□} as in the argument at the end of this section.

The apriori bounds δ and ϵ will be deduced from pointwise estimates for f and its derivative df in section 4, which in turn follow from the estimates of the second derivative d^{2}f in section 5. These latter estimates are independent on the triangulation, whereas the estimates for f and df improve under subdivision provided that the n-simplices of the triangulation do not degenerate. In this case one gets values of ϵ and δ tending to 0.

Suppose we wanted to compute the eigenvalues λ ∈ [−Λ, Λ] of a self adjoint elliptic first order differential operator P acting on sections of a Hermitian vectorbundle E over a compact Riemannian manifold M. Recall that P has discrete spectrum and that the eigenvectors are smooth by elliptic regularity. We embed E isometrically in a trivial bundle M × ℂ^{L}. Relying on the above theorem we exclude eigenvalues of P by computing eigenvalues of the finite dimensional operator P^{□}. Here we can work with δ = δ(Λ) and ϵ = ϵ(Λ) in (2.1).

Conversely we can show existence of an eigenvalue in a certain interval once we have found a unit eigenvector v^{□} ∈ V with eigenvalue λ^{□} of P^{□}. To that end we first compute ∥P v − λ^{□}v∥ =: α. Let {f_{i}}_{i∈ℕ} be an orthonormal basis of L^{2}E consisting of eigenvectors of P and such that P f_{i} = λ_{i}f_{i}, Spec(P) = {λ_{i} | i ∈ ℕ}. We expand

v^{□} = ∑_{i} a_{i}f_{i} , ∑_{i} |a_{i}|^{2} = 1

and compute

P v^{□} − λ^{□}v^{□} = ∑_{i} a_{i}(λ_{i} − λ^{□})f_{i} .

In particular

α^{2} ≥ ∑_{i} |a_{i}|^{2}|λ_{i} − λ^{□}|^{2} ≥ min_{i}{|λ_{i} − λ^{□}|^{2}}

and there is λ_{i} ∈ Spec(P) with |λ_{i} − λ^{□}| ≤ α.

3 The Finite Elements

For a simplicial complex K we denote by K_{n} the set of its n-simplices σ = (σ_{0}, …, σ_{n}) and by

Δ^{n} = {(t_{0}, …, t_{n}) ∈ ℝ^{n+1} | ∑_{i=0}^{n} t_{i} = 1, t_{i} ≥ 0}

the standard n-simplex. Let M be given as the geometric realization of a simplicial complex K (plus a smooth structure), i.e.

M = |K| = ⋃_{σ∈K_{n}} σ × Δ^{n} /∼  (3.1)

with the identifications

σ × (t_{0}, …, t_{i−1}, 0, t_{i+1}, …, t_{n}) ∼ σ′ × (t_{0}, …, t_{i−1}, 0, t_{i+1}, …, t_{n})

if σ = (σ_{0}, …, σ_{n}), σ′ = (σ′_{0}, …, σ′_{n}) ∈ K_{n} with σ_{j} = σ′_{j} for all j ≠ i. Let E ⊂ M × ℂ^{L} be a subbundle and denote by π: M × ℂ^{L} → E the Hermitian projection. The finite elements we consider are the projection to E of piecewise linear sections of M × ℂ^{L}. Thus

V′ := {v: M → ℂ^{L} | v(x) ∈ E if x ∈ K_{0} and

v((x_{0}, …, x_{n}) × (t_{0}, …, t_{n})) = ∑_{i=0}^{n} t_{i}v(x_{i}) for (x_{0}, …, x_{n}) ∈ K_{n}} ,

V := {π ∘ v | v ∈ V′} .

The dimension of these spaces is N times the cardinality of K_{0}. Since V is contained in the Sobolev space H_{1}(E) of sections of E, we can compute the L^{2} scalar product

⟨P v, w⟩ = ∫_{M}(P v, w)dvol_{g}  (3.2)

for v, w ∈ V where (·, ·) denotes Hermitian scalar product on E and dvol_{g} the measure corrsponding to Riemannian metric on M. We define the approximate operator P^{□}: V → V by (3.2), i.e. as P^{□} := pr_{V} P |_{V}, where pr_{V} denotes the Hermitian projection L^{2}(E) → V.

4 Error Estimates

An eigenvector f: M → E, P f = λf, splits in f = v +h with v ∈ V, v = πv′, v′ ∈ V′, such that v(p) = f(p) for all p ∈ K_{0} ⊂ M. We have

|P v − λv| = |λh − P h| ≤ |λ||h| + |P h| .

In the sequel we will estimate |h| and |P h| pointwise over an n-simplex σ × Δ^{n} ⊂ M.

We fix an open covering of M by charts Φ_{s}: U_{s} → ℝ^{n} covered by bundle charts Φ̂_{s}: E|_{U_{s}} → ℝ^{n} × ℂ^{N} and a function s(σ) sucht that every n-simplex σ × Δ^{n} is contained in U_{s(σ)}. From the triangulation we have maps j_{σ}: Δ^{n} → U_{s(σ)} covered by bundle maps ĵ_{σ}: Δ^{n} × ℂ^{N} → E|_{U_{s(σ)}}. Denote by Φ_{σ}, Φ̂_{σ} the compositions Φ_{σ} := Φ_{s(σ)} ∘ j_{σ}: Δ^{n} → ℝ^{n} and Φ̂_{σ} := Φ̂_{s(σ)} ∘ ĵ_{σ}: Δ^{n} × ℂ^{N} → ℝ^{n} × ℂ^{N}. In the following diagam Δ^{n}, ℝ^{n}, ℂ^{N} and ℂ^{L} carry the standard metrics.

Δ^{n} ⟶^{f} Δ^{n} × ℂ^{N} ⟶^{ω} Δ^{n} × ℂ^{L}
↓j_{σ}  ↓ĵ_{σ}  ↓j_{σ}×id
U_{s}  ⟶^{f} E|_{U_{s}}  ↪ U_{s} × ℂ^{L}
↓Φ_{s}  ↓Φ̂_{s}  ↓Φ_{s}×id
ℝ^{n} \underrightarrow{f̃} ℝ^{n} × ℂ^{N} \underrightarrow{ω̃} ℝ^{n} × ℂ^{L}

Slightly abusing notation we will write f(x) = (x, f(x)), ω(x, v) = (x, ω_{x}v) and analagously for f̃ and ω̃. From section 5, (5.4) we obtain pointwise apriori estimates

|f̃| ≤ C̃_{0} |df̃| ≤ C̃_{1} |d^{2}f̃| ≤ C̃_{2}

independent of the triangulation. Hence

|ωf(x)| ≤ ∥ω̃∥_{∞}C̃_{0} =: C_{0}
|d_{x}ωf| ≤ ∥dω̃∥_{∞}∥dΦ_{σ}∥_{∞}C̃_{0} + ∥ω̃∥_{∞}∥dΦ_{σ}∥_{∞}C̃_{1} =: C_{1}
|d_{x}^{2}ωf| ≤ (∥d^{2}ω̃∥_{∞}∥dΦ_{σ}∥2_{∞} + ∥dω̃∥_{∞}∥d^{2}Φ_{σ}∥_{∞}) C̃_{0}

+ (2∥dω̃∥_{∞}∥dΦ_{σ}∥2_{∞} + ∥ω̃∥_{∞}∥d^{2}Φ_{σ}∥_{∞}) C̃_{1}

+ ∥ω̃∥_{∞}∥dΦ_{σ}∥2_{∞}C̃_{2}

=: C_{2}

where ∥·∥_{∞} denotes the supremum of the respective fibrewise operator norm, e.g. ∥ω̃∥_{∞} := sup_{x} ∥ω̃(x)∥_{op}.

Note that C_{0}, C_{1} and C_{2} depend on the triangulation. Passing to a subdivision replaces the j_{σ} by the composition with affine linear maps Δ^{n} → Δ^{n}. Choosing the subdivision appropriately, e.g. barycentric, these affine linear maps will have differential ≤ α < 1. The constants C′_{i}, i = 0, 1, 2 for the subdivision then satisfy C′_{i} ≤ α^{i}C_{i}. For instance, the barycentric subdivision has α = \sqrt{n/2(n + 1)}.

Next we estimate h′ = ωf − v′ and h = πh′. Since d^{2}v′ = 0, d^{2}h′ = d^{2}ωf the above estimates yield |d^{2}h′| ≤ C_{2}. From the definition of v′ we also have

h′(q_{l}) = 0 for the vertices q_{l} = (0, …, 0, 1^{l}, 0, …, 0) ∈ ℝ^{n+1} of Δ^{n}.

For z ∈ ℂ^{L} = ℝ^{2L} let h′_{z}(x) := h′(x) · z be the scalar product of ℝ^{2L}. The Taylor expansion of h′_{z} at x ∈ Δ^{n} reads

0 = h′_{z}(q_{l}) = h′_{z}(x) + d_{x}h′_{z}(q_{l} − x) + \frac{1}{2}d_{ξ}^{2}h′_{z}(q_{l} − x)  (4.1)

for a some ξ ∈ Δ^{n}. Assume that |h′_{z}| attains its maximum at x ∈ Δ^{n} and also that x lies in the interior of Δ^{n}. Otherwise the ensuing argument applied to some face Δ^{k} ⊂ ∂Δ^{n} will yield even better estimates. Since d_{x}h′_{z} = 0 we immediately get

|h′_{z}(x)| = \frac{1}{2}d_{ξ}^{2}h′_{z}(q_{l} − x) ≤ min_{l=0...n} \frac{1}{2}d_{ξ}^{2}h′_{z}(q_{l} − x) ≤ \frac{1}{2}C_{2}|z| min_{l=0...n} |q_{l} − x|^{2}
≤ \frac{1}{2}C_{2}|z|\frac{n}{n + 1} .

Therefore

|h′(x)| = max_{|z|=1} h′_{z}(x) ≤ \frac{1}{2}C_{2} \frac{n}{n + 1}

and

∥h∥_{∞} = ∥πh′∥_{∞} ≤ δ := ∥π∥_{∞} \frac{1}{2}C_{2} \frac{n}{n + 1} .  (4.2)

In order to estimate the differential we define

μ_{n} := max{∥v∥ | v ∈ ℝ^{n+1}, ∑_{l=0}^{n} vq_{l} = 0,

∃η ∈ ℝ, a ∈ Δ^{n} : |η + v(q_{l} − a)| ≤ \frac{1}{2}|q_{l} − a|^{2}, l = 0 … n, \} .

Eliminating η this becomes

μ_{n} = max{∥v∥ | v ∈ ℝ^{n+1}, ∑_{l=0}^{n} vq_{l} = 0,

∃a ∈ Δ^{n} : v(q_{l} − q_{m}) ≤ 1 + a^{2} − a(q_{l} + q_{m}), l = 0 … n, }

≤ max{∥v∥ | v ∈ ℝ^{n+1}, ∑_{l=0}^{n} vq_{l} = 0, ∃a ∈ Δ^{n} ∀l, m = 0 … n :
v(q_{l} − q_{m}) ≤ 1 + a^{2} − a(q_{l} + q_{m})}

A rough estimate for this is μ_{n} ≤ \sqrt{n + 1} using that 1 + a^{2} − a(q_{l} + q_{m}) ≤ 2 for a ∈ Δ^{n}.

Because of (4.1) the differential is

|dh′_{z}| ≤ μ_{n}C_{2}|z| ,

hence

∥dh′∥_{∞} ≤ C_{2}μ_{n} ,
∥dh∥_{∞} = ∥d(πh′)∥_{∞} ≤ ∥dπ∥_{∞}∥h′∥_{∞} + ∥π∥_{∞}∥dh′∥_{∞}
≤ ∥dπ∥_{∞} \frac{1}{2}C_{2} \frac{n}{n + 1} + ∥π∥_{∞}C_{2}μ_{n} .

Over Δ^{n} the operator P takes the form

P f(x) = A(x)d_{x}f + B(x)f(x)

with functions A: Δ^{n} → Hom(Hom(ℝ^{n}, ℂ^{N}), C^{N}) and B: Δ^{n} → Hom(ℂ^{N}, C^{N}). In terms of the operator norms we estimate

|P h| ≤ ∥A∥_{∞}∥dh∥_{∞} + ∥B∥_{∞}∥h∥_{∞} ≤ ϵ

with

ϵ := C_{2} (∥A∥_{∞} (∥dπ∥_{∞} \frac{1}{2} \frac{n}{n + 1} + ∥π∥_{∞}μ_{n}) + ∥B∥_{∞}∥π∥_{∞} \frac{1}{2} \frac{n}{n + 1}) . (4.3)

5 Estimates for the Derivatives of an Eigenvector

In this section we obtain explicit pointwise apriori estimates for the up to 2nd order derivatives of a unit eigenvector f for P with eigenvalue λ. In principle these estimates are computed from the Sobolev- and Garding- inequalities, but here we want to derive explicit expressions for the error estimates. To this end we adapt the proofs of these inequalities from [G], [S], [K], for instance, to our purpose.

Let {ϕ_{s}}_{s}, ϕ_{s}: M → [0, 1] be a partition of unity corresponding to the charts Φ_{s}: U_{s} ⟶^{≅} ℝ^{n} and choose functions ψ_{s} with compact support and such that ψ_{s} = 1 on the support of ϕ_{s}. We have f = ∑_{s} ϕ_{s}f. In the sequel we will work in one chart Φ = Φ_{s} and therefore drop the subscript s in the notation. We will also identify U_{s} via Φ_{s} with a subset of U ⊂ ℝ^{n} and E|U_{s} = U × ℂ^{N}. In particular we will not distinguish between f and f̃ as in the previous section.

In order to get a sufficiently smoothing parametrix we need to work with P^{d} instead of P for some d > 1 + n. In fact, if one can perform the inverse Fourier transform of the remainder R of the parametrix for P^{d} analytically, it suffices to take d > 2 + n/2. We will use a local parametrix Q i.e. a pseudodifferential operator of order −d such that

ϕg = QψP^{d}g + Rg  (5.1)

for any g with compact support in the chart Φ. The remainder R will also be pseudodifferential of order −d. We apply (5.1) to g = ψf which gives

ϕf = ϕψf = QψP^{d}ψf + Rψf = Qψψ̃P^{d}f + Rψf
= Qψψ̃λ^{d}f + Rψf  (5.2)

where ψ̃ is a 0-order operator defined by the relation P^{d}ψ = ψ̃P^{d}. From the expressions for Q and R as pseudodifferential operators we obtain pointwise estimates for the derivatives of f.

In the subsequent calculations integration will be over ℝ^{n} with (2π)^{−n/2} times the Lebesgue measure. The Fourier transform of a Schwartz class function g on ℝ^{n} is ĝ(ξ) := ∫ e^{−iξx}g(x)dx and the Fourier inversion formula becomes g(x) = ∫ e^{iξx}ĝ(ξ)dξ.

Below we derive pointwise estimates for ϕf, d(ϕf) and d^{2}(ϕf) in terms of the eigenvalue λ and the L^{2}-norm of ϕf. The corresponding quantities for f are readily computed from these. We use the multi-index notation x^{α} = x_{1}^{α_{1}} ⋯ x_{n}^{α_{n}}, d_{x}^{α} = \frac{∂^{α_{1}}}{∂x_{1}^{α_{1}}} ⋯ \frac{∂^{αn}}{∂x_{n}^{αn}}, |α| = α_{1} + ⋯ + α_{n} for α = (α_{1}, …, α_{n}) ∈ ℕ^{n}, x = (x_{1}, …, x_{n}) ∈ ℝ^{n}.

In the chart Φ we expand the operator ψ̃P^{d} as

ψ̃P^{d} = ∑_{β} p_{β}(x)D_{x}^{β} .

It has a smooth symbol p(x, η) = ∑_{β} p_{β}(x)η^{β} which has compact x-support. The parametrix Q we will work with is a pseudodifferential operator

Qg(x) := ∫ q(x, ξ)e^{ixξ}ĝ(ξ)dξ = ∫ q(x, ξ)e^{i(x−y)ξ}g(y)dydξ

of order −d. Its symbol q(x, ξ) is given by

q(x, ξ) = ∑_{j=0}^{−d} q_{−d−j}(x, ξ) .

where the q_{−d−j}(x, ξ) are homogeneous of degree −d − j in ξ and determined by solving the equations

∑_{|β|−|α|−d−j=l} \frac{(−i)^{|α|}}{α!} d_{ξ}^{α}q_{−d−j}(x, ξ)d_{x}^{α}p_{β}(x, ξ) = \{ \begin{matrix}ϕ(x)σ(ξ) & l = 0 \\ 0 & l = −1, …, −d \end{matrix}

(5.3)

Here we have fixed once and for all a bump function σ : ℝ^{n} → [0, 1] which is 0 near 0 and 1 outside a small neighbourhood of the origin.

In order to write down an explicit formula for R we consider the Taylor expansion of the symbol p(y, η) of ψ̃P^{d} at y at x ∈ supp(ϕ):

p(y, η) = ∑_{|α|,|β|≤d} \frac{1}{α!}d_{x}^{α}p_{β}(x)(y − x)^{α}η^{β} + ∑_{|α|=d+1,|β|≤d} r_{α,β}(x, y)(y − x)^{α}η^{β}

From this and the formula

d_{x}^{α}(gh) = ∑_{β+γ=α} \frac{α!}{β!γ!}d_{x}^{β}(g)d_{x}^{γ}(h)

we obtain R in the form

Rg(x) = ∫ r(x, ξ, y)e^{i(x−y)ξ}g(y)dydξ

with

r(x, ξ, y) = ∑_{|α|≤d} \frac{(−i)^{|α|}}{α!} d_{ξ}^{α}q(x, ξ)d_{x}^{α}p(x, ξ) − ϕ(x)

+ ∑_{|α|=d+1} \frac{1}{α!}d_{ξ}^{α}q(x, ξ) ∑_{γ+δ=β} \frac{β!}{γ!δ!}(−1)^{|γ|}ξ^{γ}d_{y}^{δ}r_{α,β}(x, y)

= r_{0}(x, ξ) + r̃(x, ξ, y)

where r_{0}(x, ξ) = ϕ(x)(σ(ξ) − 1) + truncation error in (5.3), and r̃(x, ξ, y) are explicitely known symbols of order −d − 1. With these explicit expressions for R and Q we compute the derivatives of ϕf from (5.2) and obtain

d_{x}^{α}(ϕf)(x) = d_{x}^{α}(Qψψ̃λ^{d} + Rψ)f
= ∫ d_{x}^{α} ((λ^{d}q(x, ξ)ψ(x)ψ̃(x) + r_{0}(x, ξ)) e^{ixξ}) f̂(ξ)dξ
+ ∫ d_{x}^{α} (r̃(x, ξ, y))e^{i(x−y)ξ}) ψ(y)f(y)dydξ

Finally the Cauchy-Schwartz inequality yields the estimate

\frac{|d_{x}^{α}(ϕf)(x)|_{∞}}{∥ψf∥_{L^{2}}} ≤ ∥ ∑_{β+γ=α} \frac{α!}{β!γ!}d_{x}^{β} (λ^{d}q(x, ξ)ψ(x)ψ̃(x) + r_{0}(x, ξ)) i^{|γ|}ξ^{γ}∥_{L^{2},ξ}

+ ∥∫ ∑_{β+γ=α} \frac{α!}{β!γ!}d_{x}^{β}r̃(x, ξ, y)i^{|γ|}ξ^{γ}dξ∥_{L^{2},y} .  (5.4)

References

[G] Gilkey, Peter B. Invariance theory, the heat equation, and the Atiyah-Singer index theorem Mathematics Lecture Series, 11. Publish or Perish, Inc., Wilmington, DE, 1984.

[K] Kumano-go, Hitoshi Pseudodifferential operators MIT Press, Cambridge, Mass.-London, 1981.

[L] Larson, Mats G. A posteriori and a priori error analysis for finite element approximations of self-adjoint elliptic eigenvalue problems SIAM J. Numer. Anal. 38 (2000), no. 2, 608–625 (electronic).

[N] Nakao, M. T. , Yamamoto, N. , Kimura, S. On the Best Constant in the Error Bound for the H_{0}^{1}-Projection into Piecewise Polynomial Spaces J. Approx. Theory 93, 491-500 (1998)

[S] Seeley, R. T. Complex powers of an elliptic operator in Singular Integrals, Proc. Sympos. Pure Math., Chicago, Ill., 1966, pp. 288–307, Amer. Math. Soc., Providence, R.I.