Spline functions, the discrete biharmonic operator and approximate eigenvalues

TL;DR

This paper explores the connection between cubic spline functions and the discrete biharmonic operator, providing theoretical insights and demonstrating optimal eigenvalue convergence and positivity properties in numerical approximations.

Contribution

It reveals the fundamental link between cubic splines and the discrete biharmonic operator, explaining its accuracy and establishing eigenvalue convergence and positivity results.

Findings

The discrete biharmonic operator's kernel matches the grid evaluation of the inverse fourth derivative.

Discrete eigenvalues converge at an optimal rate of O(h^4) to continuous eigenvalues.

A positivity result for the kernels shows order-preserving properties in both continuous and discrete cases.

Abstract

The biharmonic operator plays a central role in a wide array of physical models, notably in elasticity theory and the streamfunction formulation of the Navier-Stokes equations. The need for corresponding numerical simulations has led, in recent years, to the development of a discrete biharmonic calculus. The primary object of this calculus is a high-order compact discrete biharmonic operator (DBO). The numerical results have been remarkably accurate, and have been corroborated by some rigorous proofs. However, there remained the "mystery" of the "underlying reason" for this success. This paper is a contribution in this direction, expounding the strong connection between cubic spline functions (on an interval) and the DBO. It is shown in particular that the (scaled) fourth-order distributional derivative of the cubic spline is identical to the action of the DBO on grid functions. The DBO is constructed in terms of the discrete Hermitian derivative. A remarkable fact is that the kernel of the inverse of the discrete operator is (up to scaling) equal to the grid evaluation of the kernel of $\Big[\Big(\frac{d}{dx}\Big)^4\Big]^{-1} .$ Explicit expressions are presented for both kernels. The relation between the (infinite) set of eigenvalues of the fourth-order Sturm-Liouville problem and the finite set of eigenvalues of the discrete biharmonic operator is studied, and the discrete eigenvalues are proved to converge (at an "optimal" $O(h^4)$ rate) to the continuous ones. Another remarkable consequence is the validity of a comparison principle. It is well known that there is no maximum principle for the fourth-order equation. However, a positivity result is derived, both for the continuous and the discrete biharmonic equation, showing that in both cases the kernels are order preserving.

Tables

Table 1. Table 1 . First 4 4 4 eigenvalues (top row) and their numerical approximations using a grid of N = 10 − 60 𝑁 10 60 N=10-60 nodes.

	k=1	k=2	k=3	k=4
true eigenvalue	500.563902	3803.537080	14617.630131	39943.799006
N=10	500.521885	3800.689969	14567.617771	39493.816015
N=20	500.561614	3803.398598	14615.468848	39926.599754
N=30	500.563462	3803.511145	14617.236978	39940.722654
N=40	500.563764	3803.529031	14617.509451	39942.881883
N=50	500.563845	3803.533813	14617.581402	39943.430972
N=60	500.563874	3803.535512	14617.606815	39943.623511

Equations

$$x_j=jh,0\le j\le N,h=\frac{1}{N}.$$

$$\mathcal{A}=\left\{u\in H_0^2(\Omega ),u_j=f_j,j=0,1,...,N\right\}.$$

$$\Vert u{\Vert }_{H_0^2(\Omega )}^2={\int }_0^1|u^{\prime \prime }(x){|}^{2}dx,$$

$$I(u)=\underset{0}{\overset{1}{\int }}|u^{\prime \prime }(x){|}^{2}dx,u\in H_0^2(\Omega ).$$

$$I(s_f)<I(g),s_f\ne g\in \mathcal{A}.$$

$$I_{min}=\underset{u\in \mathcal{A}}{\inf }I(u).$$

$$I(u_n)\underset{n\to \infty }{\to }{I}_{min}.$$

$$u_n\underset{n\to \infty }{\overset{w}{\to }}v\in H_0^2(\Omega ).$$

$$I(v)\le \underset{n\to \infty }{\mathrm{liminf}}I(u_n),$$

$$\int_{0}^{1}|r^{\prime\prime}(x)|^{2}dx=\frac{1}{4}\Big{[}\int_{0}^{1}|v^{\prime\prime}(x)|^{2}dx+\int_{0}^{1}|w^{\prime\prime}(x)|^{2}dx+2\int_{0}^{1}v^{\prime\prime}(x)w^{\prime\prime}(x)dx\Big{]}\leq I_{min}.$$

$$\int_{0}^{1}v^{\prime\prime}(x)w^{\prime\prime}(x)dx=I_{min}=\Big{[}\int_{0}^{1}|v^{\prime\prime}(x)|^{2}dx\Big{]}^{\frac{1}{2}}\Big{[}\int_{0}^{1}|w^{\prime\prime}(x)|^{2}dx\Big{]}^{\frac{1}{2}}.$$

$$\underset{0}{\overset{1}{\int }}{s}_f^{\prime \prime }(x)v^{\prime \prime }(x)dx=0,$$

$$dimV\le 4N-3(N-1)-4=N-1.$$

$$s_f(x_j)=f_j,j=1,...,N-1,s_f(x_0)=s_f^{\prime }(x_0)=s_f(x_N)=s_f^{\prime }(x_N)=0.$$

$$x_j=jh,0\le j\le N,h=\frac{1}{N}.$$

$$f_j^{\ast }=f(x_j),0\le j\le N.$$

$$(\mathfrak{v},\mathfrak{w}{)}_h=h\sum _{j=0}^N{\mathfrak{v}}_j{\mathfrak{w}}_j,$$

$$|\mathfrak{v}{|}_h^2=h\sum _{j=0}^N{\mathfrak{v}}_j^2.$$

$$|\mathfrak{v}{|}_{\infty }=\underset{0\le j\le N}{\max }\left\{|{\mathfrak{v}}_j|\right\}.$$

$$l_{h,0}^2=\left\{\mathfrak{v},{\mathfrak{v}}_0={\mathfrak{v}}_N=0\right\}.$$

$$({\delta }_x\mathfrak{v}{)}_j=\frac{1}{2h}({\mathfrak{v}}_{j+1}-{\mathfrak{v}}_{j-1}),1\le j\le N-1,$$

$$({\delta }_x^2\mathfrak{v}{)}_j=\frac{1}{h^2}({\mathfrak{v}}_{j+1}-2{\mathfrak{v}}_j+{\mathfrak{v}}_{j-1}),1\le j\le N-1,$$

$$({\sigma }_x\mathfrak{v}{)}_j=\frac{1}{6}{\mathfrak{v}}_{j-1}+\frac{2}{3}{\mathfrak{v}}_j+\frac{1}{6}{\mathfrak{v}}_{j+1},1\le j\le N-1.$$

$${\sigma }_x=I+\frac{h^2}{6}{\delta }_x^2,$$

$$({\sigma }_x{\mathfrak{v}}_x{)}_j=({\delta }_x\mathfrak{v}{)}_j,1\le j\le N-1.$$

$$({\mathfrak{v}}_x{)}_0=({\mathfrak{v}}_x{)}_N=0.$$

$${\delta }_x^4\mathfrak{v}=\frac{12}{h^2}[{\delta }_x{\mathfrak{v}}_x-{\delta }_x^2\mathfrak{v}].$$

$$({\tilde{\delta }}_x^2\mathfrak{v}{)}_j=2({\delta }_x^2\mathfrak{v}{)}_j-({\delta }_x{\mathfrak{v}}_x{)}_j,1\le j\le N-1.$$

$$-{\widetilde{{\delta }_x}}^2=-{\delta }_x^2+\frac{h^2}{12}{\delta }_x^4.$$

$$q(x_j)={\mathfrak{v}}_j,q(x_{j\pm 1})={\mathfrak{v}}_{j\pm 1},q^{\prime }(x_{j\pm 1})=({\mathfrak{v}}_x{)}_{j\pm 1}.$$

$$({\widetilde{{\delta }_x}}^2\mathfrak{v}{)}_j=q^{\prime \prime }(x_j),({\delta }_x^4\mathfrak{v}{)}_j=q^{(4)}(x_j).$$

Full text

