<figure><img src="image_1.png" /></figure>

HAL
open science

Recent Advances in the Study of a Fourth-Order
Compact Scheme for the One-Dimensional Biharmonic
Equation

D Fishelov, M Ben-Artzi, Jean-Pierre Croisille

► To cite this version:

D Fishelov, M Ben-Artzi, Jean-Pierre Croisille. Recent Advances in the Study of a Fourth-Order Compact Scheme for the One-Dimensional Biharmonic Equation. Journal of Scientific Computing, 2012, 53 (1), pp.55-79. 10.1007/s10915-012-9611-x. hal-01283643

HAL Id: hal-01283643

https://hal.science/hal-01283643

Submitted on 5 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RECENT ADVANCES IN THE STUDY OF A FOURTH-ORDER COMPACT SCHEME FOR
THE ONE-DIMENSIONAL BIHARMONIC EQUATION

D. FISHELOV*, M. BEN-ARTZI*, AND J.-P. CROISILLE*

**Abstract.** It is well-known that non-periodic boundary conditions reduce considerably the overall accuracy of an approximating scheme. In previous papers the present authors have studied a fourth-order compact scheme for the one-dimensional biharmonic equation. It relies on Hermitian interpolation, using functional values and Hermitian derivatives on a three-point stencil. However, the fourth-order accuracy is reduced to a mere first-order near the boundary. In turn this leads to an "almost third-order" accuracy of the approximate solution. By a careful inspection of the matrix elements of the discrete operator, it is shown that the boundary does not affect the approximation, and a full ("optimal") fourth-order convergence is attained. A number of numerical examples corroborate this effect.

**Keywords:** discrete biharmonic operator, nonhomogeneous boundary conditions, fourth-order convergence, Hermitian interpolation, compact schemes.

# 1. INTRODUCTION

In this paper we discuss convergence in the sup-norm of a compact approximation to the one-dimensional biharmonic equation. We consider a *boundary value problem*, so that the values of the function and its derivative are given on the boundary. While the discrete scheme under consideration is fourth-order accurate, the truncation error deteriorates to first-order at the boundary, affecting presumably the rate of convergence of the global approximation. Indeed, in a previous work [10] we have obtained a “suboptimal” convergence rate of almost $O(h^3)$. We show here that, in fact, the convergence rate is optimal, namely $O(h^4)$.

The *Discrete Biharmonic Operator* (henceforth DBO) considered here is both compact and of high order accuracy. Such schemes have recently been gaining in popularity, as may be seen in the papers (a very partial list) [2, 5, 6, 7, 8, 9, 11, 14, 15, 16].

A convergence analysis was performed in [12, 13, 1] in cases where the accuracy of the scheme deteriorates near the boundary. In particular, in [12] and [13] a hyperbolic system of first order and a parabolic problem were analyzed in the case where extra boundary conditions were given in order to "close" the numerical scheme. It was shown in [12, 13] that if the accuracy of the extra boundary conditions is one less than that of the inner scheme, then the overall accuracy of the scheme is determined by the accuracy at inner points. In [1] it was proved for a parabolic equation that if the scheme is of order $O(h^\alpha)$ at inner points and of order $O(h^{\alpha-s})$ near the boundary, then if $s = 0, 1$ the accuracy of the scheme is $O(h^\alpha)$. However, if $s \ge 2$ then the overall accuracy the scheme is $O(h^{\alpha-s+3/2})$. In some sense our approach is an extension of the convergence analysis described in [12, 13, 1]. Here we treat a differential equation of order four. Since our scheme is of fourth order at interior points and of first order at near-boundary points, we have $\alpha = 4$ and $s = 3$. We show that the overall accuracy of the scheme does not deteriorate at all due to the lower-order approximation near the boundary.

Compact high-order schemes for the biharmonic equation can be traced back to Stephenson [16], who proposed such a scheme in two dimensions. The DBO studied here may be viewed as a one-dimensional analog of Stephenson's scheme.

In our approach, the DBO is obtained as a fourth-order derivative of an interpolating polynomial. This polynomial requires not only functional values at neighboring points, but also suitable approximate derivatives. It turns out that in order to maintain accuracy at high order, the approximate derivatives need to be evaluated as fourth-order accurate Hermitian approximations.

Here we investigate in detail the various mathematical features of the discrete approximation:

<footnote>Date: May 23, 2012.</footnote>

<footnote>* This paper is dedicated to Professor Saul Abarbanel on the occasion of his 80-th birthday.</footnote>

<footnote>* Partially supported by a French-Israeli scientific cooperation grant 3-1355. The first author was also supported by a research grant of Afeka-Tel-Aviv Academic College of Engineering.</footnote>
• The truncation error of the biharmonic operator.

• Optimal, fourth-order, convergence of the discrete solution to the continuous one.

In Section 2 we construct a compact fourth-order approximation to the the biharmonic operator. In particular, the approximation to the first-order derivative, called the Hermitian derivative, is described. The operators $\delta_x^2$, $\delta_x$ as well as the Hermitian derivative are studied in Section 3. Their matrix representations are given, as well as the fourth-order accuracy of the Hermitian derivative.

Section 4 is concerned with the study of the truncation error for the approximation of the fourth-order derivative, namely the discrete biharmonic operator. It is shown that it is of fourth order at interior points and of first order at near-boundary points.

Section 5 contains the optimal convergence of the discrete approximation to the exact solution of the one-dimensional biharmonic problem. The error (see Theorem 6) is shown to be of fourth-order in the discrete $l_h^2$ norm.

In Section 6 we show numerical results, which validate the fourth-order accuracy of the discrete solution of the one-dimensional biharmonic problem and some of its generalizations.

## 2. DERIVATION OF THREE-POINT COMPACT OPERATORS

We consider here the one-dimensional biharmonic equation on the interval $[a, b]$. For the simplicity of the presentation, we choose homogeneous boundary conditions. The one-dimensional biharmonic equation is

$$ (1) \qquad \begin{cases} u^{(4)}(x) = f(x), & a < x < b, \\ u(a) = 0, & u(b) = 0, \end{cases} \quad u'(a) = 0, u'(b) = 0. $$

We look for a high-order compact approximation to (1). We lay out a uniform grid $a = x_0 < x_1 < \dots < x_{N-1} < x_N = b$. Here $x_i = ih$ for $0 \le i \le N$ and $h = (b-a)/N$.

In what follows, we shall use the notion of grid functions. A grid function is a function defined on the discrete grid $\{x_i\}_{i=0}^N$. We denote grid functions with fraktur letters such as $u, v$. We have

$$ (2) \qquad u = (u(x_0), u(x_1), \dots, u(x_{N-1}), u(x_N)). $$

In addition, we denote by $u^* = (u(x_0), u(x_1), \dots, u(x_{N-1}), u(x_N))$ the grid function, which consists of the values of $u(x)$ at grid points.

We denote by $l_h^2$ the functional space of grid functions. This space is equipped with a scalar product and an associated norm

$$ (3) \qquad (u, v)_h = h \sum_{i=0}^{N-1} u(x_i)v(x_i), \quad |u|_h = (u, u)_h^{1/2}. $$

The subspace of grid functions, having zero boundary conditions at $x_0 = a$ and $x_N = b$, is denoted by $l_{h,0}^2$. For grid functions $u, v \in l_{h,0}^2$, we have

$$ (4) \qquad (u, v)_h = h \sum_{i=1}^{N-1} u(x_i)v(x_i). $$

We also define the sup norm for a grid function $u$

$$ (5) \qquad |u|_{\infty} = \max_{0 \le i \le N} |u(x_i)|. $$

We define the difference operators $\delta_x, \delta_x^2$ on grid functions by

$$ (6) \qquad \delta_x u_i = \frac{u_{i+1} - u_{i-1}}{2h}, \quad 1 \le i \le N-1, $$

$$ (7) \qquad \delta_x^2 u_i = \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2}, \quad 1 \le i \le N-1. $$

In these definitions the boundary values $u_0, u_N$ are assumed to be known.

Suppose that we are given data $u_{i-1}^*, u_i^*, u_{i+1}^*$ at the grid points $x_{i-1}, x_i, x_{i+1}$. In addition, we are given some approximations $u_{x,i-1}^*, u_{x,i+1}^*$ for $u'(x_{i-1}), u'(x_{i+1})$.
We seek a polynomial of degree 4

$$ (8) \qquad p(x) = u_i^* + a_1(x - x_i) + a_2(x - x_i)^2 + a_3(x - x_i)^3 + a_4(x - x_i)^4, $$

which interpolates the data $u_{i-1}^*$, $u_i^*$, $u_{i+1}^*$, $u_{x,i-1}^*$, $u_{x,i+1}^*$.

The coefficients $a_1, a_2, a_3, a_4$ of the polynomial are

