Thierry Gallou
email: gallouet@cmi.univ-mrs.fr

Yannick Sire
email: sire@cmi.univ-mrs.fr.

SOME POSSIBLY DEGENERATE ELLIPTIC PROBLEMS WITH MEASURE DATA AND NON LINEARITY ON THE BOUNDARY

The goal of this paper is to study some possibly degenerate elliptic equation in a bounded domain with a nonlinear boundary condition involving measure data. We investigate two types of problems: the first one deals with the laplacian in a bounded domain with measure supported on the domain and on the boundary. A second one deals with the same type of data but involves a degenerate weight in the equation. In both cases, the nonlinearity under consideration lies on the boundary. For the first problem, we prove an optimal regularity result, whereas for the second one the optimality is not guaranteed but we provide however regularity estimates.

Résumé. Le but de cet article est l'étude d'équations elliptiques pouvant dégénérer, à données mesures, dans un domaine borné, et avec nonlinéarité au bord du domaine. On étudie deux types de problèmes : un premier est une équation elliptique non dégénérée dans un domaine borné avec des données mesures, supportées à la fois à l'intérieur du domaine et sur le bord de celui-ci. On traite dans une deuxième partie un problème elliptique dégénéré. On établit des résultat d'existence et de régularité dans les deux cas. Dans les deux problèmes considérés, la nonlinéarité est au bord du domaine.

The paper is devoted to the study of the two following problems for γ > 1 (1.1)

   ∆u = µ 1 in Ω, u = 0 on Γ 1 , ∂ ν u = µ 2 -|u| γ-1 u on Γ 2 . (1.2)    ∇ • (d(x, ∂Ω) α ∇)u = µ 1 in Ω, u = 0 on Γ 1 , d(x, ∂Ω) α ∂ ν u = µ 2 -|u| γ-1 u on Γ 2 .
where α ∈ (-1, 1) and d(x, ∂Ω) denotes the distance from a point x ∈ Ω to the boundary ∂Ω. Notice that when α = 0, problem (1.2) reduces to (1.1). The weight d(x, ∂Ω) α degenerates at the boundary either by explosion for α < 0 or to 0 for α > 0. We will see that there is a difference, as far as regularity is concerned , between the cases α = 0 and α = 0 and we will deal with these problems separately.

Several works have been devoted to the study of elliptic equations with non smooth data. Particularly, the following equation

-∆u + |u| γ-1 u = f ∈ L 1 in Ω
with Dirichlet boundary conditions has been investigated by several authors starting with the works of Brézis and Strauss [START_REF] Brézis | Semi-linear second-order elliptic equations in L 1[END_REF]. In the case of Ω = R N , we refer the reader to the works [START_REF] Benilan | A semilinear equation in L 1 (R N )[END_REF] and [START_REF] Gallouët | The equation -∆u + |u| α-1 u = f , for 0 ≤ α ≤ 1[END_REF]. In this latter work, the authors investigate the range γ ≤ 1, and prove that a growth condition on f is necessary to ensure existence and uniqueness. The case of the p-laplacian instead of the laplacian has been investigated in [START_REF] Boccardo | Nonlinear elliptic equations in R N without growth restrictions on the data[END_REF].

When f is a measure, the picture is more complicated. In [START_REF] Boccardo | Nonlinear elliptic and parabolic equations involving measure data[END_REF], the authors solved the existence and regularity problem for measure data in smooth domains for Leray-Lions operator of the type -∇ • a(x, u, ∇u). However one can prove existence and regularity if and only if f ∈ L 1 (Ω) + W -2,γ (Ω). This latter condition is equivalent to |f |(A) = 0 for every borelian susbet A of Ω having zero W 2,γ -capacity (here γ is the conjugate exponent of γ).

Under some assumptions on γ, we prove that our problems (1.1)-(1.2) admit a solution and we study their regularity.

The motivation to investigate the degenerate problem (1.2) comes from recent investigations on non local operators. Indeed, the following result has been proved by Caffarelli and Silvestre (see [START_REF] Caffarelli | An extension problem related to the fractional Laplacian[END_REF]): Given s ∈ (0, 1), let α = 1 -2s ∈ (-1, 1). Using variables (x, y) ∈ R n+1

+ := (0, +∞) × R N , the space H s (R N ) coincides with the trace on ∂R N +1 + of H 1 (x α ) := u ∈ H 1 loc (R N +1 + ) : R n+1 + x α u 2 + |∇u| 2 dxdy < +∞ .
In other words, given any function

u ∈ H 1 (x α ) ∩ C(R N +1 + ), v := u| ∂R N +1 + ∈ H s (R N ) and there exists a constant C = C(n, s) > 0 such that v H s (R N ) ≤ C u H 1 (x α )
. So, by a standard density argument (see [START_REF] Piat | Relaxation of degenerate variational integrals[END_REF]), every u ∈ H 1 (x α ) has a well-defined trace v ∈ H s (R N ). Conversely, any v ∈ H s (R N ) is the trace of a function u ∈ H 1 (x α ). In addition, the function u ∈ H 1 (x α ) defined by (1.3) u := arg min

R N +1 + x α |∇w| 2 dx : w| ∂R N +1 + = v solves the PDE (1.4) ∇ • (x α ∇u) = 0 in R N +1 + u = v on ∂R N +1 + . By standard elliptic regularity, u is smooth in R N +1 + . It turns out that x α u x (x, •) converges in H -s (R N ) to a distribution f ∈ H -s (R N ), as x → 0 + i.e. u formally solves (1.5) ∇ • (x α ∇u) = 0 in R N +1 + -x α u x = f on ∂R N +1 + .
Consider the Dirichlet-to-Neumann operator

Γ α : H s (R N ) → H -s (R N ) v → Γ α (v) = f := -x α u x | ∂R N +1 + ,
where u is the solution of (1.3)-(1.5). Given f ∈ H -s (R N ), a function v ∈ H s (R N ) solves the equation

(1.6) 1 d N,s (-∆) s v = f in R N if and only if its lifting u ∈ H 1 (x α ) solves u = v on ∂R N +1 + and (1.7) ∇ • (x α ∇u) = 0 in R N +1 + -x α u x = f on ∂R N +1 + .
Here d N,s is a normalizing constant. Equation (1.6) involves the fractional laplacian , which symbol is a Fourier multiplier |ξ| 2s . We refer the reader to [START_REF] Landkof | Foundations of modern potential theory[END_REF] for a potential-theoretic study of the fractional laplacian and Riesz kernels.

A quick look at the previous development shows that the weight x α represents the distance of a point (x, y) ∈ R N +1 + to the boundary ∂R N +1 + to the power α. It is then natural to consider a somehow "localized" version of it, namely problem (1.2), as a starting point of regularity study of fractional order operators. We postpone to future work the study of the equation

(-∆) s u + |u| γ-1 u = f in R N where f ∈ L 1 loc (R N
). This problem is more challenging since it requires local estimates, as done for instance in [START_REF] Boccardo | Nonlinear elliptic equations in R N without growth restrictions on the data[END_REF] in the case s = 1. More precisely, we would need local estimates independent of the radius R of the following boundary problem

(1.8)    ∇ • (x α ∇u) = µ 1 , in B + R , u = 0, on ∂ + B + R , x α ∂ x u = µ 2 -|u| γ-1 u, on ∂ 0 B + R where B + R = (x, y) ∈ R + × R N , |(x, y)| < R , ∂ + B + R = (x, y) ∈ R + × R N , |(x, y)| = R , ∂ 0 B + R = (0, y), y ∈ R N , |(0, y)| < R .
It has to be noticed that the weight d(x, ∂Ω) α , since α ∈ (-1, 1), is degenerate at the boundary of the smooth domain Ω. However, this function falls into a specific class of functions, namely A 2 weights introduced by Muckenhoupt [START_REF] Muckenhoupt | Weighted norm inequalities for the Hardy maximal function[END_REF]. Indeed, a function w(x) defined on R N is said to be A 2 if the following holds: there exists C > 0 such that sup

B 1 |B| B w(x) dx 1 |B| B w(x) -1 dx ≤ C, for every ball B ⊂ R N . In our context the A 2 condition writes 1 R n R 0 (R -r) α r n-1 dr 1 R n R 0 (R -r) -α r n-1 dr ≤ C,
which is clearly satisfied since α ∈ (-1, 1). The class A 2 enjoys several very nice properties that we will describe further in the paper. First, we introduce the following notion of weak solution, valid for both problems.

Definition 1.1. A function u ∈ W 1,1 (Ω) is a weak solution of (1.2) (or (1.1) taking α = 0) if for all ϕ ∈ C ∞ (R N ) with ϕ = 0 on Γ 1 (1.9) Ω d(x, ∂Ω) α ∇u • ∇ϕ + Ω ϕdµ 1 = Γ 2 ϕdµ 2 - Γ 2 |u| γ-1 uϕ.
We now state our results. Our main result is an existence and regularity result for (1.9). Theorem 1.2. For every µ 1 ∈ M(Ω), µ 2 ∈ M(∂Ω) supported on Ω and Γ 2 respectively, γ > 1, there exists a weak solution u to (1.1) such that u ∈

1<q< N N -1 W 1,q (Ω).
As a consequence, the trace T u of u on ∂Ω satisfies

T u ∈ 1≤q< N N -1 W 1-1 q ,q (∂Ω)
The previous theorem is optimal in the sense that one can construct measures such that u / ∈ W 1,N/N -1 (Ω). Indeed, consider the case µ 1 = δ. Then the solution has the singularity of the Green function, hence is not in W 1,N/N -1 (Ω).

The second one concerns the degenerate case.

Theorem 1.3. For every µ 1 ∈ M(Ω), µ 2 ∈ M(∂Ω) supported on Ω and Γ 2 respectively, γ > 1 and α ∈ (-1, 1), there exists a weak solution u to (1.2) such that

u ∈ 1≤q< 2N +2δ(N -1) 2N -1+δ(N -1) . W 1,q (Ω, d(x, ∂Ω) α )
for some δ > 0 depending on Ω and α. As a consequence, the trace T u of u on ∂Ω satisfies