It is a pleasure to thank Jean-Pierre Croisille, Dalia Fishelov and Robert Krasny for very fruitful discussions. We thank M. Hansmann for calling our attention to the paper of Markus [15].

1. INTRODUCTION

The operator \Big{(}\frac{d}{dx}\Big{)}^{4} on the interval $[0,1]$ is certainly the simplest conceivable example of a fourth-order elliptic one-dimensional operator. As such, its spectral theory is very well understood  [6, Chapter 5] or  [10]. In classical terminology, its study is labeled as a “fourth-order Sturm-Liouville theory”. The numerical computation of the eigenvalues was carried out using a “Shannon-type” sampling method in  [5] , by “matrix methods” in  [17] and by finite element methods in  [2]. On the other hand, the analogous “discrete” treatment leaves much to be desired. By this we mean the construction of approximating finite difference operators so that their eigenvalues can be shown to converge to the spectrum of the differential operator. Clearly, a significant question is the possibility of obtaining a high rate of convergence in this case.

The goal of this paper is to fill this gap, by exploring the very interesting structural similarities between \Big{(}\frac{d}{dx}\Big{)}^{4} and a suitable discrete biharmonic operator (DBO) $\delta_{x}^{4}.$ The bridge between “continuous” and “ discrete” is achieved by using the classical cubic spline functions. It is well known that the convergence of finite-dimensional approximations to an infinite-dimensional, unbounded differential operator, does not entail the convergence of the respective spectra. One of the main results here is to provide an affirmative answer to this convergence property. The proof relies on the strong connection between the DBO and differential operations on spline functions.

A basic tool is the * discrete Hermitian derivative* on an interval, that gives a fourth-order accurate approximation to the derivative of a smooth function. It has been the cornerstone in the construction of a fourth-order discrete approximation to the one-dimensional biharmonic operator  [11] and its extension to the full fourth-order Sturm-Liouville problem  [4]. In the two-dimensional case, it has been used in the construction of a compact high-order finite difference scheme for the Navier-Stokes system in the pure streamfunction formulation  [3, Part II]. In this paper the full discrete elliptic theory of the DBO is exploited in the study of the discrete spectrum and its asymptotic behavior as the number of grid points increases to infinity.

The structure of the paper is as follows.

In Section  2 we recall the basic (classical) construction of cubic spline functions on an interval. For the convenience of the reader we provide the full (and standard) proofs of the essential properties of these functions, those that are used in the sequel.

In Section  3 we recall the definitions of the discrete finite difference operators, and in particular introduce the Hermitian derivative and the discrete biharmonic operator $\delta_{x}^{4}.$

In Section  4 we establish the equality of the Hermitian derivative and the derivative of the interpolating cubic spline. This is a fundamental fact connecting the two non-local fourth-order approximations of the derivative. We were unable to locate this remarkable fact in the literature, even though we are still convinced that such a classical fact should be well-known.

Then the connection between the discrete biharmonic operator $\delta_{x}^{4}$ and the interpolating cubic spline function is established. It is in fact the main theme of this paper. Recall that the cubic spline is a $C^{2}$ functions, with finite jumps of the third-order derivatives at grid points. The result here (Proposition  4.4) is that the sizes of these jumps are determined by the DBO acting on the grid values. We have not been able to locate such a result in the literature, even though it seems to be such a fundamental fact.

This connection enables us to prove, in Section  5 , positivity results for the continuous and discrete fourth-order operators (see Proposition  5.1 and Proposition  5.3). Recall that there is no maximum principle for the fourth-order operator. Once again, it seems to us that the positivity result should exist already in the literature, but we have not been able to locate it.

In Section  6 we first give the explicit form of the kernel (Green’s function) of the continuous operator. In the first instance, this kernel acts in $L^{2}(0,1).$ We then extend it to the negative Sobolev space $H^{-2}(0,1).$ This space includes all finite measures, and in particular all grid functions. Using the connection to cubic spline functions we establish the remarkable result that the discrete resolvent (namely, the kernel of $(\delta_{x}^{4})^{-1}$) is just the grid evaluation of the continuous kernel, up to scaling. Indeed, this can be viewed as an alternative, very natural, definition of the compact discrete biharmonic operator.

Finally, Section  7 is concerned with the eigenvalues of both the continuous and discrete operators. These eigenvalues (more precisely their inverses) are studied in terms of the “kernel tools” developed in the previous sections; the established connection between the discrete and continuous kernels implies that the discrete eigenvalues are actually obtained by a “Nyström method” [18].

The highlight of this section (and one of the main results of the entire paper) is the proof of the convergence of the discrete eigenvalues to the continuous ones, at an “optimal” fourth-order rate (Theorem  7.14). This result is obtained by combining two ingredients:

• •

A suitable adaptation (Lemma  7.12) of a more general abstract convergence theorem  [14, 15]. However, we have chosen to provide a self-contained, much simpler, proof, that builds on the analytic theory of finite-dimensional perturbations, as expounded in Kato’s classical book  [13].

• •

The dependence of the eigenvalues on the respective kernels, see Proposition  7.4.

In Appendix  A we use the approach of “generating polynomials” in order to give yet another explicit construction of the kernel of the discrete resolvent $(\delta_{x}^{4})^{-1}.$ In fact, this classical method enables us to establish a totally different point-of-view concerning the compact discrete operators used here, beginning with the Hermitian derivative. This approach has the advantage of being directly related to the definitions of the discrete operators, avoiding the “mediation” of spline functions. It is potentially applicable as a computational approach to similar (discrete) problems.

2. THE BASIC SETUP for CUBIC SPLINES

In what follows we consider the interval $\Omega=[0,1]$ with a uniform grid

[TABLE]

We fix values $f=\left\{f_{j}\right\}^{N}_{j=0}$ so that $f_{0}=f_{N}=0,$ and consider the family

[TABLE]

The space $H^{2}_{0}(\Omega)$ is the space of functions having first and second (distrbutional) derivatives in $L^{2}(\Omega)$ and vanishing, with their first-order derivatives, at the endpoints.

It is well known that the norm in $H^{2}_{0}(\Omega)$ can be defined by

[TABLE]

and we shall refer henceforth to this norm.

We consider the functional

[TABLE]

We are interested in a minimizer for this functional, restricted to $\mathcal{A}.$ Since the properties of this minimizer will be essential in the rest of this paper, we provide here the details of the proof of this classical fact of the calculus of variations. A purely algebraic proof can be found in  [1, Theorem 3.4.3] or  [7, Chapter IV, Cubic Spline Interpolation].

Claim 2.1**.**

The functional has a unique minimizer on $\mathcal{A}$ , which we designate as $s_{f},$

[TABLE]

Proof.

In fact, the functional $I(u)$ is strictly convex and $\mathcal{A}$ is convex, so the existence of a unique minimizer is guaranteed by general principles.

However, for the convenience of the reader, we provide a simple, straightforward proof. Since clearly $\mathcal{A}\neq\emptyset,$ and $I(u)\geq 0,$ we can define

[TABLE]

Let $\left\{u_{n}\right\}_{n=1}^{\infty}\subseteq\mathcal{A}$ be a sequence such that

[TABLE]

The boundedness of $\left\{u_{n}\right\}_{n=1}^{\infty}$ in the Hilbert space $H^{2}_{0}(\Omega)$ implies  [9, Appendix D.4] that, again passing to a subsequence without changing index, there is a weak limit,

[TABLE]

This weak limit satisfies

[TABLE]

so that it is indeed a minimizer.

To prove the uniqueness of such a minimizer, suppose that $w\in H^{2}_{0}(\Omega)$ is another minimizer. Let $r=\frac{v+w}{2}.$ Clearly $r\in\mathcal{A},$ and, by the Cauchy-Schwarz inequality,

[TABLE]

It follows that $r$ is also a minimizer and in particular

[TABLE]

As is well known, equality in the Cauchy-Schwarz inequality implies that $w^{\prime\prime}=\pm v^{\prime\prime}.$ The boundary conditions now yield $w=\pm v,$ and the constraints at the nodes finally force $w=v.$ ∎

Claim 2.2**.**

1. (1)