$$ (9) \qquad \left\{ \begin{aligned} a_1 &= \frac{3}{4h}(u_{i+1}^* - u_{i-1}^*) - \left(\frac{1}{4}u_{x,i+1}^* + \frac{1}{4}u_{x,i-1}^*\right), \\ a_2 &= \frac{1}{h^2}(u_{i+1}^* + u_{i-1}^* - 2u_i^*) - \frac{1}{4h}(u_{x,i+1}^* - u_{x,i-1}^*) = \delta_x^2 u_i^* - \frac{1}{2}(\delta_x u_x^*)_i, \\ a_3 &= -\frac{1}{4h^3}(u_{i+1}^* - u_{i-1}^*) + \frac{1}{4h^2}(u_{x,i+1}^* + u_{x,i-1}^*), \\ a_4 &= -\frac{1}{2h^4}(u_{i+1}^* + u_{i-1}^* - 2u_i^*) + \frac{1}{4h^3}(u_{x,i+1}^* - u_{x,i-1}^*) = \frac{1}{2h^2}((\delta_x u_x^*)_i - \delta_x^2 u_i^*). \end{aligned} \right. $$

The coefficients above require the data $u_i^*$ and $u_{x,i}^*$. In the case where only the values of $u_i^*$ are given, then $\{u_{x,i}^*\}_{i=1}^{N-1}$ have to be evaluated in terms of $\{u_i^*\}_{i=0}^N$. Looking at the first equation in (9), we see that a natural candidate for $u_{x,i}^*$ is

$$ u_{x,i}^* = a_1. $$

This yields

$$ u_{x,i}^* = \frac{3}{4h}(u_{i+1}^* - u_{i-1}^*) - \left(\frac{1}{4}u_{x,i+1}^* + \frac{1}{4}u_{x,i-1}^*\right), $$

or equivalently

$$ (10) \qquad \frac{1}{6}u_{x,i}^* + \frac{2}{3}u_{x,i}^* + \frac{1}{6}u_{x,i+1}^* = \delta_x u_i^*. $$

This is by definition the Hermitian derivative. If we introduce the three-point operator $\sigma_x$ on grid functions by

$$ (11) \qquad \sigma_x v_i = \frac{2}{6}v_{i-1} + \frac{1}{3}v_i + \frac{1}{6}v_{i+1}, \quad 1 \le i \le N-1, $$

can rewrite (10) as

$$ (12) \qquad \sigma_x u_{x,i}^* = \delta_x u_i^*, \quad 1 \le i \le N-1. $$

Observe that a knowledge of

$$ (13) \qquad u_{x,0}^* = u_0', \quad u_{x,N}^* = u_N', $$

is needed in order to solve (12). In addition, we will invoke the following relation

$$ (14) \qquad \sigma_x = I + \frac{h^2}{6} \delta_x^2. $$

Natural approximations to $u''(x_i)$, $u'''(x_i)$, $u''''(x_i)$ are $2a_2, 6a_3, 24a_4$, respectively (see (9)). We use the notation $\tilde{\delta}_x^2 u^*$, $\delta_x^3 u^*$, $\delta_x^4 u^*$ for the following operators.

$$ (15) \qquad \left\{ \begin{aligned} \tilde{\delta}_x^2 u_i^* &= 2a_2 = 2\delta_x^2 u_i^* - (\delta_x u_x^*)_i, \\ \delta_x^3 u_i^* &= 6a_3 = \frac{1}{h^2}(u_{x,i+1}^* + u_{x,i-1}^* - 2u_{x,i}^*) = (\delta_x^2 u_x^*)_i, \\ \delta_x^4 u_i^* &= 24a_4 = \frac{12}{h^2}((\delta_x u_x^*)_i - \delta_x^2 u_i^*). \end{aligned} \right. $$

This suggests that $\delta_x^4 u_i^*$ is an approximation to the fourth-order derivative of $u$ at $x_i$, namely,

$$ (16) \qquad \delta_x^4 u_i^* = \frac{12}{h^2}((\delta_x u_x^*)_i - \delta_x^2 u_i^*). $$

This approximation, called the discrete biharmonic approximation, is the one-dimensional analog of the Stephen-son's scheme [16]. Note that, in the non-periodic setting, boundary values of $u_x$ should be given in order to compute $\delta_x^4$ at near boundary points $x_1, x_{N-1}$.
### 3. THE OPERATORS $\delta_x^2$, $\delta_x$ AND THE HERMITIAN DERIVATIVE

3.1. **Matrix representation of the Hermitian derivative.** Let us provide now some matrix representations of the operators appearing in the Hermitian gradient. Let $U \in \mathbb{R}^{N-1}$ be the vector corresponding to the grid function $u \in l_{h,0}^2$,

$$ (17) \qquad U = [u_1, \cdots, u_{N-1}]^T. $$

The vector corresponding to the grid function $\delta_x u$ is

$$ (18) \qquad \frac{1}{2h} KU, $$

where the matrix $K = (K_{i,m})_{1 \le i,m \le N-1}$ is the skew-symmetric matrix

$$ (19) \qquad K = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 \\ -1 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & -1 & 0 & 1 \\ 0 & \cdots & 0 & -1 & 0 \end{bmatrix}. $$

The matrix which corresponds to $\sigma_x$ is $P/6$, where $P$ is the positive definite $(N-1) \times (N-1)$ matrix

$$ (20) \qquad P = \begin{bmatrix} 4 & 1 & 0 & \cdots & 0 \\ 1 & 4 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 1 & 4 & 1 \\ 0 & \cdots & 0 & 1 & 4 \end{bmatrix}. $$

Thus, Equation (12) can be written as

$$ (21) \qquad \frac{1}{6} P U_x = \frac{1}{2h} K U, $$

where $U, U_x$ are the vectors corresponding to $u, u_x$, respectively.

In the sequel, we shall also need the matrix representation of $\delta_x^2$. The matrix $T$, which corresponds to $-h^2\delta_x^2$, is the $(N-1) \times (N-1)$ symmetric matrix

$$ (22) \qquad T = \begin{bmatrix} 2 & -1 & 0 & \cdots & 0 \\ -1 & 2 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & -1 & 2 & -1 \\ 0 & \cdots & 0 & -1 & 2 \end{bmatrix}. $$

The matrix $P$ is related to $T$ by

$$ (23) \qquad P = 6I - T. $$

Therefore, the matrix which corresponds to the operator $\sigma_x$ (restricted to $l_{h,0}^2$) is $P/6 = I - T/6$.

3.2. **The eigenvalues and the eigenvectors of $\delta_x^2$.** To simplify the notations, we assume from now on that $[a, b] = [0, 1]$, thus $Nh = 1$.

In order to prove the fourth-order accuracy of the scheme, we shall need the eigenvalues and eigenvectors corresponding to $\delta_x^2$, and thus to the matrix $T$. The eigenvalues of $T$ are

$$ (24) \qquad \lambda_j = 4 \sin^2(\frac{j\pi}{2N}), \quad j = 1, \cdots, N-1 $$

and the corresponding normalized eigenvectors are $Z^k = (Z_{1k}, \cdots, Z_{N-1,k})^T$ (with respect to the Euclidean norm in $\mathbb{R}^{N-1}$), where

$$ (25) \qquad Z_{jk} = \left(\frac{2}{N}\right)^{1/2} \sin \frac{k\pi}{N}, \quad 1 \le k, j \le N-1. $$

We denote the column vectors as $Z^k \in \mathbb{R}^{N-1}$ and the row vectors as $Z_j \in \mathbb{R}^{N-1}$.
The matrix $Z = (Z_{jk})_{1\le j,k\le N-1} \in \mathbb{M}_{N-1}(\mathbb{R})$ is an orthogonal positive-definite matrix. Thus,

$$ (26) \qquad Z^2 = ZZ^T = I_{N-1}. $$

It follows that the matrix $T$ satisfies

$$ (27) \qquad T = Z\Lambda Z^T, $$

where $\Lambda = \text{diag}(\lambda_1, \dots, \lambda_{N-1})$. The normalized vectors (with respect to $(|\cdot|_h)$), which diagonalize the operator $-\delta_x^2$, are the grid functions $z^k$, which are defined by

$$ (28) \qquad z_{jk} = Z_{jk}/h^{1/2}. $$

Equivalently, they may be written as (noting that $Nh = 1$)

$$ (29) \qquad z_{jk} = \sqrt{2} \sin \frac{k j \pi}{N}, \quad 1 \le k, j \le N-1. $$

We have

$$ (30) \qquad \begin{cases} z_{jk} = \sqrt{2} \sin(j \frac{k \pi}{N}), & j = 1, \dots, N-1, \quad k = 1, \dots, N-1 \\ z_{0k} = 0, & z_{Nk} = 0, \\ -\delta_x^2 z^k = \tilde{\lambda}_k z^k, & \tilde{\lambda}_k = \frac{4}{h^2} \sin^2\left(\frac{k \pi}{2N}\right), & k = 1, \dots, N-1. \end{cases} $$

**3.3. The accuracy of the Hermitian derivative.** Now we state a lemma, proved in [3], which indicates the fourth-order accuracy of the Hermitian derivative.

**Lemma 1.** Suppose that $u(x)$ is a smooth function on $[a, b]$ and let $u = u^*$. Then, the Hermitian derivative $u_x$, as obtained from the values $u(x_i)$, $0 \le i \le N$ by

$$ (31) \qquad (\sigma_x u_x)_i = (\delta_x u^*)_i, \quad 1 \le i \le N-1 $$

and

$$ (32) \qquad (u_x)_0 = (u')^*(x_0), \quad (u_x)_N = (u')^*(x_N), $$

has a truncation error $u_x - (u')^*$ of order $O(h^4)$. More precisely,

$$ (33) \qquad |u_x - (u')^*|_{\infty} \le C h^4 \|u^{(5)}\|_{L^{\infty}}. $$

### 4. THE DBO AND ITS TRUNCATION ERROR

As mentioned in Section 2 the approximation $\delta_x^4 u_i^*$, suggested in (16), may serve as approximation to $u^{(4)}(x_i)$. We refer to $\delta_x^4$ as the discrete biharmonic operator (DBO). We define

**Definition 2 (Discrete biharmonic operator (DBO))**. Let $u \in l_h^2$ be a given grid function. The discrete biharmonic operator is defined by

$$ (34) \qquad \delta_x^4 u_i = \frac{12}{h^2} (\delta_x u_{x,i} - \delta_x^2 u_i), \quad 1 \le i \le N-1. $$

Here $u_x$ is the Hermitian derivative of $u$ satisfying (12) with given boundary values $u_{x,0}$ and $u_{x,N}$.

Using (16) and (10), the solution of (1) may be approximated by the scheme

$$ (35) \qquad \begin{cases} (a) \quad \delta_x^4 u_i = f(x_i) \quad 1 \le i \le N-1, \\ (b) \quad \frac{1}{6} u_{x,i-1} + \frac{2}{3} u_{x,i} + \frac{1}{6} u_{x,i+1} = \delta_x u_i, \quad 1 \le i \le N-1, \\ (c) \quad u_0 = 0, u_N = 0, u_{x,0} = 0, u_{x,N} = 0. \end{cases} $$

The scheme in (35) is the one-dimensional restriction of the scheme proposed by Stephenson in [16]. In the sequel, this scheme is referred to as the one-dimensional Stephenson Scheme to the biharmonic equation. Note that it approximates both $u$ and $u'$ at the grid points.

We first study in detail its truncation error. Let $u(x)$ be a smooth function on $[a, b]$, such that $u(a) = u(b) = 0$, $u'(a) = u'(b) = 0$. We denote by $u^*$ its related grid function.
We begin by considering the action of $\sigma_x \delta_x^4$ at the interior point $x_i$.

$$ (36) \qquad \sigma_x \delta_x^4 u_i^* = \frac{1}{6} \delta_x^4 u_{i-1}^* + \frac{2}{3} \delta_x^4 u_i^* + \frac{1}{6} \delta_x^4 u_{i+1}^*, \quad 2 \le i \le N-2, $$

where $\sigma_x$ is the Simpson operator defined in (11). The right-hand side can be expressed as

$$ (37) \qquad \frac{12}{h^2} \left( \left( \frac{1}{6} \delta_x u_{x,i-1}^* + \frac{2}{3} \delta_x u_{x,i}^* + \frac{1}{6} \delta_x u_{x,i+1}^* \right) - \left( \frac{1}{6} \delta_x^2 u_{i-1}^* + \frac{2}{3} \delta_x^2 u_i^* + \frac{1}{6} \delta_x^2 u_{i+1}^* \right) \right). $$

Using the definition of $u_x^*$, the first term in this expression is

$$ (38) \qquad \begin{aligned} & \frac{1}{6} \delta_x u_{x,i-1}^* + \frac{2}{3} \delta_x u_{x,i}^* + \frac{1}{6} \delta_x u_{x,i+1}^* = \sigma_x \delta_x u_{x,i}^* = \delta_x \sigma_x u_{x,i}^* = \delta_x \delta_x u_i^* \\ & = \frac{1}{4h^2} (u_{i+2}^* - 2u_i^* + u_{i-2}^*), \quad 2 \le i \le N-2. \end{aligned} $$

The second term in (37) may be written as

$$ (39) \qquad \begin{aligned} & \frac{1}{6} \delta_x^2 u_{i-1}^* + \frac{2}{3} \delta_x^2 u_i^* + \frac{1}{6} \delta_x^2 u_{i+1}^* \\ & = \frac{1}{6h^2} (u_{i-2}^* + 2u_{i-1}^* - 6u_i^* + 2u_{i+1}^* + u_{i+2}^*), \quad 2 \le i \le N-2. \end{aligned} $$

Therefore, inserting (38)-(39) in (36), we have

$$ (40) \qquad \sigma_x \delta_x^4 u_i^* = \frac{1}{h^4} (u_{i-2}^* - 4u_{i-1}^* + 6u_i^* - 4u_{i+1}^* + u_{i+2}^*) = \delta_x^2 \delta_x^2 u_i^*, \quad 2 \le i \le N-2. $$

Thus, in the absence of boundaries, there is a strong connection between $\delta_x^4$ and $(\delta_x^2)^2$. Explicit estimates for $\sigma_x \delta_x^4 u_i^*$ at near boundary points $x_1, x_{N-1}$ are given below (see (42)). It results from this representation that $\sigma_x \delta_x^4$ actually coincides with the operator $(\delta_x^2)^2$ at points $x_i$, $2 \le i \le N-2$. Only at near boundary points, $i=1, i=N-1$, we have a “numerical boundary layer” effect. Let us now investigate the accuracy of the DBO.

The following proposition deals with the truncation error of the DBO.

**Proposition 3.** Suppose that $u(x)$ is a smooth function on $[a, b]$. Assume, in addition, that $u(a) = u(b) = 0$, $u'(a) = u'(b) = 0$. Let $u_i^* = u(x_i)$, $(u^{(4)})^*(x_i) = u^{(4)}(x_i)$ be the grid functions corresponding, respectively, to $u, u^{(4)}$. Then the DBO $\delta_x^4$ satisfies the following accuracy properties:

<figure><img src="image_3.png" /></figure>

$$ (41) \qquad |\sigma_x \delta_x^4 u_i^* - \sigma_x (u^{(4)})^*(x_i)| \le Ch^4 \|u^{(8)}\|_{L^\infty}, \quad 2 \le i \le N-2. $$

* At near boundary points $i=1$ and $i=N-1$, the fourth order accuracy of (41) drops to first order,

$$ (42) \qquad |\sigma_x \delta_x^4 u_1^* - \sigma_x (u^{(4)})^*(x_1)| \le Ch \|u^{(5)}\|_{L^\infty}, $$

with a similar estimate for $i=N-1$.

* The error in the energy norm is given by

$$ (43) \qquad |\delta_x^4 u^* - (u^{(4)})^*|_h \le Ch^{3/2} (\|u^{(5)}\|_{L^\infty} + \|u^{(8)}\|_{L^\infty}). $$

In the above estimates $C$ is a generic constant, that does not depend on $u$.

*Proof.* According to (40), we have

$$ (44) \qquad \sigma_x \delta_x^4 u_i^* = (\delta_x^2)^2 u_i^*, \quad i=2, \dots, N-2. $$

We now expand $(\delta_x^2)^2$ in Taylor series. We have

$$ (45) \qquad \delta_x^2 u_i = u''(x_i) + \frac{h^2}{12} u^{(4)}(x_i) + \frac{h^4}{360} u^{(6)}(\xi_i), \quad 1 \le i \le N-1, $$

where $\xi_i \in (x_i - h, x_i + h)$.

We note that for any smooth function $v(x)$, $x \in [a, b]$, we have for $a+h < x < b-h$

$$ (46) \qquad \frac{v(x+h) - 2v(x) + v(x-h)}{h^2} = v''(\eta), $$

where $\eta \in (x-h, x+h)$.
Applying $\delta_x^2$ to Equation (45) at the interior points $x_i$, $2 \le i \le N-2$, we obtain

$$ (47) \qquad \begin{aligned} \sigma_x \delta_x^4 u_i^* &= (\delta_x^2)^2 u_i^* \\ &= u^{(4)}(x_i) + \frac{h^2}{6} u^{(6)}(x_i) + p_i, \quad |p_i| \le C_1 h^4 \|u^{(8)}\|_{L^\infty}, \quad 2 \le i \le N-2. \end{aligned} $$

On the other hand, $\sigma_x(u^{(4)})^*(x_i)$ may be expanded around $x_i$, $2 \le i \le N-2$, as follows.

$$ \begin{aligned} \sigma_x (u^{(4)})^*(x_i) &= \left(I + \frac{h^2}{6} \delta_x^2\right) (u^{(4)})^*(x_i) \\ &= u^{(4)}(x_i) + \frac{h^2}{6} u^{(6)}(x_i) + q_i, \quad |q_i| \le C_2 h^4 \|u^{(8)}\|_{L^\infty}, \quad 2 \le i \le N-2. \end{aligned} $$

Therefore, subtracting this equation from (47), we obtain the estimate (41).

Consider now the near boundary point $x_1$. We set $\delta_x^4 u_0^* = (u^{(4)})^*(x_0)$ and then, using the definition of $\sigma_x \delta_x^4$, we have

$$ (48) \qquad \begin{aligned} \sigma_x \delta_x^4 u_1^* - \sigma_x (u^{(4)})^*(x_1) &= \left(\frac{2}{3} \delta_x^4 u_1^* + \frac{1}{6} \delta_x^4 u_2^*\right) - \left(\frac{2}{3} (u^{(4)})(x_1) + \frac{1}{6} (u^{(4)})(x_2)\right) \\ &= \frac{2}{3} \left(\delta_x^4 u_1^* - (u^{(4)})(x_1)\right) + \frac{1}{6} \left(\delta_x^4 u_2^* - (u^{(4)})(x_2)\right). \end{aligned} $$

First, we consider the terms evaluated at $x_1$. Recall that

$$ (49) \qquad \delta_x^4 u_1^* = \frac{12}{h^2} ((\delta_x u_x^*)_1 - \delta_x^2 u_1^*), $$

where $u_x^*$ is the Hermitian derivative of $u^*$. Using the boundary values $u_0^* = u_{x,0}^* = 0$, we have, in view of (33),

$$ (50) \qquad (\delta_x u_x^*)_1 = \frac{u_{x,2}^*}{2h} = u''(x_1) + \frac{h^2}{6} u^{(4)}(x_1) + r_1, \quad |r_1| \le Ch^3 \|u^{(5)}\|_{L^\infty}, $$

and

$$ (51) \qquad \delta_x^2 u_1^* = u''(x_1) + \frac{h^2}{12} u^{(4)}(x_1) + r_2, \quad |r_2| \le Ch^3 \|u^{(5)}\|_{L^\infty}. $$

Inserting the estimates (50), (51) in Equation (49), we obtain

$$ (52) \qquad \delta_x^4 u_1^* = u^{(4)}(x_1) + r_3, \quad |r_3| \le Ch \|u^{(5)}\|_{L^\infty}. $$

Next, for $x_2$ we have

$$ (53) \qquad \delta_x^4 u_2^* = \frac{12}{h^2} ((\delta_x u_x^*)_2 - \delta_x^2 u_2^*). $$

Expanding on the term $(\delta_x u_x^*)_2$ and using again (33), we have

$$ (54) \qquad (\delta_x u_x^*)_2 = \frac{u_{x,3}^* - u_{x,1}^*}{2h} = u''(x_2) + \frac{h^2}{6} u^{(4)}(x_2) + s_1, \quad |s_1| \le Ch^3 \|u^{(5)}\|_{L^\infty}. $$

For the second term $\delta_x^2 u_2^*$ we have, as in (51),

$$ (55) \qquad \delta_x^2 u_2^* = u''(x_2) + \frac{h^2}{12} u^{(4)}(x_2) + s_2, \quad |s_2| \le Ch^3 \|u^{(5)}\|_{L^\infty}. $$

Inserting the estimates (54), (55) in Equation (53), we obtain, as in (52),

$$ (56) \qquad \delta_x^4 u_2^* = u^{(4)}(x_2) + s_3, \quad |s_3| \le Ch \|u^{(5)}\|_{L^\infty}. $$

Combining the estimates for $r_3$ and $s_3$ and inserting them in (48), we obtain

$$ (57) \qquad |\sigma_x \delta_x^4 u_1^* - \sigma_x (u^{(4)})_1^*| \le Ch \|u^{(5)}\|_{L^\infty}, $$

which proves (42).

(iii) Let $t_i = \delta_x^4 u_i^* - (u^{(4)})_i^*$ be the truncation error for the fourth-order derivative approximation. We have

$$ (58) \qquad \sigma_x t = v, $$

where $v \in l_{h,0}^2$ satisfies the estimates established in the previous parts of the lemma

$$ (59) \qquad |v_1|, |v_{N-1}| \le Ch \|u^{(5)}\|_{L^\infty}, \quad |v_i| \le Ch^4 \|u^{(8)}\|_{L^\infty}, \quad 2 \le i \le N-2. $$
The representative matrix of $\sigma_x$ restricted to $l_{h,0}^2$ is $P/6 = I - T/6$. The eigenvalues of $P/6$ are

$$ (60) \qquad 1 - \frac{2}{3}\mu_k = 1 - \frac{2}{3}\sin^2\left(\frac{k\pi}{2N}\right). $$

The matrix norm of its inverse is

$$ (61) \qquad \left| (P/6)^{-1} \right|_2 = \max_{k=1, \dots, N-1} \left| \frac{1}{1 - \frac{2}{3} \sin^2\left(\frac{k\pi}{2N}\right)} \right| \le 3. $$

From (58), (59) and (61), we obtain,

$$ (62) \qquad |\mathbf{t}|_h \le C|\mathbf{v}|_h. $$

Finally, since

$$ (63) \qquad |\mathbf{v}|_h^2 \le Ch(2h^2 + \sum_{i=2}^{N-2} h^8) (\|u^{(5)}\|_{L^\infty}^2 + \|u^{(8)}\|_{L^\infty}^2) \le Ch^3 (\|u^{(5)}\|_{L^\infty}^2 + \|u^{(8)}\|_{L^\infty}^2), $$

we get (43).

<figure><img src="image_4.png" /></figure>

## 5. OPTIMAL RATE OF CONVERGENCE OF THE ONE-DIMENSIONAL STEPHENSON SCHEME

In order to prove the fourth-order convergence of the scheme, we invoke the matrix representation for the discrete biharmonic operator.

### 5.1. Matrix representation of the DBO.
We have shown in [4] that the matrix form of the DBO (see Definition 2) is obtained from the matrix form of operators $\mathbf{u} \mapsto \mathbf{u}_x$ (see (21)), $\mathbf{u} \mapsto \delta_x \mathbf{u}$ (see (18)) and $\mathbf{u} \mapsto \delta_x^2 \mathbf{u}$ (see (22)). Let $U \in \mathbb{R}^{N-1}$ be the vector corresponding to the grid function $\mathbf{u} \in l_{h,0}^2$.

Therefore, the matrix representation of $\mathbf{u} \mapsto \delta_x^4 \mathbf{u}$ is

$$ (64) \qquad SU = \frac{12}{h^2} \left[ \frac{3}{2h^2} KP^{-1}K + \frac{1}{h^2}T \right] U = \frac{6}{h^4} \left[ 3KP^{-1}K + 2T \right] U. $$

The fact that we deal with a boundary value problem, rather than a periodic one, means that $PK - KP \neq 0$. However, the commutator is non-zero only at near-boundary points. Using the precise form of this commutator, we get the following proposition.

**Proposition 4.** (i) The operator $\sigma_x \delta_x^4$ has the matrix form

$$ (65) \qquad PS = \frac{6}{h^4}T^2 + \frac{6}{h^4}[e_1(e_1 + KP^{-1}e_1)^T + e_{N-1}(e_{N-1} - KP^{-1}e_{N-1})^T], $$

where

$$ (66) \qquad e_1 = (1, 0, \dots, 0)^T, \quad e_{N-1} = (0, \dots, 0, 1)^T. $$

(ii) The symmetric positive definite operator $\delta_x^4$ (see (64)) has the matrix form

$$ (67) \qquad S = \frac{6}{h^4}P^{-1}T^2 + \frac{36}{h^4}(V_1V_1^T + V_2V_2^T), $$

where the vectors $V_1, V_2$ are

$$ (68) \qquad \left\{ \begin{array}{l} V_1 = (\alpha - \beta)^{1/2} P^{-1} \left( \frac{\sqrt{2}}{2}e_1 - \frac{\sqrt{2}}{2}e_{N-1} \right) \\ V_2 = (\alpha + \beta)^{1/2} P^{-1} \left( \frac{\sqrt{2}}{2}e_1 + \frac{\sqrt{2}}{2}e_{N-1} \right). \end{array} \right. $$

The constants $\alpha, \beta$ are

$$ (69) \qquad \left\{ \begin{array}{l} \alpha = 2(2 - e_1^T P^{-1} e_1) \\ \beta = 2e_{N-1}^T P^{-1} e_1. \end{array} \right. $$
**Remark 5.** In view of the positivity of $P^{-1}$, we have $0 \le e_1^T P^{-1}e_1 \le 1/2$ and $|e_{N-1}^T P^{-1}e_1| \le 1/2$, so that $3 \le \alpha \le 4$ and $|\beta| \le 1$. Thus, $(\alpha \pm \beta)^{1/2}$ are well defined.

5.2. **Error estimate for the one-dimensional Stephenson scheme.** In [3] we carried out an error analysis based on the coercivity of $\delta_x^4$. The analysis presented there was based on an energy ($l^2$) method and led to a "sub-optimal" convergence rate of $h^{3/2}$. In [10] we have improved this result by showing that the convergence rate is almost three (the error is bounded by $Ch^3 \log(|h|)$). Here we prove the optimal (fourth-order) convergence of the scheme.

In order to obtain an optimal convergence rate, we use the matrix structure of $\delta_x^4$ given in (67). Let $u$ be the exact solution of (1) and let $u$ be its approximation by the Stephenson scheme (35). Let $u^*$ be the grid function corresponding to $u$. We consider the error between the approximated solution $u$ and the collocated exact solution $u^*$,

$$ \epsilon = u - u^*, $$

The grid function $u^*$ satisfies

$$ (70) \qquad \delta_x^4 u_i^* = f^*(x_i) + \mathbf{r}_i, \quad 1 \le i \le N-1, $$

where $\mathbf{r}$ is by definition the truncation error. We later refer to Proposition 3 for estimates on $\mathbf{r}$.

The error $\epsilon = u - u^*$ satisfies

$$ (71) \qquad \begin{aligned} \delta_x^4 \epsilon_i &= -\mathbf{r}_i, & 1 \le i \le N-1, \\ \epsilon_0 &= 0, \quad \epsilon_N = 0, \quad \epsilon_{x,0} = 0, \quad \epsilon_{x,N} = 0. \end{aligned} $$

We prove the following error estimate.

**Theorem 6.** Let $u$ be the exact solution of (1) and assume that $u$ has continuous derivatives up to order eight on $[a, b]$. Let $u$ be the approximation to $u$, given by the Stephenson scheme (35). Let $u^*$ be the grid function corresponding to $u$. The, the error $\epsilon = u - u^*$ satisfies

$$ (72) \qquad |\epsilon|_h \le Ch^4, $$

where $C$ depends only on $f$.

*Proof.* Let $U, U^* \in \mathbb{R}^{N-1}$ be the vectors corresponding to $u, u^*$, respectively, and let $F$ be the vector corresponding to $f^*$. We denote by $E = U - U^*$ and $R$ the vectors corresponding to $\epsilon = u - u^*$ and $\mathbf{r}$, respectively.

Using the matrix representation (67), we can write Equations (1) and (70) in the form

$$ (73) \qquad SU = F, $$

and

$$ (74) \qquad SU^* = F + R. $$

We therefore have

$$ (75) \qquad SE = -R. $$

In view of (67) we have that

$$ (76) \qquad PSP = \frac{6}{h^4} T^2 P + \frac{36}{h^4} J J^T, $$

where

$$ (77) \qquad J = \frac{\sqrt{2}}{2} [(\alpha - \beta)^{1/2}(e_1 - e_{N-1}), (\alpha + \beta)^{1/2}(e_1 + e_{N-1})]. $$

Inverting PSP and multiplying by $PR$, we have

$$ (78) \qquad -P^{-1}E = P^{-1}S^{-1}R = (PSP)^{-1}PR. $$

Our goal is to bound the elements of $P^{-1}E$ by $Ch^4$. Note that by Proposition 3 we have

$$ (79) \qquad \begin{aligned} &|(PR)_1|, |(PR)_{N-1}| \le Ch, \\ &|(PR)_j| \le Ch^4, \quad 2 \le j \le N-2. \end{aligned} $$

Thus, we need to estimate $(PSP)^{-1}PR$. We decompose PSP as follows

$$ (80) \qquad PSP = GH^{-1}, $$
where

$$
(81) \qquad G = I + 6JJ^T P^{-1}T^{-2}, \quad H = \frac{h^4}{6} P^{-1}T^{-2},
$$

so that

$$
(82) \qquad (PSP)^{-1} = HG^{-1}.
$$

Note that with $L = (6/h^4)H$, $Q = 6JJ^T$, we have

$$
(83) \qquad G = I + QL.
$$

We first estimate the elements of the matrix H.

**Estimate of the elements of H.**

In what follows we use $C$ as expressing various constants that do not depend on $h$. As in (27), we can diagonalize $H$ by

$$
H = Z \Lambda' Z^T,
$$

where the *j*-th column of the matrix *Z* is *Z*<sup>*j*</sup>, as defined in (25). Recall that *P* = 6*I* − *T* (see (23)), and that the
eigenvalues *λ*<sub>*j*</sub> of *T* are given by (24). Therefore, the eigenvalues *κ*<sub>*j*</sub>, 1 ≤ *j* ≤ *N* − 1 of *P* are given by

$$
(84) \qquad \kappa_j = 6 - \lambda_j = 6 - 4 \sin^2\left(\frac{j\pi}{2N}\right), \quad 1 \le j \le N-1.
$$

The diagonal matrix $\Lambda'$ contains the eigenvalues of $H$, which can be written as

$$
\theta_j = \frac{h^4}{6} \lambda_j^{-2} \kappa_j^{-1} = \frac{h^4}{96 \sin^4\left(\frac{j\pi}{2N}\right) (6 - 4 \sin^2\left(\frac{j\pi}{2N}\right))}, \quad j = 1, \dots, N-1.
$$

The element $H_{i,k}$ of the matrix $H$ is

$$
H_{i,k} = \sum_{j=1}^{N-1} Z_{i,j} \theta_j Z_{j,k}.
$$

$$
(85) \qquad H_{i,k} = \sum_{j=1}^{N-1} \frac{h^4}{96} \frac{2}{N} \sin\left(\frac{ij\pi}{N}\right) \frac{\sin\left(\frac{jk\pi}{N}\right)}{\sin^4\left(\frac{j\pi}{2N}\right)(6-4\sin^2\left(\frac{j\pi}{2N}\right))}.
$$

We can now estimate the order of magnitude of the elements of $H$ as functions of $h$. In fact, we shall inspect separately the first and last columns of $H$ and the rest ($k = 2, \dots, N-2$). The reason is that writing

$$
(86) \qquad (HG^{-1}PR)_i = \sum_{k=1}^{N-1} H_{i,k} (G^{-1}PR)_k,
$$

we shall see that $(G^{-1}PR)_1, (G^{-1}PR)_{N-1}$ can only be estimated by $Ch^2$ (see (112) below), so that the additional accuracy should come from $H_{i,1}$, $H_{i,N-1}$. Consider first the elements $(i,k)$ of $H$ for $k=1, N-1$. It suffices to consider $k=1$.

$$
(87) \qquad H_{i,1} = \sum_{j=1}^{N-1} \frac{h^4}{96} \frac{2}{N} \sin\left(\frac{ij\pi}{N}\right) \frac{1}{\sin^4\left(\frac{j\pi}{2N}\right)(6-4\sin^2\left(\frac{j\pi}{2N}\right))} \sin\left(\frac{j\pi}{N}\right).
$$

Recall the elementary inequalities

$$
(88) \qquad \sin x \ge \frac{2}{\pi}x, \quad 0 \le x \le \frac{\pi}{2},
$$

$$
(89) \qquad |\sin x| \le |x|, \quad 2 \le 6 - 4 \sin^2\left(\frac{j\pi}{2N}\right) \le 6.
$$

Noting that $h = 1/N$ and using the estimate $|\sin(\frac{ij\pi}{N})| \le 1$, we obtain

$$
(90) \qquad |H_{i,1}| = |H_{1,i}| \le C \sum_{j=1}^{N-1} h^5 \frac{1}{(jh)^4} (jh) \le Ch^2, \quad i=2, \dots, N-2.
$$
Similarly, we have

$$ (91) \qquad C_1 h^3 \le H_{1,1} \le C \sum_{j=1}^{N-1} h^5 \frac{1}{(jh)^4} (jh)^2 \le C_2 h^3. $$

This estimate holds equally for $H_{N-1,N-1}$. For the other corner elements of $H$ we have

$$ (92) \qquad |H_{1,N-1}| = |H_{N-1,1}| \le C_2 h^3. $$

For $i, k = 2, \dots, N-2$ we have

$$ (93) \qquad |H_{i,k}| \le C \sum_{j=1}^{N-1} h^5 \frac{1}{(jh)^4} \le Ch. $$

Therefore, the orders of magnitude of the elements of $H$ are bounded by

$$ (94) \qquad \begin{bmatrix} Ch^3 & Ch^2 & \dots & Ch^2 & Ch^3 \\ Ch^2 & Ch & \dots & Ch & Ch^2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ Ch^2 & Ch & \dots & Ch & Ch^2 \\ Ch^3 & Ch^2 & \dots & Ch^2 & Ch^3 \end{bmatrix}. $$

**Estimate of the elements of $G^{-1}$**.

We show that $G$ is invertible and we estimate its elements. First note that the elements of $L$ are the elements of $H$ multiplied by $6/h^4$.

The matrix $Q$ is $(N-1) \times (N-1)$, but it has only four non-zero components at the corner positions,

$$ (95) \qquad Q_{1,1} = Q_{N-1,N-1} = 6\alpha, \quad Q_{1,N-1} = Q_{N-1,1} = 6\beta. $$

Therefore, $QL$ has only two non-zero rows - the first and the last. The first row is given by

$$ (QL)_{1,j} = 6(\alpha L_{1,j} + \beta L_{N-1,j}), \quad j = 1, \dots, N-1 $$

and the last row is given by

$$ (QL)_{N-1,j} = 6(\beta L_{1,j} + \alpha L_{N-1,j}), \quad j = 1, \dots, N-1. $$

Thus,

$$ (96) \qquad \left\{ \begin{array}{l} G_{1,1} = 1 + 6(\alpha L_{1,1} + \beta L_{N-1,1}) =: a_1, \\ G_{1,N-1} = 6(\alpha L_{1,N-1} + \beta L_{N-1,N-1}) =: a_{N-1}, \\ G_{1,j} = 6(\alpha L_{1,j} + \beta L_{N-1,j}) =: b_j, \quad j = 2, \dots, N-2 \end{array} \right. $$

and

$$ (97) \qquad \left\{ \begin{array}{l} G_{N-1,1} = 6(\beta L_{1,1} + \alpha L_{N-1,1}) = G_{1,N-1} = a_{N-1}, \\ G_{N-1,N-1} = 1 + 6(\beta L_{1,N-1} + \alpha L_{N-1,N-1}) = G_{1,1} = a_1, \\ G_{N-1,j} = 6(\beta L_{1,j} + \alpha L_{N-1,j}) = b_{N-j}, \quad j = 2, \dots, N-2, \end{array} \right. $$

where the symmetries of $L$ have been used. In rows 2, 3, ..., $N-2$ the matrix $G$ has 1 on the diagonal and otherwise it is zero.

The orders of magnitude of $a_1, a_{N-1}$ and $b_j$ ($2 \le j \le N-2$) follow from those of the elements of $L$. Namely, $|a_1|, |a_{N-1}| \le C/h$ and $|b_j| \le C/h^2$ for $j = 2, \dots, N-2$. In what follows we shall need lower bounds for $a_1$ and $a_1^2 - a_{N-1}^2$. From their definitions above it is seen that we need an inspection of the terms $L_{1,1} = (6/h^4)H_{1,1}$, $L_{1,N-1} = (6/h^4)H_{1,N-1}$. Using the definitions of $L_{1,1}$, and $L_{1,N-1}$, we obtain

$$ (98) \qquad \begin{aligned} L_{1,1} &> |L_{1,N-1}| \\ L_{1,1} \mp L_{1,N-1} &= \frac{h}{4} \sum_{\substack{j=1 \\ \text{even or odd}}}^{N-1} \frac{\sin^2 \frac{j\pi}{N}}{\sin^4 \frac{j\pi}{2N} (6-4\sin^2 \frac{j\pi}{2N})} \\ &\ge Ch \sum_{\substack{j=1 \\ \text{even or odd}}}^{N-1} \frac{(jh)^2}{(jh)^4} = \frac{C}{h}. \end{aligned} $$

In the above, take "j = even" for "-" and "j = odd" for "+".
Using (96) and the bounds $3 \le \alpha \le 4$, $|\beta| \le 1$ (Remark 5), we get in view of (98) and (91)

$$ (99) \qquad |a_1| \ge 6(3L_{1,1} - |L_{N-1,1}|) - 1 \ge 12L_{1,1} - 1 \ge \frac{C}{h}. $$

Next, we treat the difference $a_1^2 - a_{N-1}^2$. Since we have the upper bound $|a_1^2 - a_{N-1}^2| \le C_1/h^2$, we again need only a lower bound. We write the difference $a_1^2 - a_{N-1}^2$ as

$$ (100) \qquad a_1^2 - a_{N-1}^2 = [1 + 6(\alpha + \beta)(L_{1,1} + L_{1,N-1})] \cdot [1 + 6(\alpha - \beta)(L_{1,1} - L_{1,N-1})] $$

(using the symmetries of $L$). In view of (98) and $\alpha \ge 3$, $|\beta| \le 1$, we obtain

$$ (101) \qquad |a_1^2 - a_{N-1}^2| \ge C_2/h^2. $$

To compute the inverse of $G$, we apply Gaussian elimination using the following method. We perform operations on rows of $G$ and apply the same operations to the identity matrix $I$. When $G$ is transformed to the identity matrix, $I$ is transformed to $G^{-1}$.

We first divide the first and the last row of $G$ by $a_1$ and annihilate the terms $j = 2, \dots, N-2$ of both rows by subtracting suitable multiplies of rows 2, ..., $N-2$. add the result to the first row, for $j = 2, 3, \dots, N-2$. The result is $G_1$, where

$$ (102) \qquad G_1 = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 & 0 & \frac{a_{N-1}}{a_1} \\ 0 & 1 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 1 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \dots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 0 & 1 & 0 \\ \frac{a_{N-1}}{a_1} & 0 & 0 & \dots & 0 & 0 & 1 \end{bmatrix}. $$

The same operations on the identity matrix yield the matrix

$$ (103) \qquad I_1 = \begin{bmatrix} \frac{1}{a_1} & \frac{-b_2}{a_1} & \frac{-b_3}{a_1} & \dots & \frac{-b_{N-3}}{a_1} & \frac{-b_{N-2}}{a_1} & 0 \\ 0 & 1 & -0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 1 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \dots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 0 & 1 & 0 \\ 0 & \frac{-b_{N-2}}{a_1} & \frac{-b_{N-3}}{a_1} & \dots & \frac{-b_3}{a_1} & \frac{-b_2}{a_1} & \frac{1}{a_1} \end{bmatrix}. $$

In order to eliminate the non-zero element of $G_1$ in position $(N-1, 1)$, we subtract a suitable multiple of the first row and add the result to the last row, thus getting the transformed matrix $G_2$

$$ (104) \qquad G_2 = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 & 0 & \frac{a_{N-1}}{a_1} \\ 0 & 1 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 1 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \dots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 0 & 1 & 0 \\ 0 & 0 & 0 & \dots & 0 & 0 & \frac{a_1^2 - a_{N-1}^2}{a_1^2} \end{bmatrix}. $$

The corresponding matrix $I_2$ (obtained similarly from $I_1$) is

$I_2 =$

$$ (105) \qquad \begin{bmatrix} \frac{1}{a_1} & \frac{-b_2}{a_1} & \frac{-b_3}{a_1} & \dots & \frac{-b_{N-3}}{a_1} & \frac{-b_{N-2}}{a_1} & 0 \\ 0 & 1 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 1 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \dots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 0 & 1 & 0 \\ \frac{-a_{N-1}}{a_1^2} & \frac{b_2a_{N-1}-a_1b_{N-2}}{a_1^2} & \frac{b_3a_{N-1}-a_1b_{N-3}}{a_1^2} & \dots & \frac{b_{N-3}a_{N-1}-a_1b_3}{a_1^2} & \frac{b_{N-2}a_{N-1}-a_1b_2}{a_1^2} & \frac{1}{a_1} \end{bmatrix}. $$
Now we divide the last row of $G_2$ and $I_2$ by $\frac{a_1^2-a_{N-1}^2}{a_1^2}$. We get

$$ (106) \qquad G_3 = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 & 0 & \frac{a_{N-1}}{a_1} \\ 0 & 1 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & 1 & 0 \\ 0 & 0 & 0 & \cdots & 0 & 0 & 1 \end{bmatrix} $$

and

$$ I_3 = $$

$$ (107) \qquad \begin{bmatrix} \frac{1}{a_1} & \frac{-b_2}{a_1} & \frac{-b_3}{a_1} & \cdots & \frac{-b_{N-3}}{a_1} & \frac{-b_{N-2}}{a_1} & 0 \\ 0 & 1 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & 1 & 0 \\ \frac{-a_{N-1}}{a_1^2-a_{N-1}^2} & \frac{b_2 a_{N-1}-a_1 b_{N-2}}{a_1^2-a_{N-1}^2} & \frac{b_3 a_{N-1}-a_1 b_{N-3}}{a_1^2-a_{N-1}^2} & \cdots & \frac{b_{N-3} a_{N-1}-a_1 b_3}{a_1^2-a_{N-1}^2} & \frac{b_{N-2} a_{N-1}-a_1 b_2}{a_1^2-a_{N-1}^2} & \frac{a_1}{a_1^2-a_{N-1}^2} \\ \end{bmatrix} $$

Finally, we eliminate the $(1, N-1)$ element in $G_3$ by subtracting a multiple of the last row. The corresponding operation on $I_3$ yields the inverse $G^{-1}$ as

$$ G^{-1} = $$

$$ (108) \qquad \begin{bmatrix} \frac{a_1}{a_1^2-a_{N-1}^2} & \frac{a_{N-1}b_{N-2}-a_1b_2}{a_1^2-a_{N-1}^2} & \frac{a_{N-1}b_{N-3}-a_1b_3}{a_1^2-a_{N-1}^2} & \cdots & \frac{a_{N-1}b_3-a_1b_{N-3}}{a_1^2-a_{N-1}^2} & \frac{a_{N-1}b_2-a_1b_{N-2}}{a_1^2-a_{N-1}^2} & \frac{-a_{N-1}}{a_1^2-a_{N-1}^2} \\ 0 & 1 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & 1 & 0 \\ \frac{-a_{N-1}}{a_1^2-a_{N-1}^2} & \frac{a_{N-1}b_2-a_1b_{N-2}}{a_1^2-a_{N-1}^2} & \frac{a_{N-1}b_3-a_1b_{N-3}}{a_1^2-a_{N-1}^2} & \cdots & \frac{a_{N-1}b_{N-3}-a_1b_3}{a_1^2-a_{N-1}^2} & \frac{a_{N-1}b_2-a_1b_{N-2}}{a_1^2-a_{N-1}^2} & \frac{a_1}{a_1^2-a_{N-1}^2} \\ \end{bmatrix} $$

We give accurate estimates for the non-trivial elements of $G^{-1}$, those in the first and last rows. Using (99), (101) and the corresponding upper bounds, one readily observes that

$$ (109) \qquad |(G^{-1})_{1,1}|, |(G^{-1})_{N-1,N-1}|, |(G^{-1})_{N-1,1}|, |(G^{-1})_{1,N-1}| \le Ch. $$

Similarly, and using also $|b_j| \le C/h^2$, we get

$$ (110) \qquad |(G^{-1})_{1,j}|, |(G^{-1})_{N-1,j}| \le \frac{C}{h}, \quad j=2, \dots, N-2. $$

Therefore, the elements of $G^{-1}$ are bounded by

$$ G^{-1} = $$

$$ (111) \qquad \begin{bmatrix} Ch & C/h & C/h & \cdots & C/h & C/h & Ch \\ 0 & 1 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & 1 & 0 \\ Ch & C/h & C/h & \cdots & C/h & C/h & Ch \end{bmatrix}. $$
We can now bound the elements of $G^{-1}PR$ using (109)-(110) and (79).

$$ (112) \qquad \begin{aligned} & |(G^{-1}PR)_1| \le \sum_{k=1}^{N-1} |(G^{-1})_{1,k}| \cdot |(PR)_k| \\ &= |(G^{-1})_{1,1}| \cdot |(PR)_1| + \sum_{k=2}^{N-2} |(G^{-1})_{1,k}| \cdot |(PR)_k| + \\ & \qquad + |(G^{-1})_{1,N-1}| \cdot |(PR)_{N-1}| \\ & \le C_1 h \cdot h + C_2 (N-3)(1/h) \cdot h^4 \le Ch^2. \end{aligned} $$

Similarly, we have that $|(G^{-1}PR)_{N-1}| \le Ch^2$.

For $i = 2, ..., N-2$

$$ (113) \qquad |(G^{-1}PR)_i| = |(PR)_i| \le Ch^4. $$

Finally we consider the product $HG^{-1}PR$ (see (78), (82))

$$ (114) \qquad -P^{-1}E = HG^{-1}PR. $$

Combining the estimates (90)-(93) with (112)-(113), we obtain

$$ (115) \qquad \begin{aligned} & |(HG^{-1}PR)_i| \le \sum_{k=1}^{N-1} |H_{i,k}| \cdot |(G^{-1}PR)_k| \\ &= |H_{i,1}| \cdot |(G^{-1}PR)_1| + \sum_{k=2}^{N-2} |H_{i,k}| \cdot |(G^{-1}PR)_k| + |H_{i,N-1}| \cdot |(G^{-1}PR)_{N-1}| \\ & \le C_1 h^2 h^2 + C_2 (N-3) h h^4 \le Ch^4. \end{aligned} $$

Therefore,

$$ (116) \qquad |(P^{-1}E)_i| = |(HG^{-1}PR)_i| \le Ch^4, \quad 1 \le i \le N-1. $$

**Conclusion of the proof of Theorem 6.** Using (116) we obtain that the Euclidean norm of the vector $E = U - U^*$ satisfies the estimate

$$ (117) \qquad |E| = |PP^{-1}E| \le C|P^{-1}E| \stackrel{(116)}{\le} C\sqrt{\sum_{i=1}^{N-1}(h^4)^2} = Ch^{-1/2}h^4. $$

Thus, in view of the definition of the $l_2$ norm

$$ (118) \qquad |e|_h \le Ch^4. $$

This proves the fourth order error estimate result.

## 6. NUMERICAL RESULTS

In order to assess the spatial fourth-order accuracy of the scheme, we performed several numerical tests. In the tables below we show $e_{max}$ - the error in the maximum norm, and $e_2$ - the error in the $l_2$ norm.

$$ \begin{aligned} e_{max} &= \max |u_{comp} - u_{exact}|, \\ e_2 &= \|u_{comp} - u_{exact}\|_{l^2} = |u_{comp} - u_{exact}|_h. \end{aligned} $$

Here, $u_{comp}$ and $u_{exact}$ are the computed and the exact solutions, respectively.

We illustrate the numerical properties of the scheme (35) as follows.

* The scheme (35) is observed to be fourth-order accurate in the maximum and the discrete $l^2$ norms, whenever homogeneous or nonhomogeneous boundary conditions are applied. This is shown in Case 1.

* In case of highly oscillatory solutions, the scheme (35) behaves remarkably well. Case 2 describes the convergence of the scheme for such a family of solutions. In this case too fourth-order accuracy is observed. In addition, the magnitude of the errors is very small even for coarse grids.

* Finally, in Case 3 we show that the scheme can be also used for nonlinear biharmonic equations, retaining the fourth-order accuracy.

### 6.1. Case 1: Polynomial solutions.
We consider two polynomial solutions. The first corresponds to homogeneous boundary conditions and the second to inhomogeneous conditions.
6.1.1. *Homogeneous boundary conditions.* Consider the polynomial solution of $u^{(4)} = f$

$$u(x) = x^4(x-1)^2 \text{ on } [a, b] = [0, 1].$$

It satisfies

$$(119) \qquad u(0) = u'(0) = u(1) = u'(1) = 0.$$

Thus, choosing

$$(120) \qquad f(x) = u^{(4)} = 360x^2 - 240x + 24,$$

then $u(x)$ is the unique solution of the biharmonic problem

$$(121) \qquad \left\{ \begin{array}{ll} u^{(4)} = f, & 0 < x < 1, \\ u(0) = u(1) = 0, \\ u'(0) = u'(1) = 0. \end{array} \right.$$

Our objective is to recover approximations $u_i$ of $u(x_i)$ from the knowledge of the discrete data $f(x_i)$ on the grid $0 = x_0 < x_1 < \dots < x_{N-1} < x_N = 1$. The problem (121) is approximated by

$$(122) \qquad \left\{ \begin{array}{ll} \delta_x^4 u_j = f(x_j), & 1 \le j \le N-1, \\ u_0 = u_N = 0, \\ u_{x,0} = u_{x,N} = 0. \end{array} \right.$$

In Table 1 we display numerical results for the fourth-order scheme (122). Observe that fourth-order accuracy is achieved in both the maximum and the $l^2$ norms.

<table><caption>TABLE 1. Compact scheme for $u^{(4)} = f$ with exact solution: $u = x^4(x-1)^2$ on $[0, 1]$. We present $e_{max}$ the error in the maximum norm, and $e_2$ the error in the $l^2$ norm.</caption><thead><tr><th>Mesh</th><th>N = 16</th><th>Rate</th><th>N = 32</th><th>Rate</th><th>N = 64</th><th>Rate</th><th>N = 128</th></tr></thead><tbody><tr><td>e<sub>max</sub></td><td>7.8231(-6)</td><td>4.00</td><td>4.8894(-7)</td><td>4.00</td><td>3.0589(-8)</td><td>4.00</td><td>1.9106(-9)</td></tr><tr><td>e<sub>2</sub></td><td>5.6157(-6)</td><td>4.00</td><td>3.5099(-7)</td><td>4.00</td><td>2.1937(-8)</td><td>4.00</td><td>1.3739(-9)</td></tr></tbody></table>

6.1.2. *Nonhomogeneous boundary conditions.* Here we consider a polynomial solution, but with nonhomogeneous values at the two end points,

$$u(x) = x^5, \text{ on } [a, b] = [0, 1].$$

The function $u(x)$ is the solution of the biharmonic problem

$$(123) \qquad \left\{ \begin{array}{ll} u^{(4)} = f, & 0 < x < 1, \\ u(0) = 0, u(1) = 1, \\ u'(0) = 0, u'(1) = 5, \end{array} \right.$$

where the function $f(x)$ is

$$(124) \qquad f(x) = u^{(4)} = 120x.$$

Thus, we resolve numerically

$$(125) \qquad \left\{ \begin{array}{ll} \delta_x^4 u_j = f(x_j), & 1 \le j \le N-1, \\ u_0 = 0, & u_N = 1, \\ u_{x,0} = 0, & u_{x,N} = 5. \end{array} \right.$$

Our purpose is to demonstrate the fourth-order accuracy of the scheme for the case of nonhomogeneous boundary conditions. Indeed, the numerical results reported in Table 2 assess the fourth-order accuracy of the scheme in this case too.
<table><caption>TABLE 2. Compact scheme for $u^{(4)} = f$ with exact solution: $u = x^5$ on [0, 1]. We present $e_{max}$ the error in the maximum norm, and $e_2$ the error in the $l_2$ norm.</caption><thead><tr><th>Mesh</th><th>N = 16</th><th>Rate</th><th>N = 32</th><th>Rate</th><th>N = 64</th><th>Rate</th><th>N = 128</th></tr></thead><tbody><tr><td>e<sub>max</sub></td><td>9.6857(-7)</td><td>3.99</td><td>6.1118(-8)</td><td>4.00</td><td>3.8200(-9)</td><td>3.98</td><td>2.4129(-10)</td></tr><tr><td>e<sub>2</sub></td><td>7.0187(-7)</td><td>4.00</td><td>4.3873(-8)</td><td>4.00</td><td>2.7420(-9)</td><td>4.00</td><td>1.7062(-10)</td></tr></tbody></table>

6.2. **Case 2: Oscillating solutions.** We consider a family of functions defined by

$$ (126) \qquad u_{\varepsilon}(x) = p(x) \sin(1/q_{\varepsilon}(x)), $$

where the polynomial functions $p(x)$ and $q_{\varepsilon}(x)$ are given by

$$ (127) \qquad p(x) = 16x^2(1-x)^2, \quad q_{\varepsilon}(x) = 1/((x-1/2)^2 + \varepsilon), \quad \varepsilon > 0. $$

For small $\varepsilon$ the function $u_{\varepsilon}$ oscillates in the middle of the interval. The parameter $\varepsilon$ serves as a tuning parameter for the frequency of the oscillations.

In Figure 1 we display the functions $u_{\varepsilon}(x)$ corresponding to

$$ (128) \qquad \varepsilon = 7.510^{-2}, \quad \varepsilon = 5.010^{-2}, \quad \varepsilon = 2.510^{-2}. $$

As in Case 1, we consider the approximation

<figure><img src="image_5.png" /><figcaption>FIGURE 1. The oscillating function $x \mapsto 16x^2(1-x)^2 \sin(1/((x-0.5)^2 + \varepsilon))$ for $\varepsilon = 7.5 \cdot 10^{-2}$ (left), $\varepsilon = 5.0 \cdot 10^{-2}$ (center), $\varepsilon = 2.5 \cdot 10^{-2}$ (right).</figcaption></figure>

$$ (129) \qquad \left\{ \begin{array}{l} \delta_x^4 u_j = f_{\varepsilon}(x_j), \quad 1 \le j \le N-1, \\ u_0 = 0, \quad u_N = 0, \\ u_{x,0} = 0, \quad u_{x,N} = 0, \end{array} \right. $$

where the function $f_{\varepsilon}(x)$ is defined by

$$ (130) \qquad f_{\varepsilon}(x) = u_{\varepsilon}^{(4)}(x). $$

The results are reported in Figure 2 on a LogLog scale. In addition, in Table 3 we display the errors for different values of $\varepsilon$ and $N$. The results clearly demonstrate the asymptotic fourth-order convergence in both norms. The magnitude of the errors on relatively coarse grids is remarkably small. Observe that a maximum error of order $10^{-3}$ is obtained for $\varepsilon = 7.5 \cdot 10^{-2}$, $\varepsilon = 5.0 \cdot 10^{-2}$ and $\varepsilon = 2.5 \cdot 10^{-2}$ with $N = 32$, $N = 64$ and $N = 128$, respectively.

6.3. **Case 3: A nonlinear biharmonic equation.** As a final example we consider the nonlinear problem

$$ (131) \qquad \left\{ \begin{array}{l} u^{(4)} - H(u) = f, \quad 0 < x < 1, \\ u(0) = 0, u(1) = 0, \\ u'(0) = 0, u'(1) = 0, \end{array} \right. $$

where $H$ is assumed to be a $k$-lipschitz function.
<figure><img src="image_6.png" /><figcaption>FIGURE 2. Convergence rates for the problem (129) with $\epsilon = 7.5 \cdot 10^{-2}$ (left), $\epsilon = 5.0 \cdot 10^{-2}$ (center), $\epsilon = 2.5 \cdot 10^{-2}$ (right). Circles correspond to the maximum norm and squares correspond to the $l^2$ norm.</figcaption></figure>

<table><caption>TABLE 3. Compact scheme for $u^{(4)} = f_{\epsilon}$ for $\epsilon = 7.5 \cdot 10^{-2}$, $\epsilon = 5.0 \cdot 10^{-2}$, $\epsilon = 2.5 \cdot 10^{-2}$. We present $e_{max}$ the error in the maximum norm, and $e_2$ the error in the $l_2$ norm.</caption><thead><tr><th>Mesh</th><th>N = 16</th><th>Rate</th><th>N = 32</th><th>Rate</th><th>N = 64</th><th>Rate</th><th>N = 128</th><th>Rate</th><th>N = 256</th><th>Rate</th><th>N = 512</th></tr></thead><tbody><tr><td>e<sub>max</sub>, ε = 7.5 10<sup>-2</sup></td><td>7.3136(+0)</td><td>11.67</td><td>2.2523(-3)</td><td>4.62</td><td>9.1729(-5)</td><td>4.09</td><td>5.3870(-6)</td><td>4.02</td><td>3.3124(-7)</td><td>4.01</td><td>2.0501(-8)</td></tr><tr><td>ε = 5.0 10<sup>-2</sup></td><td>2.7999(+3)</td><td>14.45</td><td>1.2469(-1)</td><td>6.72</td><td>1.1839(-3)</td><td>4.27</td><td>6.1165(-5)</td><td>4.06</td><td>3.6542(-6)</td><td>4.04</td><td>2.2279(-7)</td></tr><tr><td>ε = 2.5 10<sup>-2</sup></td><td>2.9713(+5)</td><td>0.26</td><td>2.4744(+5)</td><td>18.75</td><td>5.6056(-1)</td><td>6.75</td><td>5.1981(-3)</td><td>4.44</td><td>2.3902(-4)</td><td>4.04</td><td>1.4556(-6)</td></tr><tr><td>e<sub>2</sub>, ε = 7.5 10<sup>-2</sup></td><td>4.7042(0)</td><td>12.34</td><td>9.0697(-4)</td><td>4.45</td><td>4.1464(-5)</td><td>4.11</td><td>2.4066(-6)</td><td>4.03</td><td>1.4768(-7)</td><td>4.03</td><td>9.0506(-9)</td></tr><tr><td>ε = 5.0 10<sup>-2</sup></td><td>1.7151(+3)</td><td>14.97</td><td>5.3507(-2)</td><td>7.06</td><td>3.9973(-4)</td><td>4.28</td><td>2.0575(-5)</td><td>4.07</td><td>1.2252(-6)</td><td>4.02</td><td>7.5503(-8)</td></tr><tr><td>ε = 2.5 10<sup>-2</sup></td><td>1.8588(+5)</td><td>0.27</td><td>1.5403(+5)</td><td>19.46</td><td>2.1295(-1)</td><td>7.25</td><td>1.3947(-3)</td><td>4.41</td><td>6.5468(-5)</td><td>4.11</td><td>3.8041(-6)</td></tr></tbody></table>

The approximation of (131) is obtained via the (nonlinear) scheme

$$ (132) \qquad \begin{cases} \delta_x^4 \mathbf{u}_j - H(\mathbf{u}_j) = f(x_j), & 1 \le j \le N-1, \\ \mathbf{u}_0 = 0, \quad \mathbf{u}_N = 0, \\ \mathbf{u}_{x,0} = 0, \quad \mathbf{u}_{x,N} = 0. \end{cases} $$

Equation (131) has a unique solution under the sufficient condition

$$ (133) \qquad k < \lambda_{min}, $$

where $\lambda_{min}$ is the smallest eigenvalue of the problem

$$ (134) \qquad \begin{cases} u^{(4)} = \lambda u, \\ u(0) = 0, \quad u(1) = 0, \\ u'(0) = 0, \quad u'(1) = 0. \end{cases} $$

A sufficient condition for (133) to hold is that

$$ (135) \qquad k < \left( \frac{3\pi}{4} \right)^4. $$

In Table 4 we display numerical results for the function $H(u) = 100 \sin^2 u$. Here the right-hand side $f$ is selected as $u_{\epsilon}^{(4)} - H(u)$, with $\epsilon = 5.0 \cdot 10^{-2}$ (see Case 2). Observe the fourth-order accuracy of the scheme.

<table><caption>TABLE 4. Compact scheme for $u^{(4)} - 100 \sin^2 u = f$ with exact solution: $u = u_{\epsilon}(x)$ on $[0, 1]$, with $\epsilon = 5.0 \cdot 10^{-2}$. We present $e_{max}$ the error in the maximum norm, and $e_2$ the error in the $l_2$ norm.</caption><thead><tr><th>Mesh</th><th>N = 64</th><th>Rate</th><th>N = 128</th><th>Rate</th><th>N = 256</th><th>Rate</th><th>N = 512</th></tr></thead><tbody><tr><td>e<sub>max</sub></td><td>1.2391(-3)</td><td>4.27</td><td>6.4175(-5)</td><td>4.06</td><td>3.8378(-6)</td><td>4.04</td><td>2.3401(-7)</td></tr><tr><td>e<sub>2</sub></td><td>4.0399(-4)</td><td>4.28</td><td>2.0839(-5)</td><td>4.07</td><td>1.2417(-6)</td><td>4.03</td><td>7.6153(-8)</td></tr></tbody></table>
**Acknowledgment:** This paper is dedicated to Professor Saul Abarbanel on the occasion of his 80-th birthday.

We would like to thank Professor B. Bialecki of the Colorado School of Mines, who challenged us with providing
a proof for the fourth-order accuracy of the three point biharmonic operator.

REFERENCES

[1] S. Abarbanel, A. Ditkowski, and B. Gustafsson. On error bound of finite difference approximations for partial differential equations. *J. Sci. Comput.*, 15:79–116, 2000.

[2] I. Altas, J. Dym, M. M. Gupta, and R. P Manohar. Multigrid solution of automatically generated high-order discretizations for the biharmonic equation. *SIAM J. Sci. Comput.*, 19:1575–1585, 1998.

[3] M. Ben-Artzi, J-P. Croisille, and D. Fishelov. Convergence of a compact scheme for the pure streamfunction formulation of the unsteady Navier-Stokes system. *SIAM J. Numer. Anal.*, 44:1997–2024, 2006.

[4] M. Ben-Artzi, J-P. Croisille, and D. Fishelov. A fast direct solver for the biharmonic problem in a rectangular grid. *SIAM J. Sci. Comput.*, 31:303–333, 2008.

[5] A. Brüger, B. Gustafsson, P. Lötstedt, and J. Nilsson. High order accurate solution of the incompressible Navier-Stokes equations. *J. Comput. Phys.*, 203:49–71, 2005.

[6] G.F. Carey and W.F. Spotz. High-order compact scheme for the stream-function vorticity equations. *Int. J. Numer. Meth. Eng.*, 38:3497–3512, 1995.

[7] G.F. Carey and W.F. Spotz. Extension of high-order compact schemes to time dependent problems. *Numer. Meth. Part. Diff. Eqs*, 17:657–672, 2001.

[8] M. H. Carpenter, D. Gottlieb, and S. Abarbanel. The stability of numerical boundary treatments for compact high-order finite difference schemes. *J. Comput. Phys.*, 108:272–295, 1993.

[9] W. E and J.-G. Liu. Essentially compact schemes for unsteady viscous incompressible flows. *J. Comput. Phys.*, 126:122–138, 1996.

[10] D. Fishelov, M. Ben-Artzi, and J-P. Croisille. Recent developments in the pure streamfunction formulation of the Navier-Stokes system. *J. Sci. Comput.*, 45:238–258, 2010.

[11] M. M. Gupta, R.P.Manohar, and J.W.Stephenson. Single cell high order scheme for the convection-diffusion equation with variable coefficients. *Int. J. Numer. Meth. Fluids*, 4:641–651, 1984.

[12] B. Gustafsson. The convergence rate for difference approximations to mixed initial boundary value problems. *Math. of Comput.*, 29:396–406, 1975.

[13] B. Gustafsson. The convergence rate for difference approximations to general mixed initial boundary value problems. *SIAM J. Numer. Anal.*, 18:179–190, 1981.

[14] S. K. Lele. Compact finite-difference schemes with spectral-like resolution. *J. Comput. Phys.*, 103:16–42, 1992.

[15] Ming Li and Tao Tang. A compact fourth-order finite difference scheme for unsteady viscous incompressible flows. *J. Sci. Comput.*, 16:29–45, 2001.

[16] J. W. Stephenson. Single cell discretizations of order two and four for biharmonic problems. *J. Comput. Phys.*, 55:65–80, 1984.

DALIA FISHELOV: AFEKA - TEL-AVIV ACADEMIC COLLEGE OF ENGINEERING, 218 BNEI-EFRAIM ST., TEL-AVIV 69107, ISRAEL
*E-mail address:* daliaf@post.tau.ac.il

MATANIA BEN-ARTZI: INSTITUTE OF MATHEMATICS, THE HEBREW UNIVERSITY, JERUSALEM 91904, ISRAEL
*E-mail address:* mbartzi@math.huji.ac.il

JEAN-PIERRE CROISILLE: DEPARTMENT OF MATHEMATICS, LMAM, UMR 7122, UNIVERSITY OF PAUL VERLAINE-METZ, METZ 57045, FRANCE
*E-mail address:* croisil@poncelet.univ-metz.fr