T u ∈ 1<q< 2N +2δ(N -1)
2N -1+δ(N -1) .

W 1-1+α q ,q (∂Ω)

This second theorem is not optimal, as far as regularity is concerned, because of the degeneracy of the weight close to the boundary.

In order to prove Theorems 1.2-1.3, we proceed in several steps:

• We first approximate the problem with smooth data.

• We obtain estimates with suitable dependence on the data.

• We pass to the limit. This technique is reminiscent of the technique developped in [START_REF] Boccardo | Nonlinear elliptic and parabolic equations involving measure data[END_REF], for instance. In the degenerate case, even if the techniques are similar, this requires new ingredients from the theory of degenerate elliptic equations with A 2 weights, which have been developped in [START_REF] Fabes | The local regularity of solutions of degenerate elliptic equations[END_REF]. For sake of clarity, we will first deal with the problem α = 0, for which an optimal result can be obtained. In a second part, we will consider α = 0.

2. The case α = 0: problem (1.1) 2.1. Regularization of the boundary problem. To do so we start by approximating the problem (1.1) and consider it for smooth sequences {µ n 1 } n≥0 , {µ n 2 } n≥0 ⊂ L ∞ converging to µ 1 and µ 2 in the weak * topology , i.e.

Ω ϕdµ n 1 → Ω ϕdµ 1 , ∀ ϕ ∈ C(Ω) Γ 2 ϕdµ n 2 → Γ 2 ϕdµ 2 , ∀ ϕ ∈ C(∂Ω)
The regularized problem then writes weakly (2.10)

Ω ∇u n • ∇ϕ + Ω ϕµ n 1 = Γ 2 ϕµ n 2 - Γ 2 |u n | γ-1 u n ϕ.
We have the following result, as a consequence of a standard topological degree aregument (see [START_REF] Benilan | A semilinear equation in L 1 (R N )[END_REF]).

Theorem 2.1. Let γ > 1 and µ n 1 , µ n 2 as before. There exists a weak solution H 1 (Ω) of (2.10).

2.2.

Estimates on the regularized problem. We now estimate the solution u n and its gradient ∇u n with a convenient dependence on the data µ n 1 , µ n 2 . Lemma 2.2. Let γ > 1 and u n be a weak solution of (2.10) with µ n 1 , µ n 2 ∈ L ∞ (Γ 2 ). Then for every θ > 1, there exists C > 0 depending on θ and on a bound of the L 1 norm of µ n 1 , µ n 2 such that

(2.11) u n L γ (Γ 2 ) ≤ C

and

(2.12)

Ω |∇u n | 2 (1 + |u n |) θ ≤ C.
Proof. We proof follows the same lines as the one in [START_REF] Boccardo | Nonlinear elliptic equations in R N without growth restrictions on the data[END_REF], taking into account that the nonlinearity is on the boundary. Consider the weak formulation of (2.10), i.e.

(2.13)

Ω ∇u n • ∇ϕ + Ω ϕµ n 1 = Γ 2 ϕµ n 2 - Γ 2 |u n | γ-1 u n ϕ.
for any ϕ defined as before. We consider a suitable test function for the weak formulation (2.13), by considering for θ > 1

(2.14) φ θ (r) = r 0 dt (1+t) θ if r ≥ 0, -φ θ (-r) if r < 0,
Notice that φ θ is bounded in R. We plug the following test function in (2.13) ϕ = φ θ (u n ), which is a suitable test function by Stampacchia theorem. This gives, since φ θ is bounded,

Ω |∇u n | 2 (1 + |u n |) θ + Γ 2 |u n | γ-1 u n φ θ (u n ) ≤ C( µ n 1 L 1 (Γ 2 ) + µ n 2 L 1 (Ω) ).
This gives the desired result since

|u n | γ ≥ C|u n | γ-1 u n φ θ (u n ) for some C > 0 depending on θ.
From the previous estimates, we deduce regularity estimates on the gradient of u n in some Lebesgue space. This is the object of the following lemma.

Lemma 2.3. Let µ n 1 ∈ L ∞ (Ω) and µ n 2 ∈ L ∞ (Γ 2 )
. Consider u n a weak solution of (2.13) such that

Γ 2 |u n | γ ≤ C
for some constant C > 0 and

(2.15)

Ω |∇u n | 2 (1 + |u n |) θ ≤ D
for all θ > 1 and some constant D > 0 depending on θ and on a bound of the L 1 norm of µ n 1 , µ n 2 . Then, for every q ∈ [1, N N -1 ), there exists C depending on C, D such that

u n W 1,q (Ω) ≤ C.
Proof. We write, by Hölder inequality for q ∈ [1, 2)

Ω |∇u n | q = Ω |∇u n | q (1 + |u n |) θq/2 (1 + |u n |) θq/2 ≤ Ω |∇u n | 2 (1 + |u n |) θ q/2 Ω (1 + |u n |) θq 2-q 2-q 2 ≤ C 1 Ω (1 + |u n |) θq 2-q 2-q 2 .
We have for q * = qN N -q

Ω (1 + |u n |) θq 2-q ≤ ε Ω |u n | q * + C 2
for every ε > 0 provided that θ is chosen such that θq 2-q < q * anf for some constant C 2 depending on θ, q, ε. We now use the Sobolev embeddings to estimate the last term: since u n = 0 on Γ 1 ⊂ ∂Ω and q ∈ [1, 2), we have the Sobolev embedding

u n L q * (Ω) ≤ C ∇u n L q (Ω) .
Hence, chossing ε small enough, this yields to u L q * (Ω) ≤ C 3 for some C 3 > 0. By observing that one can take θ arbitrary close to 1, one gets that

∇u n L q (Ω) ≤ C for every q ∈ [1, N N -1 ).
The previous theorem admits the following corollary by just using the trace embedding.

Corollary 2.4. Let µ n 1 , µ n 2 as before. Consider u n a weak solution of (2.13) such that

Γ 2 |u n | γ ≤ C
for some constant C > 0 and

(2.16)

Ω |∇u n | 2 (1 + |u n |) θ ≤ D
for θ > 1 and some constant D > 0 depending on θ and on a bound of the L 1 norm of µ n 1 , µ n 2 . Then the trace on ∂Ω of u n denoted T u n satisfies for every q ∈ (1, N N -1 ), the estimate

T u n W 1-1 q ,q (∂Ω) ≤ C
for some C depending C and D.

2.3. Passage to the limit in the regularized problem. We now come to the proof of Theorem 1.2 by passing to the limit n → +∞ in the weak formulation (2.13). From the previous sections, we have that up to subsequence (2.17)

u n → u weakly in W 1,q (Ω), q ∈ [1, N N -1
).

By continuity of the trace operator, we have

(2.18) T u n → T u weakly in W 1-1 q ,q (∂Ω), q ∈ (1, N N -1
).

Then T u n converges strongly to T u in L q (∂Ω) for some q and then we have up to extraction of a subsequence (2.19) T u n → T u a.e. on Γ 2 .

We have now to pass to the limit in the non linear term of the weak formulation (2.13). To do so, we adopt the strategy of [START_REF] Boccardo | Nonlinear elliptic equations in R N without growth restrictions on the data[END_REF] and have to prove equi-integrability of |u n | γ on Γ 2 . Let t > 0 and define (2.20)

ψ(s) = inf((s -t) + , 1) if s ≥ 0, -ψ(-s) if s < 0, Plugging ϕ = ψ(u n ) in (2.13) and noticing that Ω ∇u n • ∇ψ(u n ) ≥ 0,
we are led to the following estimate

E 0 n,t+1 |u n | γ ≤ E 0 n,t |µ n 2 | + En,t |µ n 1 | where E n,t = (x, t) ∈ Ω × R + | |u n (x)| ≥ t and E 0 n,t = (x, t) ∈ ∂Ω × R + | |T u n (x)| ≥ t .
Hence the quantity |u n | γ is equi-integrable and Vitali's theorem ensures that

Γ 2 |u n | γ-1 u n → Γ 2 |u| γ-1 u.
This proves the desired result.

3. The degenerate case α = 0: problem (1.2)

We now come to the proof of Theorem (1.3). We will not give all the proofs since some of them are almost identical but emphasize the place where we loose regularity.

3.1. Regularization. Consider again smooth sequences {µ n 1 } n≥0 , {µ n 2 } n≥0 converging to µ 1 and µ 2 in the weak * topology in the sense described above. The regularized problem then writes weakly

(3.21) Ω d(x, ∂Ω) α ∇u n • ∇ϕ + Ω ϕµ n 1 = Γ 2 ϕµ n 2 - Γ 2 |u n | γ-1 u n ϕ.
Using the Lax-Milgram theorem in [START_REF] Fabes | The local regularity of solutions of degenerate elliptic equations[END_REF] together with a standard topological degree argument, one gets Theorem 3.1. Let γ > 1 and µ n 1 , µ n 2 as before. There exists a weak solution H 1 (Ω, d(x, ∂Ω) α ) of (3.21).

3.2.

Estimates. The proof of the following lemma is identical to the proof of Lemma (2.2) Lemma 3.2. Let γ > 1 and u n be a weak solution of (3.21) with µ n 1 , µ n 2 ∈ L ∞ (Ω) and L ∞ (Γ 2 ) respectively. Then for every θ > 1, there exists a constant C > 0 depending on a bound of the L 1 norms of µ n 1 , µ n 2 and θ such that

(3.22) u n L γ (Γ 2 ) ≤ C and (3.23) Ω d(x, ∂Ω) α |∇u n | 2 (1 + |u n |) θ ≤ C.
Before stating the lemma concerning the regularity on the gradient of u n , we recall the Sobolev embedding proved in [START_REF] Fabes | The local regularity of solutions of degenerate elliptic equations[END_REF].

Theorem 3.3. (see [START_REF] Fabes | The local regularity of solutions of degenerate elliptic equations[END_REF]) Let Ω be an open set of R N and w an A 2 function; Let u be a function which is 0 on a portion of ∂Ω. Then there exist constants C > 0 and δ > 0 depending on Ω and the weight w such that for all u ∈ H 1 (Ω, w) and all 1 ≤ k ≤ N N -1 + δ the following holds