$s_{f}$* is a cubic polynomial in each interval $[x_{j},x_{j+1}],\quad j=0,1,...,N-1.$* 2. (2)

$s_{f}\in C^{2}_{0}(\Omega).$** 3. (3)

The previous two properties , supplemented by the constraints $s_{f}(x_{j})=f_{j},\quad j=1,...,N-1,$ and $s_{f}(x_{0})=s^{\prime}_{f}(x_{0})=s_{f}(x_{N})=s^{\prime}_{f}(x_{N})=0$ determine $s_{f}$ uniquely.

Proof.

The basic property of $s_{f}$ is that

[TABLE]

for all $v\in H^{2}_{0}(\Omega)$ vanishing at all nodes $x_{j},\quad 0\leq j\leq N.$

We obtain the first property by taking any $v\in C^{\infty}_{0}(x_{j},x_{j+1}),$ and integrating twice by parts. The second property follows by observing that the test functions can have arbitrary values $v^{\prime}(x_{j}),v^{\prime\prime}(x_{j}).$

Finally, consider the space $V\subseteq\left\{s\in C^{2}(\Omega)\cap H^{2}_{0}(\Omega)\right\},$ such that $s\in V$ is a cubic polynomial in each interval $[x_{j},x_{j+1}],\quad j=0,1,...,N-1.$ Since $s\in V$ has four parameters in each interval $[x_{j},x_{j+1}],$ subtracting the number of constraints at all interior nodes and the endpoints yields

[TABLE]

Since on the other hand we have for every $f=\left\{f_{j}\right\}^{N}_{j=0}$ so that $f_{0}=f_{N}=0,$ a corresponding $s_{f}\in V,$ it follows that $dim\,V=N-1,$ so that $s_{f}$ is the unique function in $V$ satisfying the constraints. ∎

Definition 2.3**.**

The function $s_{f}$ is called the (“type I”) cubic spline corresponding to the constraints

[TABLE]

Claim 2.4**.**

Consider the vectors $f=\left\{f_{j}\right\}^{N}_{j=0}$ such that $f_{0}=f_{N}=0.$ Then the map $f\hookrightarrow s_{f}\in H^{2}_{0}(\Omega)$ is one-to-one and linear.

Remark 2.5**.**

When we introduce discrete spaces of grid functions in the following section, we shall denote by $l^{2}_{h,0}$ the space of such vectors $f.$

Proof.

The fact that the map is one-to-one is obvious since $s_{f}$ determines $f.$ The linearity follows from the uniqueness part in Claim  2.2. Indeed, if $s_{f},s_{g}$ correspond to $f,g=\left\{f_{j},g_{j}\right\}^{N}_{j=0}$ so that $g_{0}=g_{N}=f_{0}=f_{N}=0,$ respectively, then $s_{f}+s_{g}\in V$ (the space introduced in the proof of Claim  2.2) and it satisfies the constraints corresponding to $f+g,$ hence $s_{f+g}=s_{f}+s_{g}.$

∎

3. ** SETUP and DEFINITION OF THE DISCRETE OPERATORS**

We equip the interval $\Omega=[0,1]$ with a uniform grid

[TABLE]

The approximation is carried out by grid functions $\mathfrak{v}$ defined on $\left\{x_{j},\,0\leq j\leq N\right\}.$ The space of these grid functions is denoted by $l^{2}_{h}.$ For their components we use either $\mathfrak{v}_{j}$ or $\mathfrak{v}(x_{j}).$

For every smooth function $f(x)$ we define its associated grid function

[TABLE]

The discrete $l^{2}_{h}$ scalar product is defined by

[TABLE]

and the corresponding norm is

[TABLE]

For linear operators $\mathcal{A}:l^{2}_{h}\to l^{2}_{h}$ we use $|\mathcal{A}|_{h}$ to denote the operator norm.

The discrete sup-norm is

[TABLE]

The discrete homogeneous space of grid functions is defined by

[TABLE]

Given $\mathfrak{v}\in l^{2}_{h,0}$ we introduce the basic (central) finite difference operators

[TABLE]

The cornerstone of our approach to finite difference operators is the introduction of the Hermitian derivative of $\mathfrak{v}\in l^{2}_{h,0},$ that will replace $\delta_{x}.$ It will serve not only in approximating (to fourth-order of accuracy) first-order derivatives, but also as a fundamental building block in the construction of finite difference approximations to higher-order derivatives.

First, we introduce the “Simpson operator”

[TABLE]

Note the operator relation (valid in $l^{2}_{h,0}$)

[TABLE]

so that $\sigma_{x}$ is an “approximation to the identity”.

The Hermitian derivative $\mathfrak{v}_{x}$ is now defined by

[TABLE]

Remark 3.1**.**

In the definition  (3.8), the values of $(\mathfrak{v}_{x})_{j},\,\,j=0,N,$ need to be provided , in order to make sense of the left-hand side (for $j=1,N-1$). If not otherwise specified, we shall henceforth assume that $\mathfrak{v}_{x}\in l^{2}_{h,0},$ namely

[TABLE]

In particular, the linear correspondence $l^{2}_{h,0}\ni\mathfrak{v}\to\mathfrak{v}_{x}\in l^{2}_{h,0}$ is well defined, but not onto, since $\delta_{x}$ has a non-trivial kernel.

The biharmonic discrete operator is given by (for $\mathfrak{v},\,\mathfrak{v}_{x}\in l^{2}_{h,0}),$

[TABLE]

We next introduce a fourth-order replacement to the operator $\delta^{2}_{x}$ (see  [3, Equation (10.50)(c)]),

[TABLE]

Note that, in accordance with Remark  3.1 the operator $\widetilde{\delta}_{x}^{2}$ is defined on grid functions $\mathfrak{v}\in l^{2}_{h,0},$ so that also $\mathfrak{v}_{x}\in l^{2}_{h,0}.$

The connection between the two difference operators for the second-order derivative is given by

[TABLE]

Remark 3.2**.**

Clearly the operators $\delta_{x},\,\delta_{x}^{2},\,\delta_{x}^{4}$ depend on $h,$ but for notational simplicity this dependence is not explicitly indicated.

The fact that the biharmonic discrete operator $\delta_{x}^{4}$ is positive (in particular symmetric) is proved in  [3, Lemmas 10.9, 10.10]. Therefore its inverse \Big{(}\delta^{4}_{x}\Big{)}^{-1} is also positive. In fact, it satisfies a strong coercivity property, that is also established in the aforementioned reference.

An interpretation to the finite-difference operators $~{}\widetilde{\delta_{x}}^{2}$ and $\delta_{x}^{4}$ is provided by the “polynomial approach” [3, Section 10.3], as follows.

Let $q(x)$ be a fourth-order polynomial such that

[TABLE]

Then

[TABLE]

The discrete biharmonic operator gives a very accurate approximation to the continuous one (“optimal 4-th order accuracy”) , as seen in the following claim  [3, Theorem 10.19] .

Claim 3.3**.**

Let $f(x)\in C^{4}(\Omega),\,\,\Omega=[0,1].$ Let $u(x)$ satisfy

[TABLE]

subject to homogeneous boundary conditions

[TABLE]

Then

[TABLE]

Remark 3.4**.**

The “$O(h^{4})$” here means that there exists a constant $C>0,$ depending only on $f,$ such that for all integers $N>1,$

[TABLE]

Observe that the grid functions in this estimate are defined on the grid of (the variable) mesh size $h.$

4. SPLINES , HERMITIAN DERIVATIVES and the DISCRETE BIHARMONIC OPERATOR

We use the notation of the previous section.

Let $\mathfrak{u}\in l^{2}_{h,0}$ be a grid function vanishing at the endpoints and let $s_{\mathfrak{u}}\in H^{2}_{0}(\Omega)$ be the corresponding spline function.

We use interchangeably the notation $\mathfrak{u}_{j}=\mathfrak{u}(x_{j}).$

Let $\mathfrak{u}_{x}$ be the Hermitian derivative of $\mathfrak{u},$ and we set at the endpoints

[TABLE]

Proposition 4.1**.**

For all interior nodes, $s_{\mathfrak{u}}^{\prime}(x_{j})=\mathfrak{u}_{x}(x_{j}),\quad 1\leq j\leq N-1.$

