Spline functions, the biharmonic operator and approximate eigenvalues

Abstract

The biharmonic operator plays a central role in a wide array of physical models, such as elasticity theory and the streamfunction formulation of the Navier–Stokes equations. Its spectral theory has been extensively studied. In particular the one-dimensional case (over an interval) constitutes the basic model of a high order Sturm–Liouville problem. The need for corresponding numerical simulations has led to numerous works. The present paper relies on a discrete biharmonic calculus. The primary object of this calculus is a high-order compact discrete biharmonic operator (DBO). The DBO is constructed in terms of the discrete Hermitian derivative. However, the underlying reason for its accuracy remained unclear. This paper is a contribution in this direction, expounding the strong connection between cubic spline functions (on an interval) and the DBO. The first observation is that the (scaled) fourth-order distributional derivative of the cubic spline is identical to the action of the DBO on grid functions. It is shown that the kernel of the inverse of the discrete operator is (up to scaling) equal to the grid evaluation of the kernel of $\left[ \left( \frac{d}{dx}\right) ^4\right] ^{-1}$, and explicit expressions are presented for both kernels. As an important application, 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. The discrete eigenvalues are proved to converge (at an “optimal” $O(h^4)$ rate) to the continuous ones. Another 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.

Appendix A: The discrete biharmonic operator—generating polynomials

Consider again the discrete fourth-order equation

$$\begin{aligned} \delta ^4_x\mathfrak {u}=\mathfrak {f}, \end{aligned}$$

(A.1)

where $\mathfrak {u},\,\mathfrak {u}_x\in l^2_{h,0}.$

In this section we obtain a direct proof of Corollary 5.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 [4, Section 10.6, Eq. (10.137)] and was used as the main tool in proving Claim 2.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 interpretation”, as expressed in Corollary 5.2.

By (2.9), Eq. (A.1) can be rewritten as

$$\begin{aligned} \frac{12}{h^2}\left[ \frac{(\mathfrak {u}_{x})_{j+1}-(\mathfrak {u}_{x})_{j-1}}{2h}-\frac{\mathfrak {u}_{j+1}+\mathfrak {u}_{j-1}-2\mathfrak {u}_j}{h^2} \right] =\mathfrak {f}_j,\;\;\;1\le j\le N-1, \end{aligned}$$

(A.2)

where by (2.8),

$$\begin{aligned}&\displaystyle \frac{1}{6}(\mathfrak {u}_{x})_{j-1}+\frac{2}{3}(\mathfrak {u}_{x})_{j}+\frac{1}{6}(\mathfrak {u}_{x})_{j+1}=\frac{\mathfrak {u}_{j+1}-\mathfrak {u}_{j-1} }{2h},\;\;\;1\le j\le N-1, \end{aligned}$$

(A.3)

$$\begin{aligned}&\displaystyle \mathfrak {u}_0=\mathfrak {u}_N=(\mathfrak {u}_x)_0=(\mathfrak {u}_x)_N=0. \end{aligned}$$

(A.4)

The system  (A.2), (A.3A.4) 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:

$$\begin{aligned} p(z)=\sum _{j=1}^{N-1} \mathfrak {u}_jz^j,\;\;\;q(z)=\sum _{j=1}^{N-1} (\mathfrak {u}_x)_j z^j,\;\;\;\phi (z)=\sum _{j=1}^{N-1} \mathfrak {f}_j z^j. \end{aligned}$$

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

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

$$\begin{aligned}&\frac{1}{2h}\left[ (z^{-1}-z)q(z)-(\mathfrak {u}_x)_1+(\mathfrak {u}_x)_{N-1}z^N \right] \nonumber \\&\qquad \quad -\frac{1}{h^2}\left[ (z+z^{-1}-2)p(z)-\mathfrak {u}_1-z^N \mathfrak {u}_{N-1} \right] =\frac{h^2}{12}\phi (z). \end{aligned}$$

(A.5)

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

$$\begin{aligned}&\left( \frac{1}{6}z^{-1}+\frac{2}{3}+\frac{1}{6}z \right) q(z)-\frac{1}{6}(\mathfrak {u}_x)_1-\frac{1}{6}z^N(\mathfrak {u}_x)_{N-1}\nonumber \\&\quad =\frac{1}{2h}\left[ (z^{-1}-z)p(z)-\mathfrak {u}_1+z^N\mathfrak {u}_{N-1} \right] . \end{aligned}$$

(A.6)

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