(3.24) u L 2k (Ω,w) ≤ C ∇u L 2 (Ω,w)
Remark that here the Sobolev exponent is not 2N N -2 but 2N N -1 + 2δ and this is where we loose regularity in the estimates. The constant δ depends on the weight d(x, ∂Ω) α and the domain Ω. Lemma 3.4. Let µ n 1 and µ n 2 as before. Consider u n a weak solution of (3.21) such that

Γ 2 |u n | γ ≤ C
for some constant C > 0 depending on θ and

(3.25) Ω d(x, ∂Ω) α |∇u n | 2 (1 + |u n |) θ ≤ D
for θ > 1 and some constant D > 0 depending on θ and on a bound of the L 1 norms of µ n 1 , µ n 2 . Therefore, we have that u n is bounded (in terms of q, C, D) in W 1,q (Ω, d(x, ∂Ω) α ) for every 1 ≤ q < 2N +2δ(N -1)

2N -1+δ(N -1) .. Proof. We write, by Hölder inequality for q ∈ [1, 2)

Ω d(x, ∂Ω) α |∇u n | q = Ω d(x, ∂Ω) qα/2 |∇u n | q (1 + |u n |) qθ/2 d(x, ∂Ω) α-qα/2 (1 + |u n |) θq/2 ≤ Ω d(x, ∂Ω) α |∇u n | 2 (1 + |u n |) θ q/2 Ω d(x, ∂Ω) α (1 + |u n |) θq 2-q 2-q 2 ≤ D 1 Ω d(x, ∂Ω) α (1 + |u n |) θq 2-q 2-q 2 .
We now use the Sobolev embedding (3.24) to estimate the last term.

Since u n = 0 on Γ 1 ⊂ ∂Ω, we have the Sobolev embedding

u n L 2k (Ω,d(x,∂Ω) α ) ≤ C ∇u n L 2 (Ω,d(x,∂Ω) α )
for all 1 ≤ k ≤ N N -1 + δ for some δ > 0. On the other hand from the inequality (3.25), we have that φ θ/2 (u n ) n is bounded in H 1 (Ω, d(x, ∂Ω) α ). By the previous embedding, the sequence φ θ/2 (u n ) n is also bounded in L 2k (Ω, d(x, ∂Ω) α ) and since, when r is big enough, the function φ θ/2 (r) behaves like r 1-θ/2 , we deduce that

|u n | 1-θ/2 ∈ L 2k (Ω, d(x, ∂Ω) α )
and then

|u n | 2k-kθ ∈ L 1 (Ω, d(x, ∂Ω) α ). Therefore, the sequence {u n } n is bounded in L q (Ω, d(x, ∂Ω) α ) if θq 2 -q = 2k -kθ, i.e. q = 2(2 -θ)k θ + k(2 -θ) > 1
for 1 < θ < 2. Taking θ arbitrarily close to 1, one gets q < 2N + 2δ(N -1) 2N -1 + δ(N -1) .

Since d(x, ∂Ω) α is integrable, one gets a bound for the last term in the previous range of q.

The previous theorem admits the following corollary by just using the trace embedding in [START_REF] Nekvinda | Characterization of traces of the weighted Sobolev space W 1,p (Ω, d M ) on[END_REF]. for θ > 1 and some constant D > 0 depending on θ and on a bound of the L 1 norm of µ n 1 , µ n 2 . Then the trace on ∂Ω of u n denoted T u n is bounded in W 1-1+α q ,q (∂Ω) for every 1 < q < 2N +2δ(N -1) 2N -1+δ(N -1) .

The passage to the limit works exactly as in the previous section and this proves Theorem 1.3.

Extensions of the previous result

Thanks to the approach developped in [START_REF] Boccardo | Nonlinear elliptic and parabolic equations involving measure data[END_REF], it would be possible to extend our results to more general problems, such as (4.27)    -divA(x, ∇u) = µ 1 , in Ω, u = 0, on Γ 1 ⊂ ∂Ω, A(x, ∇u) • ν = µ 2 + h(u), on Γ 2 ⊂ ∂Ω under structural assumptions on A(x, ∇u) (as to be monotone on the gradient term (like for instance the case of the p-laplacian) )and structural assumptions on h.

We refer the interested reader to [START_REF] Boccardo | Nonlinear elliptic and parabolic equations involving measure data[END_REF] for the structural assumptions on the data. Remark that we can assume, thanks to the present work, degeneracy of the operator A(x, .) close to the boundary.

Corollary 3 . 5 . 2 ( 1 +

 3521 Let µ 1 , µ n 2 as before. Consider u n a weak solution of (3.21) such thatΓ 2 |u n | γ ≤ Cfor some constant C > 0 and(3.26) Ω d(x, ∂Ω) α |∇u n | |u n |) θ ≤ D

  Let Ω be a smooth bounded open subset of R N for N ≥ 2. Let {Γ 1 , Γ 2 } be a measurable partition of ∂Ω such that |Γ 1 | > 0. Consider µ 1 ∈ M(Ω) and µ 2 ∈ M(∂Ω), two Radon measures supported on Ω and Γ 2 respectively.

	3. The degenerate case α = 0: problem (1.2)	10
	3.1. Regularization	10
	3.2. Estimates	10
	4. Extensions of the previous result	13
	References	13
	1. Introduction
	1. Introduction	2
	2. The case α = 0: problem (1.1)	6
	2.1. Regularization of the boundary problem	6
	2.2. Estimates on the regularized problem	6
	2.3. Passage to the limit in the regularized problem	9


J Bleyer

D T Pham

P De Buhan

Failure design of high rise concrete panels under fire conditions

Keywords: high rise walls, structural fire safety, yield design, periodic plates and shells, finite elements, non-linear optimization

The present contribution addresses the crucial issue of predicting the failure of high rise reinforced concrete walls subjected to fire loading conditions. A three step procedure is proposed, in which the yield design (limit analysis) method is applied on two separate levels. First, as a way for assessing the local strength properties of the wall, modeled as a plate, as a function of temperature gradient across the wall thickness. Secondly for analyzing the overall stability of the wall in its fire-induced deformed configuration on account of the previously determined strength properties. A comprehensive numerical approach is set up, based on fem-based plate yield design computations associated with non-linear optimization schemes, and a first illustrative application is performed leading to first conclusive results regarding structural fire safety.

INTRODUCTION

Evaluating the ultimate load bearing performance of reinforced concrete members under fire exposure has gained increasing attention in the last decades (see among many others Lie & Celikkol (1991); [START_REF] Franssen | Fire tests and calculation methods for circular concrete columns[END_REF]; [START_REF] El-Fitiany | Assessing the flexural and axial behaviour of reinforced concrete members at elevated temperatures using sectional analysis[END_REF]). Attention has been more particularly focused on the determination of axial force-bending moment interaction diagrams of a reinforced concrete section subjected to a fire induced temperature gradient [START_REF] Caldas | Interaction diagrams for reinforced concrete sections subjected to fire[END_REF]; Law & Gillie (2010)). In this context, the yield design (or limit analysis) approach [START_REF] Chen | Plasticity in reinforced concrete[END_REF]; Salenc ¸on (2013)) has proved being a suitable framework for determining such interaction diagrams in a rigorous way, either under ambient temperature [START_REF] Averbuch | A yield design approach to reinforced concrete structures[END_REF]; [START_REF] Koechlin | Global constitutive model for reinforced concrete plates[END_REF]), or when subjected to a fire-induced temperature gradient (Pham et al. (2015a)).

High rise concrete walls are large size reinforced concrete structures for which the analysis of fire resistance requires a more sophisticated approach than for conventional, i.e. smaller size structures. The sole degradation of the stiffness and strength properties of reinforced concrete due to severe temperature increase, cannot explain as such the collapse of these structures. The purpose of the present contribution is to extend the range of application of the yield design approach (Salenc ¸on ( 2013)) in order to analyse the global stability of high rise walls, taking the geometry changes induced by the thermal loading into account, in addition to the fire-induced decay of material stiffness and strength properties. The main original feature of the present approach lies therefore in the fact that the geometric configuration of the high rise wall on which the yield design approach is to be performed, is not known a priori, but should be determined from a preliminary calculation accounting for geometrically non-linear second order effects.

The whole procedure has already been applied in the particular case when the high rise wall can be modelled as a one dimensional beam-column (Pham et al. (2015b)). It is here extended to rectangular reinforced concrete panels, which can be adequately schematized as a two-dimensional plate (slab) or shell of low curvature in its deformed configuration. An illustrative application of the whole design procedure will be presented, pointing, as a first conclusion, to the very high sensitivity of the wall failure with respect to its width.

OUTLINE OF A THREE-STEP PROCEDURE

Due to the thermal deformations induced by the fire exposure on one of its face, high rise panels exhibit potentially important out-of-plane (horizontal) displacements, which lead to an eccentricity of the gravity load (self-weight) with respect to the initial plane (Figure 1). As a consequence, bending moments are generated in the wall in addition to the preexisting compressive membrane forces distribution, which is usually known as a second order (or P-delta) effect (see for instance the classical textbook by [START_REF] Bazant | Stability of Structures: Elastic, Inelastic, Fracture and Damage Theories[END_REF]). As the eccentricity increases, the moment due to self-weight eccentricity also increases, thus subjecting the wall to still higher bending moments and associated curvature deformations. Simultaneously, but independently, elevated temperature leads to degrading the mechanical properties of constituent materials. It is the combined effect of fire-induced material strength degradation and development of bending moments due to geometry change, which may trigger an overall failure of the structure, even before the occurrence of any buckling phenomenon.

Following such a line of reasoning, this contribution will implement a three-step computational procedure, summarized in the flowchart of Figure 2 below.