Proof.

To simplify notation we shift $x_{j}=0,$ so we need to show

[TABLE]

The quadratic part of $s_{\mathfrak{u}}$ is continuous, so the equality for this part follows from Simpson’s rule.

Thus we need only check for $s_{\mathfrak{u}}(x)=a^{\pm}x^{3}$ for $\pm x>0.$ But this can be verified directly. ∎

In addition to $\mathfrak{u}\in l^{2}_{h,0},$ let $\mathfrak{v}\in l^{2}_{h,0}$ be a grid function vanishing at the endpoints and let $s_{\mathfrak{v}}$ be the corresponding spline function. At the endpoints we impose again the boundary conditions  (4.1).

Claim 4.2**.**

The map $(\mathfrak{u},\mathfrak{v})\to\int\limits_{0}^{1}s_{\mathfrak{u}}^{\prime\prime}(x)s_{\mathfrak{v}}^{\prime\prime}(x)dx$ is a scalar product on $l^{2}_{h,0}.$

Proof.

In view of Claim  2.4 the map is bilinear. Furthermore , if $\int\limits_{0}^{1}|s_{\mathfrak{u}}^{\prime\prime}(x)|^{2}dx=0,$ then $s_{\mathfrak{u}}^{\prime\prime}\equiv 0$ and since $s_{\mathfrak{u}}\in H^{2}_{0}$ it follows that also $s_{\mathfrak{u}}\equiv 0,$ which implies $\mathfrak{u}=0.$ ∎

We denote by $\delta^{4}_{x}\mathfrak{u}$ the Stephenson fourth-order derivative of $\mathfrak{u}.$ It is interesting that the scalar product of the previous claim can be expressed in terms of this fourth-order derivative.

Proposition 4.3**.**

*Let $\mathfrak{u},\mathfrak{u}_{x},\mathfrak{v},\mathfrak{v}_{x}\in l^{2}_{h,0}.$

The discrete scalar product of $\delta^{4}_{x}\mathfrak{u}$ and $\mathfrak{v}$ satisfies*

[TABLE]

Proof.

Pick $j\in\left\{1,2,...,N-1\right\}$ and let $Q_{j}(x)$ be the fourth-order polynomial used in the construction of $(\delta^{4}_{x}\mathfrak{u})_{j},$ namely,

[TABLE]

Observe that the second line above follows from Proposition 4.1.

Consider the polynomial $Q_{j}-s_{\mathfrak{u}}$ in the interval $[x_{j},x_{j+1}].$ It is a fourth-order polynomial with double zeros at $x_{j},x_{j+1},$ so it must have the form

[TABLE]

and similarly

[TABLE]

However,

[TABLE]

by definition of the discrete biharmonic operator.

Let us now compute

[TABLE]

since the fourth-order derivative of $s_{\mathfrak{u}}$ vanishes identically in the interval.

By summation, and recalling that $s_{\mathfrak{u}}\in C^{2},$ we get

[TABLE]

From Equations  (4.4), (4.5) we get

[TABLE]

and inserting this in Equation  (4.7) yields

[TABLE]

∎

Proposition 4.4**.**

The jump of the third order derivatives of the cubic splines at the nodes is given by

[TABLE]

Proof.

Combine Equations  (4.8) and  (4.6). ∎

Remark 4.5**.**

