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

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.


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 α ) at inner points and of order O(h α−s ) near the boundary, then if s = 0, 1 the accuracy of the scheme is O(h α ). However, if s ≥ 2 then the overall accuracy the scheme is O(h α−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 α = 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: • 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 δ 2 x , δ 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 onedimensional biharmonic problem. The error (see Theorem 6) is shown to be of fourth-order in the discrete l 2 h 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.

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 We look for a high-order compact approximation to (1). We lay out a uniform grid a = In what follows, we shall use the notion of grid functions. A grid function is a function defined on the discrete grid {x i } N i=0 . We denote grid functions with fraktur letters such as u, v. We have In addition, we denote by u * = (u(x 0 ), u(x 1 ), · · · , 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 2 h the functional space of grid functions. This space is equipped with a scalar product and an associated norm The subspace of grid functions, having zero boundary conditions at x 0 = a and x N = b, is denoted by l 2 h,0 . For grid functions u, v ∈ l 2 h,0 , we have We also define the sup norm for a grid function u We define the difference operators δ x , δ 2 x on grid functions by 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 ). 2 We seek a polynomial of degree 4 The coefficients a 1 , a 2 , a 3 , a 4 of the polynomial are 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 * . Looking at the first equation in (9), we see that a natural candidate for u * This is by definition the Hermitian derivative. If we introduce the three-point operator σ x on grid functions by can rewrite (10) as Observe that a knowledge of (13) u * x,0 = u ′ 0 , u * x,N = u ′ N , is needed in order to solve (12). In addition, we will invoke the following relation (14) σ 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 notatioñ δ 2 x u * , δ 3 x u * , δ 4 x u * for the following operators.
x u * i is an approximation to the fourth-order derivative of u at x i , namely, This approximation, called the discrete biharmonic approximation, is the one-dimensional analog of the Stephenson's scheme [16]. Note that, in the non-periodic setting, boundary values of u x should be given in order to compute δ 4 x at near boundary points x 1 , x N −1 .
3. The operators δ 2 x , δ 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 ∈ R N −1 be the vector corresponding to the grid function u ∈ l 2 h,0 , The vector corresponding to the grid function δ x u is where the matrix K = (K i,m ) 1≤i,m≤N −1 is the skew-symmetric matrix The matrix which corresponds to σ x is P/6, where P is the positive definite (N − 1) × (N − 1) matrix Thus, Equation (12) can be written as where U, U x are the vectors corresponding to u, u x , respectively. In the sequel, we shall also need the matrix representation of δ 2 x . The matrix T , which corresponds to −h 2 δ 2 x , is the (N − 1) × (N − 1) symmetric matrix The matrix P is related to T by Therefore, the matrix which corresponds to the operator σ x (restricted to l 2 h,0 ) is P/6 = I − T /6. 3.2. The eigenvalues and the eigenvectors of δ 2 x . To simplify the notations, we assume from now on that [a, b] = [0, 1], thus N h = 1.
In order to prove the fourth-order accuracy of the scheme, we shall need the eigenvalues and eigenvectors corresponding to δ 2 x , and thus to the matrix T . The eigenvalues of T are We denote the column vectors as Z k ∈ R N −1 and the row vectors as Z j ∈ R N −1 .
It follows that the matrix T satisfies where Λ = diag(λ 1 , · · · , λ N −1 ). The normalized vectors (with respect to (| · | h ), which diagonalize the operator −δ 2 x , are the grid functions z k , which are defined by (28) Equivalently, they may be written as (noting that N h = 1) We have The accuracy of the Hermitian derivative. Now we state a lemma, proved in [3], which indicates the fourth-order accuracy of the Hermitian derivative.

The DBO and its truncation error
As mentioned in Section 2 the approximation δ 4 x u * i , suggested in (16), may serve as approximation to u (4) (x i ). We refer to δ 4 x as the discrete biharmonic operator (DBO). We define Definition 2 (Discrete biharmonic operator (DBO)). Let u ∈ l 2 h be a given grid function. The discrete biharmonic operator is defined by 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 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 σ x δ 4 x at the interior point x i .
where σ x is the Simpson operator defined in (11). The right-hand side can be expressed as Using the definition of u * x , the first term in this expression is The second term in (37) may be written as Thus, in the absence of boundaries, there is a strong connection between δ 4 x and (δ 2 x ) 2 . Explicit estimates for σ x δ 4 x u * i at near boundary points x 1 , x N −1 are given below (see (42)). It results from this representation that σ x δ 4 x actually coincides with the operator (δ 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.
be the grid functions corresponding, respectively, to u, u (4) . Then the DBO δ 4 x satisfies the following accuracy properties: • At near boundary points i = 1 and i = N − 1, the fourth order accuracy of (41) drops to first order, • The error in the energy norm is given by In the above estimates C is a generic constant, that does not depend on u.
Proof. According to (40), we have We now expand (δ 2 x ) 2 in Taylor series. We have We note that for any smooth function where η ∈ (x − h, x + h). 6 Applying δ 2 x to Equation (45) at the interior points Therefore, subtracting this equation from (47), we obtain the estimate (41). Consider now the near boundary point x 1 . We set δ 4 x u * 0 = (u (4) ) * (x 0 ) and then, using the definition of σ x δ 4 x , we have . First, we consider the terms evaluated at x 1 . Recall that x is the Hermitian derivative of u * . Using the boundary values u * 0 = u * x,0 = 0, we have, in view of (33), Inserting the estimates (50), (51) in Equation (49), we obtain Expanding on the term (δ x u * x ) 2 and using again (33), we have For the second term δ 2 x u * 2 we have, as in (51), Inserting the estimates (54), (55) in Equation (53), we obtain, as in (52), Combining the estimates for r 3 and s 3 and inserting them in (48), we obtain (4) ) * i be the truncation error for the fourth-order derivative approximation. We have (58) where v ∈ l 2 h,0 satisfies the estimates established in the previous parts of the lemma (59) The representative matrix of σ x restricted to l 2 h,0 is P/6 = I − T /6. The eigenvalues of P/6 are The matrix norm of its inverse is From (58), (59) and (61), we obtain, Finally, since we get (43).

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 fo the discrete biharmonic operator.

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 u → u x (see (21)), u → δ x u (see (18)) and u → δ 2 x u (see (22)). Let U ∈ R N −1 be the vector corresponding to the grid function u ∈ l 2 h,0 . Therefore, the matrix representation of u → δ 4 x u is The fact that we deal with a boundary value problem, rather than a periodic one, means that P K − KP = 0. However, the commutator is non-zero only at near-boundary points. Using the precise form of this commutator, we get the following proposition. (ii) The symmetric positive definite operator δ 4 x (see (64)) has the matrix form The constants α, β are

5.2.
Error estimate for the one-dimensional Stephenson scheme. In [3] we carried out an error analysis based on the coercivity of δ 4 x . 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 δ 4 x 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 * , e = u − u * , The grid function u * satisfies where r is by definition the truncation error. We later refer to Proposition 3 for estimates on r.
The error e = u − u * satisfies 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 e = u − u * satisfies where C depends only on f .
Proof. Let U, U * ∈ 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 e = u − u * and r, respectively. Using the matrix representation (67), we can write Equations (1)  We therefore have In view of (67) we have that Inverting P SP and multiplying by P R, we have Our goal is to bound the elements of P −1 E by Ch 4 . Note that by Proposition 3 we have Thus, we need to estimate (P SP ) −1 P R. We decompose P SP as follows Note that with L = (6/h 4 )H, Q = 6JJ T , we have 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Λ ′ Z T , where the j-th column of the matrix Z is Z j , as defined in (25). Recall that P = 6I − T (see (23)), and that the eigenvalues λ j of T are given by (24). Therefore, the eigenvalues κ j , 1 ≤ j ≤ N − 1 of P are given by The diagonal matrix Λ ′ contains the eigenvalues of H, which can be written as The element H i,k of the matrix H is .
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, ..., N − 2). The reason is that writing we shall see that (G −1 P R) 1 , (G −1 P R) 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) Recall the elementary inequalities Noting that h = 1/N and using the estimate | sin( ijπ N )| ≤ 1, we obtain Similarly, we have This estimate holds equally for H N −1,N −1 . For the other corner elements of H we have For i, k = 2, ..., N − 2 we have Therefore, the orders of magnitude of the elements of H are bounded by 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) × (N − 1), but it has only four non-zero components at the corner positions, Therefore, QL has only two non-zero rows -the first and the last. The first row is given by  In the above, take "j = even" for "−" and "j = odd" for "+".
We first divide the first and the last row of G by a 1 and annihilate the terms j = 2, ..., 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, ..., N − 2. The result is G 1 , where The same operations on the identity matrix yield the matrix (103) . . .