• Step n • 1. Preliminary determination of the wall deformed configuration. This step consists in evaluating the equilibrium configuration of the wall under the combined action of thermal gradient (due to a temperature increase of one of the wall face from 20C to 1050 C) and self-weight, through the use of a classical fem-based thermo-elastic computational tool accommodating for second order effects.

•

Step n • 2. Determination of temperature dependent interaction diagrams. The objective of this phase, which is completely independent from the first one, is to determine the membrane forces-bending moments yield strength capacities of any wall cross-section as a function of the prescribed temperature gradient, as it has been done for heterogeneous periodic plates [START_REF] Bleyer | A computational homogenization approach for the yield design of periodic thin plates-Part I: Construction of the macroscopic strength criterion[END_REF]).

•

Step n • 3. Yield analysis and design of the wall in its deformed configuration determined in step n • 1, on account of its reduced strength properties evaluated in step n • 2. This will be performed by means of a finite element formulation of the yield design plate (or more precisely shallow shell) problem [START_REF] Bleyer | On the performance of non-conforming finite elements for the upper bound limit analysis of plates[END_REF]), making use of the yield design static and kinematic approaches, thereby leading to the formulation of a non linear optimization problem. One of the important differences with the case where the panel is modelled as a one-dimensional beam (see Pham et al. (2015b)), is that steps 2 and 3 will not be treated separately but simultaneously within the same global optimization procedure.

PREDICTING THE THERMAL-INDUCED EQUILIBRIUM CONFIGURATION OF THE PANEL

Prior to the application of thermal loading, the panel is modelled as a plane rectangular vertical plate of height H and width L contained in the Oxz-plane as shown in Figure 3, where the different boundary conditions corresponding to situations encountered in practice, are specified. In short, the plate is simply supported along its four edges (free rotation along these edges), while the other prescribed kinematic conditions indicate that the panel is free to move both horizontally and vertically along its upper edge, horizontally and vertically upwards along the lower edge, and vertically and horizontally to the right (resp. to the left) along the left (resp. right) hand lateral side.

The wall being subjected to its own weight only, characterized by a uniform in-plane surface vertical density w, the thermal loading is applied by progressively increasing a) Thermo-elastic strains as well as rotations remain small, but the change of geometry produced by horizontal out-of-plane displacements is taken into account in an explicit way, through the emergence of bending moment distributions and related elastic bending curvatures due to self weight eccentricity. b) The influence of high temperature on the material stiffness degradation is accounted for through experimentally derived curves, giving for instance the variation of the concrete Youngs modulus as a function of temperature (EN1992-1-2 ( 2004)).

The finite element simulation has been carried out with the following sets of data:

L = 12 m and L = 60 m, H = 12 m, h = 0.15 m, w = 3.68 kN/m 2 , Ec = 19.2 GPa, Θmax = 1050 • C (1)
where h is the wall thickness and Ec the concrete Youngs modulus at ambient temperature, one of the walls face being subjected to an ISO 834 fire (EN1991-1-2 ( 2002)).

The result of the non-linear finite element simulation for L = 12 m is represented in Figure 4, in the form of the contours of the out-of-plane displacements experienced by the wall for a 120mn fire exposure (values in meters). The displacement profile of the side view has been magnified by a factor 10. As expected, the displacement reaches its maximum value of 0.30m in the central part of the panel.

TEMPERATURE DEPENDENT YIELD STRENGTH CAPACITIES OF THE PANEL

Following the same kind of upscaling procedure as that adopted for determining the interaction diagram of a reinforced Referring to the lower bound static approach of yield design for this particular problem, plane stress fields of the form:

σ = σ αβ (ξ 3 )e α ⊗ e β α, β = 1, 2
(2) are considered in the concrete material, while each of the four reinforcing bars, modelled as 1D structural elements embedded in the concrete material, undergoes a purely axial force denoted by n ± α,β (see Figure 5).

The concrete material obeys a Mohr-Coulomb tension cutoff strength condition which can be written as:

F (σ αβ ) ≤ 0 ⇔      -2fc ≤ σ 11 + σ 22 ≤ 0 σ 11 σ 22 -σ 2 12 ≤ 0 (σ 11 + fc)(σ 22 + fc) -σ 2 12 ≤ 0 ( 3
)
where fc is the uniaxial compressive strength of the concrete material, while each reinforcement obeys the following strength condition:

|n ± α,β | ≤ Afy ( 4
)
where fy is the reinforcing material (steel) uniaxial yield strength and A = πφ 2 /4 the cross sectional area of each individual bar of diameter φ.

Under such conditions, it can be proved that the panel strength criterion expressed in terms of in plane membrane forces and bending moments, is defined as:

F hom (N αβ , M αβ ) ≤ 0 ⇔                                        ∃σ αβ (ξ 3 ), n ±
α,β such that:

N = h/2 -h/2 σ αβ (ξ 3 )e α ⊗ e β dξ 3 + n + 1 + n - 1 d e 1 ⊗ e 1 + n + 2 + n - 2 d e 2 ⊗ e 2 M = - h/2 -h/2 ξ 3 σ αβ (ξ 3 )e α ⊗ e β dξ 3 - h/2 -a d (n + 1 -n - 1 )e 1 ⊗ e 1 + (n + 2 -n - 2 )e 2 ⊗ e 2 subject to F σ αβ (ξ 3 ) ≤ 0 ∀ξ 3 ∈ [-h/2; h/2] and |n ± α,β | ≤ Afy
(5) When a thermal loading is applied, characterized by a temperature increase profile θ(ξ 3 ) across the wall thickness, the material yield strength properties are modified through the introduction of non dimensional reduction coefficients, the variation of which with temperature, is shown in Figure 6: fc(ξ 3 ) = kc(θ(ξ 3 ))fc and fy(ξ 3 ) = ky(θ(ξ 3 ))fy (6)

The numerical yield design calculations presented in the next section have been performed with the following geometrical characteristics and strength properties at ambient temperature (θ = 20 • C): d = 10 cm, a = 3 cm, φ = 6 mm, fy = 500 MPa, fc = 32 MPa (7)

STABILITY ANALYSIS OF THE PANEL IN ITS DEFORMED CONFIGURATION: A FINITE ELEMENT APPROACH

The stability analysis of the wall is assessed using finite element-based yield design static and kinematic approaches in the deformed configuration by estimating a safety factor λ + corresponding to the non-dimensional multiplicative factor of the loading at which the wall collapses in the sense of yield design. The wall is discretized as an assembly of triangular facets and each triangular finite element is a combination of a plate bending and a membrane finite element.

Implementation of the static approach

In each finite element, a linear interpolation of the in-plane membrane tensor N and a quadratic interpolation of the bending moment tensore M are considered, these fields being a priori discontinuous from one element to another. The local equilibrium inside each facet, the continuity conditions between two adjacent facets and the boundary conditions can be written as a global linear equality of the form HΣ + F = 0 where Σ = {N e αβ , M e αβ , V e α } T is the generalized vector of static unknowns and F the global vector of external loads (weight). The panel strength criterion (5) cannot be written as such in practice. A lower bound approximation F n int (N αβ , M αβ ) ≤ 0 is thus obtained by considering piecewise constant stress fields σ αβ in n different layers of the wall thickness. This implicit discretization also enables to account for the evolution of the concrete strength in the wall thickness. Finally, for a given discretization, a lower bound estimate λ lb of the safety factor λ + can be obtained by solving the following non-linear convex optimization problem:

λ + ≥ λ lb = max Σ λ s.t. HΣ + λF = 0 F n int (N e αβ , M e αβ ) ≤ 0 ∀e (8)
Due to the particular form of the tension cut-off Mohr-Coulomb criterion in plane stress (3), the non-linear constraints

F n int (N e αβ , M e
αβ ) ≤ 0 can be expressed as a combination of linear constraints and second-order conic (SOC) constraints, so that the previous optimization problem can be formulated as a second-order cone programming (SOCP) problem for which very efficient interior-point solvers, like the MOSEK software package Mosek (2002), are available.

Implementation of the kinematic approach

The kinematic approach is based upon the consideration of a virtual shell kinematics consisting of a virtual velocity field u and an angular velocity field θ which are linked in each facet by a Love-Kirchhoff kinematic constraint, due to the fact that the shear resistance of the wall is assumed to be infinite. The stability of the wall in the sense of the yield design kinematic approach is then assured if :

∀(u, θ) kinematically admissible W ext (u, θ) ≤ Wrm(u, θ) = Ω π hom (δ, χ)dΩ (9) Prepared using ECMAuth.cls
where W ext denotes the work of external loads, Wrm denotes the maximum resisting work computed using the wall strength criterion support function defined as :

π hom (δ, χ) = sup F hom (N ,M )≤0 {N : δ + M : χ} (10)
δ being the in-plane membrane strain tensor and χ the bending curvature strain associated to the virtual velocity fields.

A quadratic interpolation for u and a linear interpolation for θ are considered in each finite element, these fields being a priori discontinuous from one element to another, the contribution of the potential discontinuities ("yield lines") being taken into account when computing the maximum resisting work.

The strain-velocity compatibility conditions are written as

δ = {δ e αβ } = B δ U et χ = {χ e αβ } = BχU
where U denotes all kinematic unknowns. The work of external loads is then given by W ext = F T • U. Finally, the maximum resisting work is approximated by a numerical integration quadrature in each triangle and on each edge of the mesh.

Due to the specific form of the strength criterion (5), it can be shown that its support function is given by the sum of the concrete contribution and the reinforcing bars contribution: In practice, the remaining integral is estimated by excess using a n-points trapezoidal quadrature rule.

hom (δ, χ) = h/2 -h/2 Π(δ -ξ 3 χ; ξ 3 )dξ 3
Finally, an upper bound estimate λ ub of the safety factor is obtained by solving the following non-linear convex optimization problem:

λ + ≤ λ ub = min U e n i=1 ω e,i Π(d i ; ξ i 3 ) + 4 k=1 N k y |d k a | s.t. F T • U = 1 d i = δ -ξ i 3 χ = (B δ -ξ i 3 Bχ)U ∀i d 1 a = δ 11 -aχ 11 , d 2 a = δ 22 -aχ 22 d 3 a = δ 11 + aχ 11 , d 4 a = δ 22 + aχ 22 (12)
where ω e,i are integration weights and ξ i 3 are integration points in the thickness direction (note that the discontinuity term has not been written for the sake of concision). Again, due to the particular form of the Mohr-Coulomb criterion, its support function can be expressed using secondorder cone constraints, so that the previous problem is a SOCP problem as well and can be solved using the same solver.

RESULTS AND COMMENTS

The finite element computations have been performed on a half wall (with appropriate symmetry conditions) discretized into 900 triangular finite elements and using n = 10 integration points for the strength criterion and support function computation. The optimization problems have been solved by MOSEK in a few minutes only.

The values of the computed safety factor estimates have been reported in Table 1 for the considered geometry. One can observe that the relative gap between the lower bound static approach and the upper bound kinematic approach is of 12% which can be considered as a quite good estimation from an engineering point of view. One interesting aspect of the yield design, in particular the kinematic approach, is that an optimal collapse mechanism is obtained and can provide valuable information for the designer. For example, the collapse mechanism represented in Figure 7 shows that most of the collapse takes place in the bottom part of the wall, the out-of-plane collapse velocity reaching its maximum at approximately 3m from the ground. This is confirmed by the contour plot of the support function of Figure 8 which shows that the highest resisting work ("dissipation") takes place in the bottom corners and in a central region near the bottom of the wall.

The above described three-step procedure has been then applied to a high rise wall having exactly the same geometrical, loading and strength characteristics as the previous one, except that the width was equal to 60 m instead of 12 m. This example may not seem very realistic but serves as a limit case example to better illustrate the influence of the lateral boundary conditions on the collapse mechanism and the safety factor. Figure 9 shows the contours of the out-of-plane displacements produced by the same thermal gradient. It is worth noting that the maximum value is almost equal to 90cm, that is three times the corresponding value of the square panel (Figure 4). The associated lower and upper bound estimates for the stability factor, which remain relatively close to each other (relative difference of 11%), are reported in the bottom line of Table 1. They show a considerable reduction of the stability factor estimates, which are decreased by a factor almost equal to 14. Figure 10 displays the failure mechanism resulting from the finite element optimization of the upper bound kinematic approach, where it could be observed that it involves an in-plane horizontal velocity jump along the vertical central line. The corresponding distribution of the support function (dissipation or maximum resisting work per unit area) is displayed in Figure 11. Note that the contribution of the different velocity jump lines between adjacent elements, and notably the above mentioned central vertical velocity jump line, is not accounted for in the latter figure.

CONCLUDING REMARKS

It has been shown in this paper, that the combination of a yield design homogenization method for periodic plates, with a fem-based numerical approach, involving non-linear optimization schemes implemented through dedicated solvers, makes it possible to assess the stability of a high rise reinforced concrete panel under fire loading in a rather quick way, leading to fairly accurate bounds on the panel overall stability factor. The computational procedure can be easily generalized to other situations where the geometrical, material, as well as thermal loading characteristics can be varied, allowing for parametric studies.

One of the most important input data necessary to perform such a stability analysis, is to obtain a reliable evaluation of the wall deformed configuration due to thermal loading. Indeed, as it has been proved on the illustrative example investigated in this paper, the amplitude of the out-of-plane displacements has a considerable influence on the wall stability, which is highly sensitive to the lateral boundary conditions.

In this respect, the current evaluation method, based on a second order thermo-elastic analysis, probably provides a conservative estimate for the wall deflection, since only material stiffness, and not yield, properties are accounted for in the analysis. This point certainly requires further investigations in the near future, without it being necessary to resort to sophisticated and complex incremental computational methods, hardly compatible with an engineering design approach.

Fig. 1 .

 1 Fig. 1. Schematic view of a high rise wall subject to fire loading: (left) initial and (right) deformed configurations

Fig. 2 .Fig. 3 .

 23 Fig. 2. A three-step evaluation procedure for the yield design of fire loaded high rise vertical walls

Fig. 4 .

 4 Fig. 4. Contours of the panel out-of-plane displacements for L = 12 m: face and side views.

Fig. 5 .

 5 Fig. 5. Yield design auxiliary problem for determining the reinforced concrete panel strength capacities

Fig. 6 .

 6 Fig. 6. Evolution of concrete and reinforcement strength reduction factors with temperature increase (EN1992-1-2 (2004))

  +Ny(a) (|δ 11 -aχ 11 | + |δ 22 -aχ 22 |) +Ny(-a) (|δ 11 + aχ 11 | + |δ 22 + aχ 22 |) 3 ), Π(d; ξ 3 ) = sup F (σ αβ )≤0 {σ : d} is the support function of the Mohr-Coulomb strength criterion at ξ 3 .

Fig. 7 .Fig. 8 .

 78 Fig. 7. Optimal failure mechanism obtained in the kinematic approach for L = 12 m (contours represent the relative amplitude of the Uy velocity)

Fig. 9 .

 9 Fig. 9. Contours of the panel out-of-plane displacements for L = 60 m: face view.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Optimal failure mechanism obtained in the kinematic approach for L = 60 m (contours represent the relative amplitude of the Uy velocity)

Table 1 .

 1 Failure design of high rise concrete panels under fire loading Values of the obtained safety factor estimates for different geometries

	Geometry	λ lb	λ ub	Relative gap
	12 m×12 m 32.11 36.23	12.1%
	60 m×12 m	2.39	2.66	10.7%

Prepared using ECMAuth.cls


D Fishelov

Matania Ben-Artzi

J.-P Croisille
email: croisil@poncelet.univ-metz.fr

Jean-Pierre Croisille

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

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

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 [START_REF] Fishelov | Recent developments in the pure streamfunction formulation of the Navier-Stokes system[END_REF] 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) [START_REF] Altas | Multigrid solution of automatically generated high-order discretizations for the biharmonic equation[END_REF][START_REF] Brüger | High order accurate solution of the incompressible Navier-Stokes equations[END_REF][START_REF] Carey | High-order compact scheme for the stream-function vorticity equations[END_REF][START_REF] Carey | Extension of high-order compact schemes to time dependent problems[END_REF][START_REF] Carpenter | The stability of numerical boundary treatments for compact high-order finite difference schemes[END_REF][START_REF] Liu | Essentially compact schemes for unsteady viscous incompressible flows[END_REF][START_REF] Gupta | Single cell high order scheme for the convection-diffusion equation with variable coefficients[END_REF][START_REF] Lele | Compact finite-difference schemes with spectral-like resolution[END_REF][START_REF] Li | A compact fourth-order finite difference scheme for unsteady viscous incompressible flows[END_REF][START_REF] Stephenson | Single cell discretizations of order two and four for biharmonic problems[END_REF].