In the literature (e.g.  [1, 7] one can find various expressions for the jump of the third order derivatives of the cubic spline. However Proposition  4.4 provides a new expression, that can be interpreted as a “fourth-order derivative” of the function at the node.

We can also interpret the second derivative of $s_{\mathfrak{u}}$ in terms of the finite difference operators. Recall that this derivative is continuous at the nodes.

Corollary 4.6**.**

The value of $s_{\mathfrak{u}}^{\prime\prime}(x_{j})$ is given by

[TABLE]

Proof.

From Equation  (4.4) we get

[TABLE]

By definition, $Q_{j}^{\prime\prime}(x_{j})=(\widetilde{\delta_{x}}^{2}\mathfrak{u})_{j}$ and from  (4.6) we have $A_{j}=\frac{1}{24}Q_{j}^{(4)}(x_{j})=\frac{1}{24}(\delta^{4}_{x}\mathfrak{u})_{j},$ hence

[TABLE]

∎

Remark 4.7**.**

Note that invoking the relation  (3.11) we obtain from  (4.11)

[TABLE]

4.1. EVALUATING the INTEGRAL $\int\limits_{0}^{1}|s_{\mathfrak{u}}^{\prime\prime}(x)|^{2}dx$

We first compute over a grid interval

[TABLE]

To simplify notation, we set $y=x-x_{j},$ so that $s(y)=s_{\mathfrak{u}}(x)$ is a cubic polynomial in $y\in[0,h].$ Writing

[TABLE]

we get readily

[TABLE]

and

[TABLE]

hence

[TABLE]

Since $s^{\prime}(y)$ is a quadratic polynomial, we have

[TABLE]

and

[TABLE]

Turning now back to the variable $x,$ and taking into account the equalities

[TABLE]

[TABLE]

we obtain

[TABLE]

and

[TABLE]

Remark 4.8**.**

Equation  (4.3) can then be used to define the discrete fourth-order derivative $\delta^{4}_{x}\mathfrak{u}$ when $\mathfrak{u},\mathfrak{u}_{x}\in l^{2}_{h,0}.$ From equation  (4.14) we obtain an explicit expression for $\delta^{4}_{x}\mathfrak{u}_{j},$ which is actually the Stephenson expression.

5. POSITIVITY

It is well known that there is (in general) no maximum principle for elliptic partial differential operators of order higher than two. For the biharmonic equation in multi-dimensional domains there exist versions of the principle that involve estimates of the gradient of the solution, see  [16] and references therein. Under Dirichlet boundary conditions (the only ones considered here) the preservation of positivity property means that $\Delta^{2}u\geq 0\Rightarrow u\geq 0.$ It is actually a property of the domain. The maximum principle implies preservation of positivity but of course not vice versa. In the multi-dimensional case (excluding the one-dimensional case) we refer to  [12] and references therein.

In our one-dimensional case we have the following proposition.

Proposition 5.1**.**

Let

[TABLE]

where $u\in H^{4}(\Omega)\cap H^{2}_{0}(\Omega).$ Then the following comparison principle holds.

If $f(x)\geq 0,\,\,x\in\Omega,$ then also $u(x)\geq 0,\,\,x\in\Omega.$

Proof.

Suppose to the contrary that for some $y\in(0,1)$ we have $u(y)<0.$ We can assume that $y$ is a minimum point for $u,$ so that

[TABLE]

Since $u^{\prime}$ vanishes at the endpoints, we infer that there are points

[TABLE]

such that

[TABLE]

Let

[TABLE]

Consider the function $v(x)=u^{\prime\prime}(x).$ It satisfies in the interval $[a,b]$ the inequality

[TABLE]

as well as $v(a)=v(b)=0$ and $v(y)\geq 0.$

The standard maximum principle now yields

[TABLE]

hence also $u^{\prime}(x)\equiv u^{\prime}(y)=0,\,x\in[a,b].$

If $a>0$ we get a contradiction since there is a point $\xi\in(0,a)$ with $u^{\prime\prime}(\xi)=0.$ Similarly if $b<1.$ We conclude that $u^{\prime}(x)\equiv 0,\,x\in[0,1],$ hence $u(x)\equiv u(y)<0,\,x\in[0,1].$ However this contradicts the boundary condition $u(0)=u(1)=0.$

∎

Remark 5.2**.**

In Section  6 below we derive an expression for the resolvent kernel  (6.3). Since it is easy to see that the kernel is nonnegative, we obtain another proof of Proposition  5.1.

5.1. POSITIVITY of the DISCRETE BIHARMONIC OPERATOR

We now show that the same positivity property holds also for the discrete biharmonic operator.

Proposition 5.3**.**

Let

[TABLE]

where $\mathfrak{u},\,\mathfrak{u}_{x}\in l^{2}_{h,0}.$ Then the following comparison principle holds.

If $\mathfrak{f}_{j}\geq 0,\,\,0\leq j\leq N,$ then also $\mathfrak{u}_{j}\geq 0,\,\,0\leq j\leq N.$

Proof.

Suppose to the contrary that $\mathfrak{u}_{j_{0}}<0$ for some index $1\leq j_{0}\leq N-1.$

Let $s_{\mathfrak{u}}\in C^{2}_{0}(\Omega)$ be the corresponding spline function. Since $s_{\mathfrak{u}}(x_{j_{0}})=\mathfrak{u}_{j_{0}}<0$ it follows that there exists a minimum point $y\in\Omega$ so that

[TABLE]

We have

[TABLE]

Since $s_{\mathfrak{u}}^{\prime}$ vanishes at the endpoints, we infer that there are points

[TABLE]

such that

[TABLE]

Let

[TABLE]

Let $w(x)=s_{\mathfrak{u}}^{\prime\prime}(x).$ The function $w$ is continuous and linear in grid intervals. In view of Proposition  4.4 we get, in the sense of distributions,

[TABLE]

where $\delta_{y}$ is the Dirac measure at $y.$

Since $w(a)=w(b)=0,$ the standard maximum principle yields

[TABLE]

hence

[TABLE]

and in particular $s_{\mathfrak{u}}^{\prime}(a)=s_{\mathfrak{u}}^{\prime}(b)=0.$

As in the proof of Proposition  5.1 we conclude that $a=0$ and $b=1,$ and therefore

[TABLE]

which is a contradiction to the boundary conditions.

∎

Corollary 5.4**.**

Let $\mathfrak{u}$ satisfy the conditions of Proposition  5.3. Let $s_{\mathfrak{u}}$ be the corresponding spline function. Then

[TABLE]

Proof.

The assumption that there exists a point $y\in(0,1)$ such that $s_{\mathfrak{u}}(y)<0$ leads to a contradiction; this follows from the proof of Proposition  5.3 . ∎

6. THE CONTINUOUS and DISCRETE RESOLVENT KERNEL

The operator $\mathcal{L}=d^{4}/dx^{4},$ with homogeneous boundary conditions ($\phi\in D(\mathcal{L})\Rightarrow\phi(0)=\phi^{\prime}(0)=\phi(1)=\phi^{\prime}(1)=0$) is positive definite (in particular self adjoint) with domain $D(\mathcal{L})=H^{4}([0,1])\cap H^{2}_{0}([0,1].$ We now consider the kernel of $\mathcal{L}^{-1},$ namely, Green’s function of the biharmonic problem

[TABLE]

where $u\in H^{4}(\Omega)\cap H^{2}_{0}(\Omega).$ A standard computation leads to the following

Claim 6.1**.**

The solution of  (6.1) is given by

[TABLE]

where

[TABLE]

Proof.

By the general theory, we verify that in the sense of distributions, for each fixed $y,$ as a function of $x,$

[TABLE]

where $\delta_{y}$ is the Dirac measure at $y.$ In addition, $K(x,y)$ is symmetric in $x,\,y$ and satisfies the homogeneous boundary conditions (as a function of $x$). ∎

6.1. EXTENDING the KERNEL to $H^{-2}(\Omega)$

The domain of \Big{(}\frac{d}{dx}\Big{)}^{4} (as a self-adjoint operator in $L^{2}(\Omega),$ subject to homogeneous boundary conditions) is $H^{2}_{0}(\Omega)\cap H^{4}(\Omega).$ When extended (in the sense of distributions) to $H^{2}_{0}(\Omega),$ it maps it to its dual $H^{-2}(\Omega)$  [9, Chapter 5]. On the other hand, the general theory (or a direct inspection of the expression  (6.3)) ensures that, for every fixed $x\in\Omega,$ we have $K(x,\cdot)\in H^{2}_{0}(\Omega).$ It follows that Equation  (6.2) can be extended to all $u\in H^{2}_{0}(\Omega)$ (or, alternatively, to all $f\in H^{-2}(\Omega)$) as

[TABLE]

where $<\cdot,\cdot>$ is the \Big{(}H^{2}_{0}(\Omega),\,H^{-2}(\Omega)\Big{)} coupling.

We now fix a mesh size $h=\frac{1}{N}$ and consider the grid functions $\mathfrak{u}\in l^{2}_{h,0}$ vanishing at the endpoints. As in Section  4 we let $s_{\mathfrak{u}}\in H^{2}_{0}(\Omega)$ be the corresponding spline function.

Let

[TABLE]

We note that $SP_{h}$ is a finite-dimensional subspace of $H^{2}_{0}(\Omega).$ However, it is not fully contained in $H^{4}(\Omega).$ Therefore, as observed above, we can extend the differential operator \Big{(}\frac{d}{dx}\Big{)}^{4} to the union \Big{[}H^{4}(\Omega)\cap H^{2}_{0}(\Omega)\Big{]}\cup SP_{h}.

As was shown in Proposition  4.4, the action of the operator on $SP_{h}$ is given by a combination of Dirac delta-functions at the nodes $x_{j},$ that can be written as an equality of grid functions

[TABLE]

The right-hand side in this equation is a finite measure, and we recall that, owing to the Sobolev embedding theorem, all finite measures are contained in $H^{-2}(\Omega).$

Thus, Equation  (6.4) takes here the form

[TABLE]

Corollary 6.2**.**

The discrete operator $(\delta^{4}_{x})^{-1}:l^{2}_{h,0}\to l^{2}_{h,0}$ is represented by a matrix $\left\{K^{h}_{i,j}\right\}_{1\leq i,j\leq N-1},$ explicitly given by

[TABLE]

where $K(x,y)$ is the resolvent kernel of \Big{(}\frac{d}{dx}\Big{)}^{4}, as in Equation  (6.3).

7. CONTINUOUS and DISCRETE EIGENVALUES

7.1. THE CONTINUOUS OPERATOR

We now consider the eigenvalues of the operator $\mathcal{L},$ introduced in Section  6.

The operator has a compact resolvent, and the kernel $K$ of $\mathcal{L}^{-1}$ is given in Claim  6.1. The spectrum of $\mathcal{L}$ consists of an increasing sequence of positive simple eigenvalues, which we designate as $\left\{0<\lambda_{1}<\lambda_{2}<...<\lambda_{k}<...\right\}.$

Since these eigenvalues play an important role in the sequel, we provide below the details of their evaluation, repeating the proof of  [6, Lemma 5.5.4].

Let $\phi\in H^{4}([0,1])\cap H^{2}_{0}([0,1])$ be a real eigenfunction

[TABLE]

Clearly, this function must be of the form

[TABLE]

where $\beta$ is real and $\beta^{4}=\lambda.$

The conditions $\phi(0)=\phi^{\prime}(0)=0$ clearly imply

[TABLE]

and $\phi(1)=0$ yields

[TABLE]

The remaining condition $\phi^{\prime}(1)=0$ yields

[TABLE]

Multiplying the two equations and invoking standard identities we get

[TABLE]

which is to be considered as the equation determining the discrete eigenvalues.

Changing $\beta\to-\beta$ we can keep $A,C$ unmodified but reverse the signs of $B,D.$ It therefore follows that for $-\beta<0$ (solution of  (7.3)) we get the same eigenfunction  (7.1) as for $\beta>0,$ and we can consider only positive $\beta.$

We therefore get the full set of eigenfunctions (for $\beta>0$ solving  (7.3)),

[TABLE]

where $A,B$ satisfy  (7.2).

In order to estimate the location of the eigenvalues it therefore suffices to consider the positive solutions of  (7.3). The following claim is easy to verify.

Claim 7.1**.**

Equation  (7.3) has a sequence of positive solutions as follows.

[TABLE]

The corresponding eigenvalues $\lambda_{0}=\beta_{0}^{4},\,\lambda^{(1)}_{k}=(\beta^{(1)}_{k})^{4},\,\lambda^{(2)}_{k}=(\beta^{(2)}_{k})^{4}$ of $\mathcal{L}$ are all simple.

We denote by

[TABLE]

the orthonormal set of the associated eigenfunctions.

7.2. THE DISCRETE OPERATOR

We simplify the notation above and denote by $\left\{0<\lambda_{1}<\lambda_{2}<\ldots<\lambda_{k}<\ldots\right\}$ the (infinite) sequence of eigenvalues of \mathcal{L}=\Big{(}\frac{d}{dx}\Big{)}^{4}.

Given $h=\frac{1}{N},$ let

[TABLE]

be the finite sequence of eigenvalues of $\delta^{4}_{x}.$

We denote by $\Gamma$ the sum

[TABLE]

and let

[TABLE]

Proposition 7.2**.**

There exists a constant $C>0,$ independent of $h,$ so that

[TABLE]

Proof.

We introduce the (infinite) set of reciprocals of the eigenvalues of $\mathcal{L},$ namely, the eigenvalues of the kernel $K(x,y)$  (6.3),

[TABLE]

while

[TABLE]

is the set of eigenvalues of $(\delta^{4}_{x})^{-1},$ corresponding to the discrete kernel $K^{h}$  (6.6).

By the standard trace formula, it follows that

[TABLE]

Since $K(x,x)=\frac{1}{3}x^{3}(1-x)^{3},$ the numerical values of $\Gamma$ and $C$ can easily be calculated, and it turns out that

[TABLE]

On the other hand

[TABLE]

so that  (7.6) is established (and even with an explicit constant). ∎

Remark 7.3**.**

Observe that $\Gamma_{h}$ is the discrete trapezoidal approximation to the integral for $\Gamma.$ By the standard estimate for the trapezoidal rule, we obtain

[TABLE]

with $C=\frac{1}{12}\max\limits_{0\leq x\leq 1}|(\frac{d}{dx})^{2}K(x,x)|=\frac{1}{96}.$

The fourth-order estimate  (7.6) is clearly a result of a closer inspection of the kernel $K.$

The “collective” estimate  (7.6) does not imply that an estimate of the form $\lambda_{i}^{-1}-\lambda_{h,i}^{-1}=O(h^{4})$ is valid, for any fixed value of the index $i.$ However, the next proposition provides a weaker statement in this direction. It will play a key role in the final, stronger Theorem  7.14 below.

Proposition 7.4**.**

For any fixed integer $i\geq 1$ there exist positive constants $C,\,h_{0}>0$ such that for any $0<h=\frac{1}{N}<h_{0}$ we have

[TABLE]

*where $\Lambda_{h}^{-1}$ is the set of reciprocals introduced in  (7.8). *

Proof.

Let $\phi_{i}(x)\in H^{2}_{0}(\Omega)$ be a normalized eigenfunction of \Big{(}\frac{d}{dx}\Big{)}^{4}, corresponding to $\lambda_{i}.$ Recall that $\phi_{i}\in C^{\infty}$ and \Big{(}\frac{d}{dx}\Big{)}^{-4}\phi_{i}=\lambda_{i}^{-1}\phi_{i}. Hence

[TABLE]

For simplicity, we denote by $\left\{x_{j}=jh,\,0\leq j\leq N\right\}$ the grid points , omitting the obvious dependence on $h.$

Let $\phi_{i}^{\ast}=\left\{\phi_{i}(x_{0}),\ldots,\phi_{i}(x_{k}),\ldots,\phi_{i}(x_{N})\right\}$ be the corresponding grid function.

In view of Claim  3.3 and Corollary  6.2 we have for all $0\leq k\leq N,$

[TABLE]

where here and below $C>0$ is a constant depending only on $\phi_{i}$ that changes from one estimate to the next. Using the notation  (6.6) this can be rewritten as

[TABLE]

that is

[TABLE]

On the other hand, the smoothness of the normalized $\phi_{i}$ yields

[TABLE]

The last two estimates imply the following estimate of the operator norm,

[TABLE]

for $h<h_{0}.$ By a standard result concerning resolvents of self-adjoint operators we conclude that

[TABLE]

which concludes the proof of the proposition. ∎

Remark 7.5**.**

Proposition  7.4 shows that in any neighborhood of $\lambda_{i}^{-1}$ there is a discrete eigenvalue $\lambda_{h,k}^{-1},$ provided $h>0$ is sufficiently small. Observe, however, that we cannot infer that, even the largest eigenvalue (of $\mathcal{L}^{-1}$) $\lambda_{1}^{-1}$ is the limit, as $h\to 0,$ of the largest discrete eigenvalue $\lambda_{h,1}^{-1}$ (of ($\delta_{x}^{4})^{-1}$). This is done in Theorem  7.7 below.

Remark 7.6**.**

In view of Corollary  6.2 the discrete eigenvalues in $\Lambda_{h}^{-1}$ are obtained by a “Nyström method” [18], namely, eigenvalues of the discretized kernel. The fact that for any fixed integer $i\geq 1$

[TABLE]

follows from  [18, Theorem 3]. Proposition  7.4 establishes an “optimal” $O(h^{4})$ rate to this convergence.

7.3. CONVERGENCE OF THE FIRST DISCRETE EIGENVALUE

For the first discrete eigenvalue $\lambda_{h,1}$ we can establish its convergence (as $h\downarrow 0$) to $\lambda_{1}$ as follows.

Theorem 7.7**.**

The sequence of the discrete first eigenvalues of $\delta_{x}^{4}$ converges to the first eigenvalue of the continuous operator $\mathcal{L}:$

[TABLE]

Proof.

We prove in fact that

[TABLE]

We first prove that

[TABLE]

Given $\varepsilon>0,$ it suffices to prove that there exists $h_{0}>0$ so that for any $0<h<h_{0},$

[TABLE]

Since $\lambda_{1}^{-1}$ is the greatest eigenvalue of the kernel $K,$ we have

[TABLE]

Remark that (see the proof of Proposition  7.4) the maximum is attained by $\phi_{1},$ the normalized eigenfunction corresponding to $\lambda_{1}.$ However we shall need an approximating compactly supported function.

Now let $u^{\varepsilon}\in C^{\infty}_{0}(0,1)$ be a normalized function , $\|u^{\varepsilon}\|_{L^{2}(0,1)}=1$ and such that

[TABLE]

Take $h_{0}>0$ sufficiently small, so that $u^{\varepsilon}$ vanishes in a neighborhood of the “edge” intervals $[0,h_{0}]\cup[1-h_{0},1].$

Let $h=\frac{1}{N}<h_{0}.$

For simplicity, we denote by $\left\{x_{j}=jh,\,0\leq j\leq N\right\}$ the grid points , omitting the obvious dependence on $h.$

Define a nonnegative step function

[TABLE]

Clearly $\|U^{\varepsilon}\|_{L^{2}(0,1)}=1.$

The continuity of $K(x,y)$ implies that (decreasing $h_{0}$ if necessary)

[TABLE]

Let $\mathfrak{u}^{\varepsilon}=(\mathfrak{u}^{\varepsilon}_{1},\ldots,\mathfrak{u}^{\varepsilon}_{N-1})\in l^{2}_{h,0}$ be the grid function defined by

[TABLE]

so that $|\mathfrak{u}^{\varepsilon}|_{h}=1.$

Employing the notation  (6.6), the inequality  (7.22) can be rewritten as

[TABLE]

From the maximum principle (see the notation introduced in Corollary  6.2),

[TABLE]

we infer that

[TABLE]

Combining  (7.21), (7.23) and  (7.25) we obtain

[TABLE]

The estimate  (7.18) is therefore established.

We now proceed to establish the reverse inequality

[TABLE]

Given $\varepsilon>0,$ it suffices to prove that there exists $h_{0}>0$ so that for any $0<h<h_{0},$

[TABLE]

Let $\mathfrak{u}^{h}\in l^{2}_{h,0},\,\,|\mathfrak{u}^{h}|_{h,0}=1,$ be an eigenvector corresponding to $\lambda_{h,1}$, so that

[TABLE]

Since the kernel $K^{h}$ is positive, we can assume that $\mathfrak{u}^{h}_{i}\geq 0,\,\,0\leq i\leq N.$

Let $u^{h}(x)$ be the nonnegative piecewise constant function defined by

[TABLE]

Clearly $\|u^{h}\|_{L^{2}(0,1)}=1$ so in view of  (7.20)

[TABLE]

We now replace the kernel $K(x,y)$ by the piecewise constant kernel

[TABLE]

By increasing $N$ if needed, the continuity of $K(x,y)$ implies that

[TABLE]

so that, by the Cauchy-Schwarz inequality,

[TABLE]

Observe that when changing $N$ we must also change $\mathfrak{u}^{h}$ (hence $u^{h}$), but since they are normalized this change does not affect the above estimate.

Combining  (7.31) and  (7.33) we obtain

[TABLE]

Now

[TABLE]

Thus  (7.28) is established and the proof is complete.

∎

Theorem  7.7 does not give any convergence rate for the difference $|\lambda_{1}-\lambda_{h,1}|.$ In what follows we consider this issue, using the basic variational tools.

We begin with a more general discussion.

Pick $\phi\in\left\{\phi_{1},...,\phi_{k}...\right\}$ a normalized eigenfunction of $\mathcal{L},$ with associated eigenvalue $\lambda\in\left\{0<\lambda_{1}<\lambda_{2}<\ldots<\lambda_{k}<\ldots\right\}.$

Applying the operator $\mathcal{L}$ to

[TABLE]

we get

[TABLE]

Since $\phi$ is normalized, we have

[TABLE]

and continuing in this fashion we see that all derivatives of $\phi$ are bounded by some power of $\lambda,$ and therefore in the estimates below we have a generic constant $C>0$ depending only on $\lambda.$

Let $\phi^{*}$ be the corresponding grid function, $\phi^{*}(x_{i})=\phi(x_{i}),\quad 0\leq i\leq N.$

Let $\mathfrak{v}\in l^{2}_{h,0}$ satisfy

[TABLE]

where also $\mathfrak{v}_{x}\in l^{2}_{h,0}.$

By the fourth order accuracy  (3.15) we know

[TABLE]

where $C$ is independent of $N=h^{-1},$ but depends of course on $\phi.$

It follows that

[TABLE]

Since $\phi$ is normalized, the truncation error for the trapezoid integration gives

[TABLE]

hence also

[TABLE]

Let $\bar{\mathfrak{v}}=\frac{\mathfrak{v}}{|\mathfrak{v}|_{h}},$ then it follows from  (7.38)

[TABLE]

Regarding the first eigenvalue, we can now show that $\lambda_{h,1}$ can exceed $\lambda_{1}$ by at most $O(h^{4}).$

Claim 7.8**.**

Let $\lambda_{1}$ be the first eigenvalue of $\mathcal{L}$ ( by  (7.5), $\lambda_{1}=\beta_{0}^{4}$). Then there exists a constant $C>0,$ depending on the eigenfunction $\phi_{1},$ but not on $h,$ such that

[TABLE]

Proof.

Consider  (7.41) with $\lambda=\lambda_{1}.$ By the variational minimum principle for the first eigenvalue we know that

[TABLE]

hence

[TABLE]

which proves the claim.

∎

Remark 7.9**.**

The exact first eigenvalue is $\lambda_{1}=500.5639017404.$ Numerical calculations actually show that $\lambda_{h,1}\leq\lambda_{1},$ and that $\lambda_{h,1}$ increases as $h$ decreases. This is shown in Figure  1 . We are still unable to prove this monotonicity.

Remark 7.10**.**

Observe that in Claim  7.8 we do not have a corresponding lower limit, namely, that $\lambda_{h,1}$ is above $\lambda_{1}-O(h^{4}).$ This is evident in the numerical results displayed in Figure  2. The proof of this fact is postponed to Theorem  7.14 below, where we show that the convergence of all discrete eigenvalues to the corresponding continuous ones is “optimal”, namely, at an $O(h^{4})$ rate.

7.4. CONVERGENCE OF THE DISCRETE EIGENVALUES $\lambda_{h,k},\,\,k>1$

[TABLE]

We now consider the convergence of all discrete eigenvalues to their continuous counterparts.

Numerical simulations indicate that, if we *fix an index $k,$ * then

[TABLE]

with $C>0$ depending on $k.$ This is demonstrated in Figure  3 (for $N=16$) and Figure  4 (for $N=64$). We thank Jean-Pierre Croisille for both figures. Thus, even a very coarse resolution produces excellent approximation of the eigenvalues.

The convergence result in Theorem  7.7, that dealt with the first eigenvalue, did not yield an “optimal” convergence rate, as noted in Remark  7.10.

Using a very different approach, we shall now extend the convergence to all eigenvalues, and, furthermore, obtain the optimal $O(h^{4})$ convergence rate.

Let $K_{h}(x,y)$ be the piecewise constant (positive definite) kernel introduced in  (7.32). We denote by $\mathcal{L}_{h}^{-1}$ the operator (on $L^{2}[0,1]$) whose kernel is $K_{h}.$ Clearly this operator is compact and positive definite. In fact, the following claim asserts that it has only finitely many positive eigenvalues (depending on $h,$ of course).

Claim 7.11**.**

The set of eigenvalues of $\mathcal{L}_{h}^{-1}$ is the finite set $\Lambda_{h}^{-1},$ defined in  (7.8).

Proof.

Let $u\in L^{2}[0,1]$ be an eigenfunction of $\mathcal{L}_{h}^{-1}.$ Thus, for some $\mu>0,$

[TABLE]

In particular, $u$ is piecewise constant

[TABLE]

hence (with $K^{h}$ as in Corollary  6.2)

[TABLE]

where the boundary values $u(x_{0})=u(x_{N})=0$ are included.

Thus $\mu$ is an eigenvalue of $(\delta_{x}^{4})^{-1},$ hence $\mu=\lambda_{h,k}^{-1}$ for some $1\leq k\leq N-1.$ ∎

We now proceed to establish the convergence of all discrete eigenvalues to the corresponding continuous ones. In fact, the following lemma is a special case of a theorem of Markus  [15, Corollary 5.3] concerning differences of eigenvalues of self-adjoint operators. A similar general theorem was proved (much later) by Kato [14]. However the generality of Kato’s theorem required an “extended enumeration” of the eigenvalues, adding values of boundary points of the essential spectra.

For the convenience of the reader we provide here a simple proof of the lemma, following the proof of (the finite-dimensional) Theorem 6.11 in  [13, Section II.6].

Lemma 7.12**.**

Let $h=\frac{1}{N},$ and let

[TABLE]

[TABLE]

be the sets introduced in  (7.7), (7.8), respectively.

Then there exists a constant $C>0,$ independent of $h,$ so that

[TABLE]

Proof.

Note that both $\mathcal{L}^{-1},\,\mathcal{L}_{h}^{-1},$ are Hilbert-Schmidt (hence compact) positive operators.

For $t\in[0,1]$ let

[TABLE]

which is also compact, positive self-adjoint operator. In particular, its spectrum (apart from [math]) consists of a descending sequence of positive eigenvalues

[TABLE]

In view of the discussion in  [13, Chapter VII.3.2] the functions $\mu_{k}^{-1}(t),\,1\leq k<\infty,$ are continuous, piecewise analytic functions of $t,$ and satisfy

[TABLE]

and

[TABLE]

In addition, there exists (for every fixed $t\in[0,1]$) a corresponding set of orthonormal functions (in $L^{2}(0,1)$)

[TABLE]

Pick an index $k\geq 1.$ The eigenvalue $\mu_{k}^{-1}(t)$ is continuous (in $t\in[0,1]$) and piecewise analytic, with finitely many singularities. The associated eigenfunction $\phi_{k}(x;t)$ is piecewise analytic in $t,$ with the same (finitely many) singularities. Thus, the equation

[TABLE]

can be differentiated with respect to $t$ (excluding the singularities) and we obtain

[TABLE]

Taking the scalar product with $\phi_{k}(x;t)$ we conclude that

[TABLE]

Integrating this equation and taking  (7.46) and  (7.47) into account we get

[TABLE]

The self-adjoint operator $\mathcal{A}=\mathcal{L}_{h}^{-1}-\mathcal{L}^{-1}$ is Hilbert-Schmidt, hence compact. Let $\left\{\gamma_{1},\gamma_{2},\ldots\right\}$ be the sequence of its non-zero eigenvalues (repeated according to multiplicity) with a corresponding orthonormal sequence of eigenfunctions $\left\{\chi_{1}(x),\chi_{2}(x),\ldots\right\}\subseteq L^{2}(0,1).$

Since $\phi_{k}(x;t)=\sum\limits_{j=1}^{\infty}(\phi_{k}(x;t),\chi_{j}(x))_{L^{2}(0,1)}\chi_{j}(x),$ Equation  (7.51) entails

[TABLE]

where $\sigma_{j,k}=\int_{0}^{1}(\phi_{k}(x;t),\chi_{j})_{L^{2}(0,1)}^{2}dt,\,\,1\leq j,k<\infty.$

By the orthonormality of the functions (in $x$)

[TABLE]

Let $\Phi$ be a real convex function on the real line, with $\Phi(0)=0.$ From Jensen’s inequality we get

[TABLE]

and summation over $k$ yields

[TABLE]

In particular , taking $\Phi(\xi)=\xi^{2}$ and noting  (7.52) we obtain

[TABLE]

The sum on the right-hand side is the square of the Hilbert-Schmidt norm of $\mathcal{A},$ which is $\int_{0}^{1}\int_{0}^{1}|K(x,y)-K_{h}(x,y)|^{2}dxdy,$ thus proving  (7.45).

∎

Remark 7.13**.**

Note that we obtained in particular

[TABLE]

This estimate is valid simultaneously for all $N-1$ eigenvalues. Fixing an index $k,$ we get in particular

[TABLE]

In view of Claim  7.1 we have $\lambda_{k}\approx k^{4}.$ Thus  (7.54) yields only an $O(h)$ convergence.

However it is seen in Table  1, that even with a small number of grid points, the first discrete eigenvalues approximate very well the continuous ones. We shall prove below that indeed the convergence is “optimal”.

We now proceed to prove the “optimal” estimate. Compare  (7.43) and Remark  7.10 in what concerns the first eigenvalue.

Theorem 7.14** (Optimal rate of convergence of discrete eigenvalues).**

Fix an integer $k\geq 1$ and consider the discrete eigenvalue $\lambda_{h,k}$ as a function of $h=\frac{1}{N},\,N=k+1,k+2,\ldots.$ Then there exists a constant $C>0,$ depending only on $k,$ such that

[TABLE]

Proof.

In view of  (7.54) we have

[TABLE]

The sequence $\Lambda^{-1}$ of the reciprocals of the exact eigenvalues (see  (7.7)) is monotone decreasing, so there exist $h_{0}>0$ and $\eta>0,$ such that

[TABLE]

Combined with Proposition  7.4 we infer that the only (inverse) eigenvalue that can be “close” to $\lambda_{k}^{-1}$ is $\lambda_{h,k}^{-1},$ and that

[TABLE]

thus concluding the proof of the theorem. ∎

Remark 7.15**.**

Observe that in the proof of Theorem  7.14 we relied on special properties of the kernel, via Proposition  7.4. Without using such information we obtain “sub-optimal” estimates. For example,  (7.45) implies

[TABLE]

which is not optimal, in view of Claim  7.1. Compare also to the estimate in  (7.6) which can be written as

[TABLE]

Remark 7.16**.**

The $O(h^{4})$ rate of convergence, as stated in Theorem  7.14, can be compared to the method of collocation approximation  [8]. In the case of the latter , achieving a similar rate of convergence requires the construction of an interpolating $C^{3}$ piecewise fifth-order polynomial function, and then using collocation at Gaussian points. The results here were obtained by using the discretized kernel (of the inverse operator). Owing to the observed connection between this kernel and the classical ($C^{2}$) cubic splines, the approximating eigenvalues are in fact those of the fourth-order (distributional) derivative of the interpolating cubic spline at the grid points (Proposition  4.4).

Appendix A ** THE DISCRETE BIHARMONIC OPERATOR: GENERATING POLYNOMIALS**

Consider again the discrete fourth-order equation

[TABLE]

where $\mathfrak{u},\,\mathfrak{u}_{x}\in l^{2}_{h,0}.$ In this section we obtain a direct proof of Corollary  6.2. In other words, we compute the matrix corresponding to the operator $(\delta^{4}_{x})^{-1},$ without recourse to the theory of cubic spline functions involved in the previous proof. In fact, an expression for this matrix has already been given in  [3, Section 10.6, Eq. (10.137)] and was used as the main tool in proving Claim  3.3. However, the expression there was a product of three matrices, based on the matrix representation of the Hermitian derivative. Thus, while allowing to obtain the aforementioned estimates, it did not yield an “explicit” form (such that can be used in a computer code in a straightforward way).

Remarkably, the methodology expounded here uses the discrete operators in a totally different way; it employs generating functions, and is a systematic approach that can also be applied to other problems. Although the computations involved require some work, it has the advantage of being a straightforward application of the definitions of the discrete operators. It should be mentioned that we first carried out the computation here, and it motivated our search for a parallel “functional interpreation”, as expressed in Corollary  6.2.

By (3.9), Equation (A.1) can be rewritten as

[TABLE]

where by (3.8),

[TABLE]

[TABLE]

The system  (A.2), (A.3) must be solved for $\left\{\mathfrak{u}_{j},(\mathfrak{u}_{x})_{j}\right\}_{j=1}^{N-1}.$

To do this, we introduce generating functions, which are polynomials of degree $N-1$ in the variable $z$:

[TABLE]

We know $\phi(z)$ and want to find $p(z),q(z)$.

Equations (A.2) can be encoded as the following equality of polynomials,

[TABLE]

Similarly, Equations (A.3) are equivalent to the following polynomial equality

[TABLE]

Multiplying (A.5),(A.6) by $z$, and rearranging, we have

[TABLE]

[TABLE]

We now solve the system of two linear equations (A.7),(A.8) for $p(z),q(z)$. It suffices to write the solution for $p(z)$, which is

[TABLE]

where

[TABLE]

It should be noted that the expression (A) contains the unknown quantities $\mathfrak{u}_{1},\mathfrak{u}_{N-1},(\mathfrak{u}_{x})_{1},(\mathfrak{u}_{x})_{N-1}$. Once we determine these quantities, (A.9) will give us the solution to the system (A.2),(A.3). To find these quantities, we exploit the following fact: since $p(z)$ is a polynomial, while the expression (A.9) contains $(z-1)^{4}$ in the denominator, it must be the case that $z=1$ is a root of $r(z)$ of multiplicity $4$, that is

[TABLE]

By differentiating $r(z)$ three times and then substituting $z=1$, we obtain $4$ equations for $\mathfrak{u}_{1},\mathfrak{u}_{N-1},(\mathfrak{u}_{x})_{1},(\mathfrak{u}_{x})_{N-1}$.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We set

[TABLE]

solve the linear system (A.12)-(A) for $\mathfrak{u}_{1},\mathfrak{u}_{N-1},(\mathfrak{u}_{x})_{1},(\mathfrak{u}_{x})_{N-1}$, and then substitute these values into (A),(A.9), to obtain the expression

[TABLE]

Note that since $p(z)$ is a polynomial of degree $N-1$, all terms $z^{j}$ ($j\geq N$) in fact cancel. We explicitly compute the coefficient of the term $z^{j}$ ($1\leq j\leq N-1$) in $p(z)$, which gives us $\mathfrak{u}_{j}$. Using

[TABLE]

[TABLE]

[TABLE]

we have

[TABLE]

so that the coefficient of $z^{j}$ ($1\leq j\leq N-1$) in $p(z)$ is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We now note that

[TABLE]

[TABLE]

so that

[TABLE]

[TABLE]

and using these results the expression for $\mathfrak{u}_{j}$ simplifies to

[TABLE]

[TABLE]

[TABLE]

We have thus obtained

Proposition A.1**.**

Defining the matrix elements

[TABLE]

we have that the solution of (A.1) is given by

[TABLE]

This expression is seen to be identical to  (6.5), so that Proposition  A.1 is a re-statement of Corollary  6.2.