−b3 a1
−b2 a1 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) The corresponding matrix I 2 (obtained similarly from I 1 ) is (105) . . .
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 (108) 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 Similarly, and using also |b j | ≤ C/h 2 , we get Therefore, the elements of G −1 are bounded by (111) 13 We can now bound the elements of G −1 P R using (109)-(110) and (79). (112) Similarly, we have that |(G −1 P R) N −1 | ≤ Ch 2 .

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. e max = max |u comp − u exact |, 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.
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 < . . . < x N −1 < x N = 1. The problem (121) is approximated by 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.
The function u(x) is the solution of the biharmonic problem where the function f (x) is Thus, we resolve numerically 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 assesses the fourth-order accuracy of the scheme in this case too.
where the polynomial functions p(x) and q ε (x) are given by For small ε the function u ǫ oscillates in the middle of the interval. The parameter ε serves as a tuning parameter for the frequency of the oscillations.
As in Case 1, we consider the approximation where the function f ε (x) is defined by (130) f ε (x) = u (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 ε 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 ε = 7.5 10 −2 , ε = 5.0 10 −2 and ε = 2.5 10 −2 with N = 32, N = 64 and N = 128, respectively. where H is assumed to be a k-lipschitz function.  Table 3. Compact scheme for u (4) = f ε for ε = 7.5 10 −2 , ε = 5.0 10 −2 , ε = 2.5 10 −2 . We present e max the error in the maximum norm, and e 2 the error in the l 2 norm.
The approximation of(131) is obtained via the (nonlinear) scheme A sufficient condition for (133) to hold is that (135) k < 3π 4 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  Table 4. Compact scheme for u (4) − 100 sin 2 u = f with exact solution: u = u ε (x) on [0, 1], with ε = 5.0 10 −2 . We present e max the error in the maximum norm, and e 2 the error in the l 2 norm.