A convergence analysis was performed in [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value prolems[END_REF][START_REF] Gustafsson | The convergence rate for difference approximations to general mixed initial boundary value prolems[END_REF][START_REF] Abarbanel | On error bound of finite difference approximations for partial differential equations[END_REF] in cases where the accuracy of the scheme deteriorates near the boundary. In particular, in [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value prolems[END_REF] and [START_REF] Gustafsson | The convergence rate for difference approximations to general mixed initial boundary value prolems[END_REF] 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 [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value prolems[END_REF][START_REF] Gustafsson | The convergence rate for difference approximations to general mixed initial boundary value prolems[END_REF] 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 [START_REF] Abarbanel | On error bound of finite difference approximations for partial differential equations[END_REF] 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 [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value prolems[END_REF][START_REF] Gustafsson | The convergence rate for difference approximations to general mixed initial boundary value prolems[END_REF][START_REF] Abarbanel | On error bound of finite difference approximations for partial differential equations[END_REF]. 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

(1) u (4) (x) = f (x), a < x < b, u(a) = 0, u(b) = 0, u ′ (a) = 0, u ′ (b) = 0.
We look for a high-order compact approximation to [START_REF] Abarbanel | On error bound of finite difference approximations for partial differential equations[END_REF]. We lay out a uniform grid a =

x 0 < x 1 < ... < x N -1 < x N = b. Here x i = ih for 0 ≤ i ≤ 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 } N i=0 . We denote grid functions with fraktur letters such as u, v. We have

(2) u = (u(x 0 ), u(x 1 ), • • • , u(x N -1 ), u(x N )).
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

(3) (u, v) h = h N i=0 u(x i )v(x i ), |u| h = (u, u) 1/2 h .
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

(4) (u, v) h = h N -1 i=1 u(x i )v(x i ).
We also define the sup norm for a grid function u

(5) |u| ∞ = max 0≤i≤N |u(x i )|.
We define the difference operators δ x , δ 2 x on grid functions by ( 6)

δ x u i = u i+1 -u i-1 2h , 1 ≤ i ≤ N -1, (7) δ 2 x u i = u i+1 -2u i + u i-1 h 2 , 1 ≤ i ≤ 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) 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)                    a 1 = 3 4h (u * i+1 -u * i-1 ) -( 1 4 u * x,i+1 + 1 4 u * x,i-1 ), a 2 = 1 h 2 (u * i+1 + u * i-1 -2u * i ) -1 4h (u * x,i+1 -u * x,i-1 ) = δ 2 x u * i -1 2 (δ x u * x ) i , a 3 = -1 4h 3 (u * i+1 -u * i-1 ) + 1 4h 2 (u * x,i+1 + u * x,i-1 ), a 4 = -1 2h 4 (u * i+1 + u * i-1 -2u * i ) + 1 4h 3 (u * x,i+1 -u * x,i-1 ) = 1 2h 2 (δ x u * x ) i -δ 2 x u * i .
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 N -1 i=1 have to be evaluated in terms of {u * i } N i=0 . Looking at the first equation in [START_REF] Liu | Essentially compact schemes for unsteady viscous incompressible flows[END_REF], we see that a natural candidate for u *

x,i is u *

x,i = a 1 . This yields

u * x,i = 3 4h (u * i+1 -u * i-1 ) -( 1 4 u * x,i+1 + 1 4 u * x,i-1 ), or equivalently (10) 1 6 u * x,i + 2 3 u * x,i + 1 6 u * x,i+1 = δ x u * i .
This is by definition the Hermitian derivative. If we introduce the three-point operator σ x on grid functions by ( 11)

σ x v i = 1 6 v i-1 + 2 3 v i + 1 6 v i+1 , 1 ≤ i ≤ N -1,
can rewrite (10) as ( 12)

σ x u * x,i = δ x u * i , 1 ≤ i ≤ N -1.
Observe that a knowledge of [START_REF] Gustafsson | The convergence rate for difference approximations to general mixed initial boundary value prolems[END_REF] u * x,0 = u ′ 0 , u * x,N = u ′ N , is needed in order to solve [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value prolems[END_REF]. In addition, we will invoke the following relation ( 14)

σ x = I + h 2 6 δ 2 x .
Natural approximations to u ′′ (x i ), u ′′′ (x i ), u ′′′′ (x i ) are 2a 2 , 6a 3 , 24a 4 , respectively (see [START_REF] Liu | Essentially compact schemes for unsteady viscous incompressible flows[END_REF]). We use the notation δ2 x u * , δ 3

x u * , δ 4 x u * for the following operators.

(

)            δ2 x u * i = 2a 2 = 2δ 2 x u * i -(δ x u * x ) i , δ 3 x u * i = 6a 3 = 1 h 2 (u * x,i+1 + u * x,i-1 -2u * x,i ) = (δ 2 x u * x ) i , δ 4 x u * i = 24a 4 = 12 h 2 (δ x u * x ) i -δ 2 x u * i . This suggests that δ 4 15
x u * i is an approximation to the fourth-order derivative of u at x i , namely, (16)

δ 4 x u * i = 12 h 2 (δ x u * x ) i -δ 2 x u * i .
This approximation, called the discrete biharmonic approximation, is the one-dimensional analog of the Stephenson's scheme [START_REF] Stephenson | Single cell discretizations of order two and four for biharmonic problems[END_REF]. 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 , (17

) U = [u 1 , • • • , u N -1 ] T .
The vector corresponding to the grid function

δ x u is (18) 1 2h KU,
where the matrix K = (K i,m ) 1≤i,m≤N -1 is the skew-symmetric matrix

(19) K =        0 1 0 . . . 0 -1 0 1 . . . 0 . . . . . . . . . . . . . . . 0 . . . -1 0 1 0 . . . 0 -1 0       
.

The matrix which corresponds to σ x is P/6, where P is the positive definite (N -1) × (N -1) matrix where U, U x are the vectors corresponding to u, u x , respectively.

(20) P =        4
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

(22) T =        2 -1 0 . . . 0 -1 2 -1 . . . 0 . . . . . . . . . . . . . . . 0 . . . -1 2 -1 0 . . . 0 -1 2       
.

The matrix P is related to T by (23) P = 6I -T.

Therefore, the matrix which corresponds to the operator σ x (restricted to l 2 h,0 ) is P/6 = I -T /6.

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 (24)

λ j = 4 sin 2 ( jπ 2N ), j = 1, • • • , N -1
and the corresponding normalized eigenvectors are

Z k = (Z 1k , • • • , Z N -1,k ) T (with respect to the Euclidean norm in R N -1 ), where (25)
Z jk = 2 N 1/2 sin kjπ N , 1 ≤ k, j ≤ N -1.
We denote the column vectors as Z k ∈ R N -1 and the row vectors as

Z j ∈ R N -1 . The matrix Z = (Z jk ) 1≤j,k≤N -1 ∈ M N -1 (R) is an orthogonal positive-definite matrix. Thus, (26)
Z 2 = ZZ T = I N -1 .
It follows that the matrix T satisfies

(27) T = ZΛZ T , 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)

z jk = Z jk /h 1/2 .
Equivalently, they may be written as (noting that N h = 1) (29)

z jk = √ 2 sin kjπ N , 1 ≤ k, j ≤ N -1.
We have

(30)    z jk = √ 2 sin(j kπ N ), j = 1, • • • , N -1, k = 1, ..., N -1 z 0k = 0, z N k = 0, -δ 2 x z k = λk z k , λk = 4 h 2 sin 2 ( kπ 2N ), k = 1, • • • , N -1. 3.3.
The accuracy of the Hermitian derivative. Now we state a lemma, proved in [START_REF] Ben-Artzi | Convergence of a compact scheme for the pure streamfunction formulation of the unsteady Navier-Stokes system[END_REF], 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 ≤ i ≤ N by (31) (σ x u x ) i = (δ x u * ) i , 1 ≤ i ≤ N -1 and (32) (u x ) 0 = (u ′ ) * (x 0 ), (u x ) N = (u ′ ) * (x N ), has a truncation error u x -(u ′ ) * of order O(h 4 ). More precisely, (33) |u x -(u ′ ) * | ∞ ≤ Ch 4 u (5) L ∞ .

The DBO and its truncation error

As mentioned in Section 2 the approximation δ 4 x u * i , suggested in [START_REF] Stephenson | Single cell discretizations of order two and four for biharmonic problems[END_REF], 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

(34) δ 4 x u i = 12 h 2 (δ x u x,i -δ 2 x u i ), 1 ≤ i ≤ N -1.
Here u x is the Hermitian derivative of u satisfying [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value prolems[END_REF] 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)              (a) δ 4 x u i = f (x i ) 1 ≤ i ≤ N -1, (b) 1 6 u x,i-1 + 2 3 u x,i + 1 6 u x,i+1 = δ x u i , 1 ≤ i ≤ N -1, (c) u 0 = 0, u N = 0, u x,0 = 0, u x,N = 0.
The scheme in (35) is the one-dimensional restriction of the scheme proposed by Stephenson in [START_REF] Stephenson | Single cell discretizations of order two and four for biharmonic problems[END_REF]. 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 .

(36)

σ x δ 4 x u * i = 1 6 δ 4 x u * i-1 + 2 3 δ 4 x u * i + 1 6 δ 4 x u * i+1 , 2 ≤ i ≤ N -2,
where σ x is the Simpson operator defined in [START_REF] Gupta | Single cell high order scheme for the convection-diffusion equation with variable coefficients[END_REF]. The right-hand side can be expressed as

(37) 12 h 2 1 6 δ x u * x,i-1 + 2 3 δ x u * x,i + 1 6 δ x u * x,i+1 - 1 6 δ 2 x u * i-1 + 2 3 δ 2 x u * i + 1 6 δ 2 x u * i+1 .
Using the definition of u * x , the first term in this expression is (38)

1 6 δ x u * x,i-1 + 2 3 δ x u * x,i + 1 6 δ x u * x,i+1 = σ x δ x u * x,i = δ x σ x u * x,i = δ x δ x u * i = 1 4h 2 (u * i+2 -2u * i + u * i-2 ), 2 ≤ i ≤ N -2.
The second term in (37) may be written as (39)

1 6 δ 2 x u * i-1 + 2 3 δ 2 x u * i + 1 6 δ 2 x u * i+1 = 1 6h 2 (u * i-2 + 2u * i-1 -6u * i + 2u * i+1 + u * i+2 ), 2 ≤ i ≤ N -2. Therefore, inserting (38)-(39) in (36), we have (40) σ x δ 4 x u * i = 1 h 4 (u * i-2 -4u * i-1 + 6u * i -4u * i+1 + u * i+2 ) = δ 2 x δ 2 x u * i , 2 ≤ i ≤ N -2.
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

x ) 2 at points x i , 2 ≤ i ≤ 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 δ 4

x satisfies the following accuracy properties:

• (41) |σ x δ 4 x u * i -σ x (u (4) ) * (x i )| ≤ Ch 4 u (8) L ∞ , 2 ≤ i ≤ N -2.
• At near boundary points i = 1 and i = N -1, the fourth order accuracy of (41) drops to first order, (42)

|σ x δ 4 x u * 1 -σ x (u (4) ) * (x 1 )| ≤ Ch u (5)
L ∞ , with a similar estimate for i = N -1.

• The error in the energy norm is given by

(43) |δ 4 x u * -(u (4) ) * | h ≤ Ch 3/2 ( u (5) L ∞ + u (8) L ∞ ).
In the above estimates C is a generic constant, that does not depend on u.

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

σ x δ 4 x u * i = (δ 2 x ) 2 u * i , i = 2, ..., N -2.
We now expand (δ 2 x ) 2 in Taylor series. We have

(45) δ 2 x u i = u ′′ (x i ) + h 2 12 u (4) (x i ) + h 4 360 u (6) (ξ i ), 1 ≤ i ≤ N -1,
where

ξ i ∈ (x i -h, x i + h).
We note that for any smooth function v(x), x ∈ [a, b], we have for a

+ h < x < b -h (46) v(x + h) -2v(x) + v(x -h) h 2 = v ′′ (η),
where η ∈ (xh, x + h).

Applying δ 2

x to Equation (45) at the interior points x i , 2 ≤ i ≤ N -2, we obtain (47)

σ x δ 4 x u * i = (δ 2 x ) 2 u * i = u (4) (x i ) + h 2 6 u (6) (x i ) + p i , |p i | ≤ C 1 h 4 u (8) L ∞ , 2 ≤ i ≤ N -2.
On the other hand, σ x (u (4) ) * (x i ) may be expanded around x i , 2 ≤ i ≤ N -2, as follows.

σ x (u (4) ) * (x i ) = (I + h 2 6 δ 2 x )(u (4) ) * (x i ) = u (4) (x i ) + h 2 6 u (6) (x i ) + q i , |q i | ≤ C 2 h 4 u (8) L ∞ , 2 ≤ i ≤ N -2.
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 (48)

σ x δ 4 x u * 1 -σ x (u (4) ) * (x 1 ) = 2 3 δ 4 x u * 1 + 1 6 δ 4 x u * 2 - 2 3 (u (4) )(x 1 ) + 1 6 (u (4) )(x 2 ) = 2 3 δ 4 x u * 1 -(u (4) )(x 1 ) + 1 6 δ 4 x u * 2 -(u (4) )(x 2 )
. First, we consider the terms evaluated at x 1 . Recall that (49)