$$\begin{aligned} \frac{1}{h^2}(z^2-2z+1)p(z)+\frac{1}{2h}(z^2-1)q(z)= & {} -\frac{h^2}{12}z\phi (z)+\frac{1}{h^2}\left[ \mathfrak {u}_1 z+z^{N+1}\mathfrak {u}_{N-1} \right] \nonumber \\&+\,\frac{1}{2h}\left[ -(\mathfrak {u}_x)_1z+(\mathfrak {u}_x)_{N-1}z^{N+1} \right] , \qquad \end{aligned}$$

(A.7)

$$\begin{aligned} \frac{1}{2h}(z^2-1)p(z) +\left( \frac{1}{6}z^2+\frac{2}{3}z+\frac{1}{6} \right) q(z)= & {} \frac{1}{2h}\left[ -\mathfrak {u}_{1}z+z^{N+1}\mathfrak {u}_{N-1}\right] \nonumber \\&+\frac{1}{6}(\mathfrak {u}_x)_1z+\frac{1}{6}z^{N+1}(\mathfrak {u}_x)_{N-1}. \end{aligned}$$

(A.8)

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

$$\begin{aligned} p(z)=\frac{12h^2}{(z-1)^4}\cdot r(z), \end{aligned}$$

(A.9)

where

$$\begin{aligned} r(z)= & {} \left( \frac{1}{6}z^2+\frac{2}{3}z+\frac{1}{6} \right) \cdot \frac{h^2}{12}\cdot z\phi (z)\nonumber \\&- \left( \frac{1}{6}z^2+\frac{2}{3}z+\frac{1}{6} \right) \left( \frac{1}{h^2}\left[ \mathfrak {u}_{1}z+z^{N+1}\mathfrak {u}_{N-1} \right] +\frac{1}{2h}\left[ -(\mathfrak {u}_x)_1z+(\mathfrak {u}_x)_{N-1}z^{N+1} \right] \right) \nonumber \\&+\frac{1}{2h}(z^2-1)\left[ \frac{1}{2h}\left[ -\mathfrak {u}_{1}z+z^{N+1}\mathfrak {u}_{N-1}\right] +\frac{1}{6}(\mathfrak {u}_x)_1z+\frac{1}{6}z^{N+1}(\mathfrak {u}_x)_{N-1} \right] \end{aligned}$$

(A.10)

It should be noted that the expression (A.10) 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

$$\begin{aligned} r(1)=r'(1)=r''(1)=r'''(1)=0. \end{aligned}$$

(A.11)

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}$.

$$\begin{aligned} r(1)= & {} \frac{h^2}{12}\cdot \phi (1)- \frac{1}{h^2}\left[ \mathfrak {u}_{1}+\mathfrak {u}_{N-1} \right] -\frac{1}{2h}\left[ -(\mathfrak {u}_x)_1+(\mathfrak {u}_x)_{N-1} \right] =0 \qquad \end{aligned}$$

(A.12)

$$\begin{aligned} r'(1)= & {} h^2\cdot \phi (1)+ \frac{h^2}{12}\cdot \phi '(1)-\frac{1}{h^2}\left[ \frac{5}{2}\mathfrak {u}_{1}+\left( N+\frac{3}{2}\right) \mathfrak {u}_{N-1} \right] \nonumber \\&-\frac{1}{2h}\left[ -\frac{7}{3}(\mathfrak {u}_x)_1+\left( N+\frac{5}{3}\right) (\mathfrak {u}_x)_{N-1} \right] =0 \end{aligned}$$

(A.13)

$$\begin{aligned} r''(1)= & {} \frac{7h^2}{36}\cdot \phi (1)+ \frac{h^2}{3}\cdot \phi '(1)+ \frac{h^2}{12}\cdot \phi ''(1)\nonumber \\&-\frac{1}{h^2}\left[ \frac{23}{6}\mathfrak {u}_1 +\left( \frac{5}{6}+2N+N^2 \right) \mathfrak {u}_{N-1} \right] \nonumber \\&-\frac{1}{2h}\left[ -\frac{10}{3}(\mathfrak {u}_x)_1 +\left( \frac{4}{3}+\frac{7}{3}N+N^2\right) (\mathfrak {u}_x)_{N-1}\right] =0 \end{aligned}$$

(A.14)