δ 4 x u * 1 = 12 h 2 (δ x u * x ) 1 -δ 2 x 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) (δ x u * x ) 1 = u * x,2 2h = u ′′ (x 1 ) + h 2 6 u (4) (x 1 ) + r 1 , |r 1 | ≤ Ch 3 u (5) L ∞ , and (51)
δ 2 x u * 1 = u ′′ (x 1 ) + h 2 12 u (4) (x 1 ) + r 2 , |r 2 | ≤ Ch 3 u (5) L ∞ .
Inserting the estimates (50), (51) in Equation ( 49), we obtain

(52) δ 4 x u * 1 = u (4) (x 1 ) + r 3 , |r 3 | ≤ Ch u (5) L ∞ . Next, for x 2 we have (53) δ 4 x u * 2 = 12 h 2 (δ x u * x ) 2 -δ 2 x u * 2 .
Expanding on the term (δ x u * x ) 2 and using again (33), we have

(54) (δ x u * x ) 2 = u * x,3 -u * x,1 2h = u ′′ (x 2 ) + h 2 6 u (4) (x 2 ) + s 1 , |s 1 | ≤ Ch 3 u (5) L ∞ .
For the second term δ 2 x u * 2 we have, as in (51),

(55) δ 2 x u * 2 = u ′′ (x 2 ) + h 2 12 u (4) (x 2 ) + s 2 , |s 2 | ≤ Ch 3 u (5) L ∞ .
Inserting the estimates (54), (55) in Equation (53), we obtain, as in (52), (56)

δ 4 x u * 2 = u (4) (x 2 ) + s 3 , |s 3 | ≤ Ch u (5) L ∞ .
Combining the estimates for r 3 and s 3 and inserting them in (48), we obtain (57)

|σ x δ 4 x u * 1 -σ x (u (4) ) * 1 | ≤ Ch u (5) L ∞ , which proves (42). (iii) Let t i = δ 4 x u * i -(u ( 4
) ) * i be the truncation error for the fourth-order derivative approximation. We have (58)

σ x t = v,
where v ∈ l 2 h,0 satisfies the estimates established in the previous parts of the lemma (59)

|v 1 |, |v N -1 | ≤ Ch u (5) L ∞ , |v i | ≤ Ch 4 u (8) L ∞ , 2 ≤ i ≤ N -2.
The representative matrix of σ x restricted to l 2 h,0 is P/6 = I -T /6. The eigenvalues of P/6 are (60) 1 -

2 3 µ k = 1 - 2 3 sin 2 kπ 2N .
The matrix norm of its inverse is (61)

|(P/6) -1 | 2 = max k=1,...,N -1 1 1 -2 3 sin 2 kπ 2N ≤ 3.
From (58), ( 59) and (61), we obtain, (62)

|t| h ≤ C|v| h .
Finally, since

(63) |v| 2 h ≤ Ch(2h 2 + N -2 i=2 h 8 )( u (5) 2 L ∞ + u (8) 2 L ∞ ) ≤ Ch 3 ( u (5) 2 L ∞ + u (8) 2 L ∞ ),
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. 5.1. Matrix representation of the DBO. We have shown in [START_REF] Ben-Artzi | A fast direct solver for the biharmonic problem in a rectangular grid[END_REF] 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

(64) SU = 12 h 2 3 2h 2 KP -1 K + 1 h 2 T U = 6 h 4 3KP -1 K + 2T U.
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

(67) S = 6 h 4 P -1 T 2 + 36 h 4 V 1 V T 1 + V 2 V T 2 , where the vectors V 1 , V 2 are (68)            V 1 = (α -β) 1/2 P -1 √ 2 2 e 1 - √ 2 2 e N -1 V 2 = (α + β) 1/2 P -1 √ 2 2 e 1 + √ 2 2 e N -1 .
The constants α, β are

(69) α = 2(2 -e T 1 P -1 e 1 ) β = 2e T N -1 P -1 e 1 .
Remark 5. In view of the positivity of P -1 , we have 0 ≤ e T 1 P -1 e 1 ≤ 1/2 and |e T N -1 P -1 e 1 | ≤ 1/2, so that 3 ≤ α ≤ 4 and |β| ≤ 1. Thus, (α ± β) 1/2 are well defined.

5.2.

Error estimate for the one-dimensional Stephenson scheme. In [START_REF] Ben-Artzi | Convergence of a compact scheme for the pure streamfunction formulation of the unsteady Navier-Stokes system[END_REF] 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 [START_REF] Fishelov | Recent developments in the pure streamfunction formulation of the Navier-Stokes system[END_REF] 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 = uu * , The grid function u * satisfies (70)

δ 4 x u * i = f * (x i ) + r i , 1 ≤ i ≤ N -1,
where r is by definition the truncation error. We later refer to Proposition 3 for estimates on r.

The error e = uu * satisfies (71) δ 4

x e i = -r i , 1 ≤ i ≤ N -1, e 0 = 0, e N = 0, e x,0 = 0, e x,N = 0.

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 = uu * satisfies

(72) |e| h ≤ Ch 4 ,
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 = uu * and r, respectively. Using the matrix representation (67), we can write Equations ( 1) and ( 70 We therefore have

(75) SE = -R.
In view of (67) we have that (76)

P SP = 6 h 4 T 2 P + 36 h 4 JJ T , where (77) J = √ 2 2 [(α -β) 1/2 (e 1 -e N -1 ), (α + β) 1/2 (e 1 + e N -1 )].
Inverting P SP and multiplying by P R, we have

(78) -P -1 E = P -1 S -1 R = (P SP ) -1 P R.
Our goal is to bound the elements of P -1 E by Ch 4 . Note that by Proposition 3 we have (79)

|(P R) 1 |, |(P R) N -1 | ≤ Ch, |(P R) j | ≤ Ch 4 , 2 ≤ j ≤ N -2.
Thus, we need to estimate (P SP ) -1 P R. We decompose P SP as follows (80)

P SP = GH -1 ,
where (81)

G = I + 6JJ T P -1 T -2 , H = h 4 6 P -1 T -2 , so that (82) (P SP ) -1 = HG -1 .
Note that with L = (6/h 4 )H, Q = 6JJ T , we have

(83) 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Λ ′ 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 ( 84)

κ j = 6 -λ j = 6 -4 sin 2 ( jπ 2N ), 1 ≤ j ≤ N -1.
The diagonal matrix Λ ′ contains the eigenvalues of H, which can be written as

θ j = h 4 6 λ -2 j κ -1 j = h 4 96 1 sin 4 ( jπ 2N )(6 -4 sin 2 ( jπ 2N )) , j = 1, • • • , N -1.
The element H i,k of the matrix H is

H i,k = N -1 j=1 Z i,j θ j Z j,k . (85) H i,k = N -1 j=1 h 4 96 2 N sin( ijπ N ) sin( jkπ N ) sin 4 ( jπ 2N )(6 -4 sin 2 ( jπ 2N ))
.

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

(86) (HG -1 P R) i = N -1 k=1 H i,k (G -1 P R) k ,
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)

H i,1 = N -1 j=1 h 4 96 2 N sin( ijπ N ) 1 sin 4 ( jπ 2N )(6 -4 sin 2 ( jπ 2N )) sin( jπ N ).
Recall the elementary inequalities

(88) sin x ≥ 2 π x, 0 ≤ x ≤ π 2 , (89) | sin x| ≤ |x|, 2 ≤ 6 -4 sin 2 ( jπ 2N ) ≤ 6.
Noting that h = 1/N and using the estimate

| sin( ijπ N )| ≤ 1, we obtain (90) |H i,1 | = |H 1,i | ≤ C N -1 j=1 h 5 1 (jh) 4 (jh) ≤ Ch 2 , i = 2, ..., N -2.
Similarly, we have

(91) C 1 h 3 ≤ H 1,1 ≤ C N -1 j=1 h 5 1 (jh) 4 (jh) 2 ≤ C 2 h 3 .
This estimate holds equally for H N -1,N -1 . For the other corner elements of H we have

(92) |H 1,N -1 | = |H N -1,1 | ≤ C 2 h 3 .
For i, k = 2, ..., N -2 we have

(93) |H i,k | ≤ C N -1 j=1 h 5 1 (jh) 4 ≤ Ch.
Therefore, the orders of magnitude of the elements of H are bounded by ( 94)

       Ch 3 Ch 2 . . .
       . 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 . matrix Q is (N -1) × (N -1), but it has only four non-zero components at the corner positions, (95)

Q 1,1 = Q N -1,N -1 = 6α, Q 1,N -1 = Q N -1,1 = 6β.
Therefore, QL has only two non-zero rows -the first and the last. The first row is given by

(QL) 1,j = 6(αL 1,j + βL N -1,j ), j = 1, ..., N -1
and the last row is given by

(QL) N -1,j = 6(βL 1,j + αL N -1,j ), j = 1, ...N -1.
Thus, (96)

   G 1,1 = 1 + 6(αL 1,1 + βL N -1,1 ) =: a 1 , G 1,N -1 = 6(αL 1,N -1 + βL N -1,N -1 ) =: a N -1 , G 1,j = 6(αL 1,j + βL N -1,j ) =: b j , j = 2, ..., N -2 and (97)
   G N -1,1 = 6(βL 1,1 + αL N -1,1 ) = G 1,N -1 = a N -1 , G N -1,N -1 = 1 + 6(βL 1,N -1 + αL N -1,N -1 ) = G 1,1 = a 1 , G N -1,j = 6(βL 1,j + αL N -1,j ) = b N -j , j = 2, ..., N -2,
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 ≤ j ≤ N -2) follow from those of the elements of L. Namely,

|a 1 |, |a N -1 | ≤ C/h and |b j | ≤ C/h 2 for j = 2, ..., N -2.
In what follows we shall need lower bounds for a 1 and a 2 1a 2 N -1 . 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) L 1,1 > |L 1,N -1 | L 1,1 ∓ L 1,N -1 = h 4 N -1 j=1 j even or odd sin 2 jπ N sin 4 jπ 2N (6-4 sin 2 jπ 2N ) ≥ Ch N -1 j=1 j even or odd (jh) 2 (jh) 4 = C h .
In the above, take "j = even" for "-" and "j = odd" for "+".