$$\begin{aligned} r'''(1)= & {} \frac{h^2}{12}\cdot \phi (1)+ \frac{7h^2}{12}\cdot \phi '(1)+ \frac{h^2}{2}\cdot \phi ''(1)+ \frac{h^2}{12}\cdot \phi '''(1)\nonumber \\&-\frac{1}{h^2}\left[ \frac{5}{2}\mathfrak {u}_1 + \left( -\frac{1}{2}+\frac{3}{2}N^2+N^3 \right) \mathfrak {u}_{N-1}\right] \nonumber \\&-\frac{1}{2h}\left[ -2(\mathfrak {u}_x)_1 +\left( N+2N^2+N^3 \right) (\mathfrak {u}_x)_{N-1}\right] =0. \end{aligned}$$

(A.15)

We set

$$\begin{aligned} m_0=\phi (1),\quad m_1=\phi '(1),\quad m_2=\phi ''(1),\quad m_3=\phi '''(1), \end{aligned}$$

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

$$\begin{aligned} p(z)= & {} \frac{1}{6N^4}\Big [\frac{z(z^2+4z+1)}{(z-1)^4}\cdot \phi (z)\nonumber \\&+\frac{z}{(z-1)^4}\Big ((-m_0+3m_1)-\frac{1}{N}(6m_1+6m_2)\nonumber \\&+\frac{1}{N^2}(6m_1+12m_2+3m_3)-\frac{1}{N^3}(2m_1+6m_2+2m_3) \Big )\nonumber \\&+\frac{4z^2}{(z-1)^4}\Big (-m_0+\frac{1}{N^2}(3m_1+3m_2)-\frac{1}{N^3}(2m_1+6m_2+2m_3)\Big )\nonumber \\&+\frac{z^3}{(z-1)^4}\Big (-(m_0+3m_1)+\frac{1}{N}(6m_1+6m_2)\nonumber \\&-\frac{1}{N^2}(6m_2+3m_3)-\frac{1}{N^3}(2m_1+6m_2+2m_3) \Big )\nonumber \\&+\frac{z^{N+1}}{(z-1)^4}\Big (-\frac{1}{N} (3m_1+3m_2) +\frac{1}{N^2}(6m_2+3m_3) +\frac{1}{N^3}(2m_1+6m_2 +2m_3) \Big )\nonumber \\&+\frac{z^{N+2}}{(z-1)^4}\Big (\frac{1}{N^2}(-12m_1-12m_2) +\frac{1}{N^3}(8m_1+24m_2+8m_3) \Big )\nonumber \\&+\frac{z^{N+3}}{(z-1)^4}\Big (\frac{1}{N}(3m_1+3m_2)-\frac{1}{N^2}(6m_1+12m_2+3m_3)\nonumber \\&+\frac{1}{N^3}(2m_1+6m_2+2m_3) \Big )\Big ]. \end{aligned}$$

(A.16)

Note that since p(z) is a polynomial of degree $N-1$, all terms $z^j$ ($j\ge N$) in fact cancel.

We explicitly compute the coefficient of the term $z^j$ ($1\le j\le N-1$) in p(z), which gives us $\mathfrak {u}_j$.

Using

$$\begin{aligned}&\frac{1}{(z-1)^4}= \frac{1}{6}\sum _{j=0}^\infty (j+1)(j+2)(j+3) z^{j},\quad \frac{z}{(z-1)^4}= \frac{1}{6}\sum _{j=0}^\infty j(j+1)(j+2) z^{j},\\&\frac{z^2}{(z-1)^4}= \frac{1}{6}\sum _{j=0}^\infty (j-1)j(j+1) z^{j},\quad \frac{z^3}{(z-1)^4}= \frac{1}{6}\sum _{j=0}^\infty (j-2)(j-1)j z^{j},\\&\quad \frac{z(z^2+4z+1)}{(z-1)^4} =\sum _{j=0}^\infty j^3 z^{j}, \end{aligned}$$

we have

$$\begin{aligned} \frac{z(z^2+4z+1)}{(z-1)^4}\phi (z)=\sum _{j=1}^{N-1} z^j \sum _{l=1}^{j-1} (j-l)^3 \mathfrak {f}_l+ ({\text{ terms } \text{ of } \text{ order } }\ge N), \end{aligned}$$

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

$$\begin{aligned} \mathfrak {u}_j= & {} \frac{j}{6N^4}\Bigg [ \frac{1}{j}\sum _{l=1}^{j-1} (j-l)^3 \mathfrak {f}_l +\frac{1}{6}(j+1)(j+2)\Big ((-m_0+3m_1)-\frac{1}{N}(6m_1+6m_2)\nonumber \\&+\frac{1}{N^2}(6m_1+12m_2+3m_3)-\frac{1}{N^3}(2m_1+6m_2+2m_3) \Big )\nonumber \\&+\frac{2}{3}(j-1)(j+1)\Big (-m_0+\frac{1}{N^2}(3m_1+3m_2)-\frac{1}{N^3}(2m_1+6m_2+2m_3)\Big )\nonumber \\&+\frac{1}{6}(j-2)(j-1)\Bigg (-(m_0+3m_1)+\frac{1}{N}(6m_1+6m_2) -\frac{1}{N^2}(6m_2+3m_3)\nonumber \\&-\frac{1}{N^3}(2m_1+6m_2+2m_3) \Bigg )\Bigg ]\nonumber \\= & {} \frac{j}{6N^4}\Bigg [ \frac{1}{j}\sum _{l=1}^{j-1} (j-l)^3 \mathfrak {f}_l+j\left( -j m_0+3 m_1\right) -\frac{1}{N}\cdot 6j(m_1+m_2)\nonumber \\&+\frac{3j}{N^2}\left( (j+1)m_1+(j+3)m_2+m_3 \right) -\frac{2j^2}{N^3}\left( m_1+3m_2+m_3 \right) \Bigg ]. \end{aligned}$$

We now note that

$$\begin{aligned} m_0= & {} \phi (1)=\sum _{k=1}^{N-1} \mathfrak {f}_k,\quad m_1=\phi '(1)=\sum _{k=1}^{N-1} k\mathfrak {f}_k,\\ m_2= & {} \phi ''(1)=\sum _{k=2}^{N-1} k(k-1)\mathfrak {f}_k,\quad m_3=\phi '''(1)=\sum _{k=2}^{N-1} k(k-1)(k-2)\mathfrak {f}_k, \end{aligned}$$

so that

$$\begin{aligned}&m_1+m_2= \sum _{k=1}^{N-1} [k+k(k-1)]\mathfrak {f}_k=\sum _{k=1}^{N-1} k^2\mathfrak {f}_k\\&m_1+3m_2+m_3=\sum _{k=1}^{N-1} k^3\mathfrak {f}_k, \end{aligned}$$

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

$$\begin{aligned} \mathfrak {u}_j= & {} \frac{1}{6N}\Big [ \sum _{k=1}^{j-1} \left( \frac{j}{N}{-}\frac{k}{N}\right) ^3 \mathfrak {f}_k{-} \left( \frac{j}{N}\right) ^2\cdot \sum _{k=1}^{N-1} \left( 1-\frac{k}{N} \right) \left( 2\frac{k}{N}\cdot \frac{j}{N}+ \frac{j}{N}-3\frac{k}{N} \right) \mathfrak {f}_k\Big ]\\= & {} \frac{1}{6N}\Big [ \sum _{k=1}^{j-1} \left( x_j-x_k\right) ^3 \mathfrak {f}_k+ x_j^2\cdot \sum _{k=1}^{N-1} \left( 1-x_k \right) ^2\left( 2(1-x_j)x_k+x_k-x_j \right) \mathfrak {f}_k\Big ]\\= & {} \frac{1}{6N}\Big [ (1-x_j)^2 \cdot \sum _{k=1}^{j-1} x_k^2\left( 2x_j(1-x_k)+x_j-x_k \right) \mathfrak {f}_k+\\&x_j^2\cdot \sum _{k=j}^{N-1} \left( 1-x_k \right) ^2\left( 2x_k(1-x_j)+x_k-x_j \right) \mathfrak {f}_k\Big ]. \end{aligned}$$

We have thus obtained

Proposition A.1

Defining the matrix elements

$$\begin{aligned} K^h_{j,k}={\left\{ \begin{array}{ll}\frac{1}{6N}\cdot (1-x_j)^2\cdot x_k^2 \left( 2x_j(1-x_k)+x_j-x_k \right) ,&{} 1\le k \le j\le N-1,\\ \frac{1}{6N}\cdot (1-x_k)^2\cdot x_j^2 \left( 2x_k(1-x_j)+x_k-x_j \right) ,&{}1\le j \le k\le N-1, \end{array}\right. } \end{aligned}$$

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

$$\begin{aligned} \mathfrak {u}_j= \sum _{k=1}^{N-1} K^h_{j,k} \mathfrak {f}_k. \end{aligned}$$

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