Using (96) and the bounds 3 ≤ α ≤ 4, |β| ≤ 1 (Remark 5), we get in view of ( 98) and (91)

(99) |a 1 | ≥ 6(3L 1,1 -|L N -1,1 |) -1 ≥ 12L 1,1 -1 ≥ C h .
Next, we treat the difference a 2 1a 2 N -1 . Since we have the upper bound |a 2 1a 2 N -1 | ≤ C 1 /h 2 , we again need only a lower bound. We write the difference a 2 1a 2 N -1 as (100)

a 2 1 -a 2 N -1 = [1 + 6(α + β)(L 1,1 + L 1,N -1 )] • [1 + 6(α -β)(L 1,1 -L 1,N -1 )
] (using the symmetries of L). In view of (98) and α ≥ 3 , |β| ≤ 1, we obtain

(101) |a 2 1 -a 2 N -1 | ≥ 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, ..., 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

(102) G 1 =          1 0 0 . . . 0 0 aN-1 a1 0 1 0 . . . 0 0 0 0 0 1 . . . 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 . . . 0 1 0 aN-1 a1 0 0 . . . 0 0 1         
.

The same operations on the identity matrix yield the matrix (103)

I 1 =           1 a1 -b2 a1 -b3 a1 . . . -bN-3 a1 -bN-2 a1 0 0 1 -0 . . . 0 0 0 0 0 1 . . . 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 . . . 0 1 0 0 -bN-2 a1 -bN-3 a1 . . . -b3 a1 -b2 a1
1 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) G 2 =           1 0 0 . . . 0 0 aN-1 a1 0 1 0 . . . 0 0 0 0 0 1 . . . 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 . . . 0 1 0 0 0 0 . . . 0 0 a 2 1 -a 2 N -1 a 2 1           .
The corresponding matrix I 2 (obtained similarly from I 1 ) is (105)

I 2 =           1 a1 -b2 a1 -b3 a1
. . .

I 3 =           1 a1 -b2 a1 -b3 a1
. . .

-aN-1 a 2 1 -a 2 N -1 b2aN-1-a1bN-2 a 2 1 -a 2 N -1 b3aN-1-a1bN-3 a 2 1 -a 2 N -1 . . . bN-3aN-1-a1b3 a 2 1 -a 2 N -1 bN-2aN-1-a1b2 a 2 1 -a 2 N -1 a1 a 2 1 -a 2 N -1          
.

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)

G -1 =           a1 a 2 1 -a 2 N -1 aN-1bN-2-a1b2 a 2 1 -a 2 N -1 aN-1bN-3-a1b3 a 2 1 -a 2 N -1 . . . aN-1b3-a1bN-3 a 2 1 -a 2 N -1 aN-1b2-a1bN-2 a 2 1 -a 2 N -1 -aN-1 a 2 1 -a 2 N -1 0 1 0 . . . 0 0 0 0 0 1 . . . 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 . . . 0 1 0 -aN-1 a 2 1 -a 2 N -1 aN-1b2-a1bN-2 a 2 1 -a 2 N -1 aN-1b3-a1bN-3 a 2 1 -a 2 N -1 . . . aN-1bN-3-a1b3 a 2 1 -a 2 N -1 aN-1bN-2-a1b2 a 2 1 -a 2 N -1 a1 a 2 1 -a 2 N -1          
.

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)

|(G -1 ) 1,1 |, |(G -1 ) N -1,N -1 |, |(G -1 ) N -1,1 | |(G -1 ) 1,N -1 | ≤ Ch.
Similarly, and using also |b j | ≤ C/h 2 , we get

(110) |(G -1 ) 1,j |), |(G -1 ) N -1,j | ≤ C h , j = 2, ..., N -2.
Therefore, the elements of G -1 are bounded by (111) We can now bound the elements of G -1 P R using (109)-( 110) and ( 79). ( 112)

G -1 =          Ch C/h C/
|(G -1 P R) 1 | ≤ N -1 k=1 |(G -1 ) 1,k | • |(P R) k | = |(G -1 ) 1,1 | • |(P R) 1 | + N -2 k=2 |(G -1 ) 1,k | • |(P R) k |+ +|(G -1 ) 1,N -1 | • |(P R) N -1 | ≤ C 1 h • h + C 2 (N -3)(1/h) • h 4 ≤ Ch 2 .
Similarly, we have that

|(G -1 P R) N -1 | ≤ Ch 2 . For i = 2, ..., N -2 (113) |(G -1 P R) i | = |(P R) i | ≤ Ch 4 .
Finally we consider the product HG -1 P R (see (78), ( 82))

(114) -P -1 E = HG -1 P R.

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

|(HG -1 P R) i | ≤ N -1 k=1 |H i,k | • |(G -1 P R) k | = |H i,1 | • |(G -1 P R) 1 | + N -2 k=2 |H i,k | • |(G -1 P R) k | + |H i,N -1 | • |(G -1 P R) N -1 | ≤ C 1 h 2 h 2 + C 2 (N -3)hh 4 ≤ Ch 4 .
Therefore, ( 116)

|(P -1 E) i | = |(HG -1 P R) i | ≤ Ch 4 , 1 ≤ i ≤ 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) |E| = |P P -1 E| ≤ C|P -1 E| (116) ≤ C N -1 i=1 (h 4 ) 2 = Ch -1/2 h 4 . Thus, in view of the definition of the l 2 norm (118) |e| h ≤ Ch 4 .
This proves the fourth order error estimate result.

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 |, e 2 = u comp -u exact l 2 = |u comp -u exact | h .
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 on [a, b] = [0, 1]. It satisfies (119) u(0) = u ′ (0) = u(1) = u ′ (1) = 0.
Thus, choosing (120) f (x) = u (4) = 360x 2 -240x + 24, then u(x) is the unique solution of the biharmonic problem (121)

   u (4) = f, 0 < x < 1, u(0) = u(1) = 0, u ′ (0) = u ′ (1) = 0.
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 (122)    δ 4 x u j = f (x j ), 1 ≤ j ≤ N -1, u 0 = u N = 0, u x,0 = u x,N = 0.
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.

Mesh

   u (4) = f, 0 < x < 1, u(0) = 0, u(1) = 1, u ′ (0) = 0, u ′ (1) = 5,
where the function f (x) is (124) f (x) = u (4) = 120x.

Thus, we resolve numerically (125)

   δ 4
x u j = f (x j ), 1 ≤ j ≤ N -1, u 0 = 0, u N = 1, u x,0 = 0, u x,N = 5.

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

6.2. Case 2: Oscillating solutions. We consider a family of functions defined by (126) u ε (x) = p(x) sin 1/q ε (x) ,

where the polynomial functions p(x) and q ε (x) are given by (127) p(x) = 16x 2 (1x) 2 , q ε (x) = 1/ (x -1/2) 2 + ε , ε > 0.

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.

In Figure 1 we display the functions u ε (x) corresponding to (128) ε = 7.510 -2 , ε = 5.010 -2 , ε = 2.510 -2 .

As in Case 1, we consider the approximation x u j = f ε (x j ), 1 ≤ j ≤ N -1, u 0 = 0, u N = 0, u x,0 = 0, u x,N = 0, 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.

Proposition 4 .h 4 e 1 (e 1 +

 411 (i) The operator σ x δ 4 x has the matrix form KP -1 e 1 ) T + e N -1 (e N -1 -KP -1 e N -1 ) T , where(66) e 1 = (1, 0, • • • , 0) T , e N -1 = (0, • • • , 0, 1) T .

  SU * = F + R.

.

  Now we divide the last row of G 2 and I 2 by

Figure 1 .

 1 Figure 1. The oscillating function x → 16x 2 (1x) 2 sin(1/((x -0.5) 2 + ε) for ε = 7.5 10 -2 (left), ε = 5.0 10 -2 (center), ε = 2.5 10 -2 (right) .

6. 3 .

 3 Case 3: A nonlinear biharmonic equation. As a final example we consider the nonlinear problem -H(u) = f, 0 < x < 1, u(0) = 0, u(1) = 0, u ′ (0) = 0, u ′ (1) = 0,where H is assumed to be a k-lipschitz function.

Figure 2 .

 2 Figure 2. Convergence rates for the problem (129) with ε = 7.5 10 -2 (left), ε = 5.0 10 -2 (center), ε = 2.5 10 -2 (right). Circles correspond to the maximum norm and squares correspond to the l 2 norm.

  Ch 2 Ch 3 Ch 2 Ch . . . Ch Ch 2 Ch 2 Ch . . . Ch Ch 2 Ch 3 Ch 2 . . . Ch 2 Ch 3

	. . .	. . .	. . . . . .	. . .

Table 1 .

 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.

		N = 16	Rate	N = 32	Rate	N = 64	Rate N = 128
	e max 7.8231(-6) 4.00 4.8894(-7) 4.00 3.0589(-8) 4.00 1.9106(-9)
	e 2	5.6157(-6) 4.00 3.5099(-7) 4.00 2.1937(-8) 4.00 1.3739(-9)

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 , on [a, b] = [0, 1]. The function u(x) is the solution of the biharmonic problem (123)

  Table 2 assesses the fourth-order accuracy of the scheme in this case too. max 9.6857(-7) 3.99 6.1118(-8) 4.00 3.8200(-9) 3.98 2.4129(-10 ) e 2 7.0187(-7) 4.00 4.3873(-8) 4.00 2.7420(-9) 4.00 1.7062(-10 )

	Mesh	N = 16	Rate	N = 32	Rate	N = 64	Rate	N = 128
	e

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.

This paper is dedicated to Professor Saul Abarbanel on the occasion of his 80-th birthday. * 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.

Mesh

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

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

where λ min is the smallest eigenvalue of the problem (134)

A sufficient condition for (133) to hold is that

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.