Matania Ben-Artzi

J.-P Croisille

D Fishelov

TIME EVOLUTION OF DISCRETE FOURTH-ORDER ELLIPTIC OPERATORS

Keywords: 2010 Mathematics Subject Classification. Primary 65M06; Secondary 65N12, 65M15 fourth order parabolic equation, discrete elliptic calculus, discrete evolution, biharmonic operator, Hermitian derivative, optimal convergence, boundary values, Kuramoto-Sivashinsky equation. 1

The evolution equation

is considered. A discrete parabolic methodology is developed, based on a discrete elliptic (fourth-order) calculus.

The main ingredient of this calculus is a discrete biharmonic operator (DBO). In the general case, it is shown that the approximate solutions converge to the continuous one. An "almost optimal" convergence result (O(h 4-ε )) is established in the case of constant coefficients, in particular in the pure biharmonic case. Several numerical test cases are presented, that not only corroborate the theoretical accuracy result, but also demonstrate high-order accuracy of the method in nonlinear cases. The nonlinear equations include the well studied Kuramoto-Sivashinsky equation. Numerical solutions for this equation are shown to approximate remarkably well the exact solutions.

The numerical examples demonstrate the great improvement achieved by using the DBO instead of the standard (five-point) discrete bilaplacian.

INTRODUCTION

Let Ω = [0, 1] and consider the linear evolution equation associated with the general one-dimensional fourth-order elliptic operator where we assume that (1.3) A(x), B(x), f (x) are real continuous functions, and A(x) ∈ C 1 (Ω).

The equation is supplemented by the initial data u(x, 0) = u 0 (x), and homogeneous boundary conditions

(1.4) u(0, t) = ∂ ∂x u(0, t) = u(1, t) = ∂ ∂x u(1, t) = 0, t ≥ 0.
The initial function u 0 is assumed to be sufficiently smooth, and for simplicity we take u 0 ∈ C ∞ (Ω). To justify our interest in this equation, we can gratefully quote from the Introduction of [START_REF] Russell | MOVCOL4: A moving mesh code for fourth-order time-dependent partial differential equations[END_REF]: "In recent years, fourthorder problems have also become important in image processing and for modeling diffusion processes in physics and material sciences." This paper deals with high-order discrete approximations to (1.2), and is a sequel to the "elliptic" paper [START_REF] Ben-Artzi | Discrete fourth-order Sturm-Liouville problems[END_REF]. It builds on the premise of a discrete functional calculus. The approximation procedure consists of functions and operators that are defined on a fixed grid; every ingredient of the continuous theory (not only functions and operators but also fundamental properties such as coercivity) has its discrete counterpart. There is no reference to any piecewise continuous functions, as would be the case with finite element or discontinuous Galerkin methods. An important feature of the discrete operators is their compactness, namely, depending only on nearest neighbors. The convergence of the discrete solutions to the exact one is measured solely on the grid, as it is being refined. A complete convergence proof in the general linear case is presented. In the constant coefficient case it is shown to be "almost optimal", as is further discussed below. This discretization methodology enables us to go beyond the linear case and consider the approximation to various well-known nonlinear evolution equations.

Since it has already been extensively implemented for the 2 -D Navier-Stokes system (see further below), we have chosen here to test it on the Kuramoto-Sivashinsky equation. This is carried out in the last section of the paper.

There is vast literature pertaining to the discrete approximation of higher order (in the spatial coordinates) evolution equations. The studies are divided to "categories" reflecting either interest in specific physical application or accounting for a particular numerical method. It would take a very comprehensive survey, beyond the scope of this paper, to give a full account of the literature. We shall just mention six representative examples (all including extensive bibliographies).

• The papers [START_REF] Bertozzi | A biharmonic-modified forward time stepping method for fourth order nonlinear diffusion equations[END_REF][START_REF] Boscarini | High order semi-implicit schemes for time dependent partial differential equations[END_REF][START_REF] Kassam | Fourth-order time-stepping for stiff PDEs[END_REF] focus mainly on the issue of an accurate discrete time stepping. The first deals with the biharmonic modified forward time-stepping; the highest order biharmonic operator is treated implicitly.

The second deals with extensions of the IMEX approach to high order semi-implicit schemes, mostly to second-order (in space) reaction-diffusion equations. The third introduces the modified exponential timedifferencing (ETD) method; the linear (say, biharmonic) part is eliminated by an "integrating factor", so that the nonlinear term can be separately handled by a suitable Runge-Kutta method. • In [START_REF] Bertozzi | A biharmonic-modified forward time stepping method for fourth order nonlinear diffusion equations[END_REF] the physical models of interest are the evolution of thin films and phase fields. The underlying assumption is that the fourth-order (spatial) elliptic problem at each step is accurately handled. In particular, periodic boundary conditions are imposed, avoiding the issue of boundary layers. The same type of physical applications is considered in [2], which provides an overview of finite volume and finite element methods, in one or two space dimensions. • The adaptive mesh methodology in [START_REF] Russell | MOVCOL4: A moving mesh code for fourth-order time-dependent partial differential equations[END_REF] is applied to the Cahn-Hilliard equation modeling phase separation.

The spatial part is handled using collocation by (piecewise) seventh-order Hermite polynomial. • In the area of image processing various models involving fourth-order evolution equations are extensively studied. In [START_REF] Greer | H 1 solutions of a class of fourth order nonlinear equations for image processing[END_REF] a non-linear fourth-order equation, which is used for noise reduction and simplification of two dimensional images, was investigated for the existence and uniqueness of solutions of the problem.

In [START_REF] Fishelov | Image inpainting via fluid equations[END_REF] a Navier-Stokes system in streamfunction formulation with a non-linear viscosity term was proposed and tested for image impainting. In [START_REF] Lysaker | Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time[END_REF] fourth-order partial differential equations were used for noise removal of Magnetic Resonance Images. • In [32, Section 3] a local discontinuous Galerkin (LDG) method is developed for biharmonic type equations, primarily in one space dimension. Their method is based on introducing new unknowns (related to successive derivatives of the solution) and then approximating the resulting first-order system. In [START_REF] Cheng | A discontinuous Galerkin finite element method for time dependent partial differential equations with higher order derivatives[END_REF] this reduction was circumvented by using fluxes (at cell boundaries) associated with higher order derivatives. • Optimal error analysis was carried out on finite element methods in the following papers. In [START_REF] Danumjayau | Finite element methods for one dimensional fourth order semilinear differential equation[END_REF] the author approximates a fourth-order time dependent differential equation by a finite element method and establishes optimal error estimates for the semi-discrete case. In [START_REF] Yan | Local discontinuous Galerkin methods for partial differential equations with higher order derivatives[END_REF]Section 3], an optimal convergence rate for the LDG method mentioned above was achieved. In [START_REF] Gregory | Local Mesh Refinement with Finite Elements for Elliptic Problem[END_REF] the authors treat the biharmonic problem using a conformal bicubic Hermite polynomial. Using interpolation theory it may be shown that if the solution is smooth enough then the error is bounded by Ch 4 , where h is the mesh size.

It should be emphasized that all treatments of the time dependent problems rely on a careful high order approximation of the biharmonic operator. Recently various LDG methods were proposed. A main tool is the interior penalty procedure, that must be properly chosen [START_REF] Süli | hp-version interior penalty DGFEMs for the biharmonic equation[END_REF]. It is beyond the scope of this paper to compare the compact scheme used here for ∂ 4

x with any variant of the DG methodology. We refer to [START_REF] Ben-Artzi | Discrete fourth-order Sturm-Liouville problems[END_REF] for a further discussion on this topic.

One of the oldest instances of physical models in which the biharmonic operator plays a crucial role is the 2 -D Navier-Stokes system in pure streamfunction formulation [START_REF] Ben-Artzi | Navier-Stokes Equations in Planar Domains[END_REF]. Indeed, the present study grew out of earlier work on high order approximation of this system [START_REF] Ben-Artzi | A Pure-Compact Scheme for the streamfunction formulation of Navier-Stokes equations[END_REF][START_REF] Ben-Artzi | Convergence of a compact scheme for the pure streamfunction formulation of the unsteady Navier-Stokes system[END_REF][START_REF] Kupferman | A central-difference scheme for a pure stream function formulation of incompressible viscous flow[END_REF].

The one-dimensional scheme was already extended to a two dimensional one in [START_REF] Ben-Artzi | A High Order Compact Scheme for the Pure-Streamfunction Formulation of the Navier-Stokes Equations[END_REF]. A non-optimal analysis was carried out for a second order compact scheme for the full nonlinear Navier-Stokes system in [START_REF] Ben-Artzi | Convergence of a compact scheme for the pure streamfunction formulation of the unsteady Navier-Stokes system[END_REF]. It was shown that the convergence rate is at least 3/2. That is suboptimal compared to the actual convergence rate (order 2), that was verified in numerical results.

In [START_REF] Ben-Artzi | A fast direct solver for the biharmonic problem in rectangular grid[END_REF] we developed a fast direct solver in two dimensions for the biharmonic problem in a rectangular grid. It is based on FFT expansions and on the Sherman-Morrison formula.

Our convergence analysis here is different from standard Summation by Parts arguments. Instead, it is based on a detailed discrete functional calculus. It runs parallel to the classical PDE theory and allows us to establish in a fully discrete context the coercivity properties of the operators (see Section 3). In turn these properties enable us to prove convergence, not only for the solution but also for its derivatives. In fact, the optimal (4-th order) convergence of the discrete biharmonic operator ( see Claim 4.1 below) was established in [START_REF] Ben-Artzi | Navier-Stokes Equations in Planar Domains[END_REF]Chapter 10] by invoking matrix methods. We refer to [START_REF] Ben-Artzi | Discrete fourth-order Sturm-Liouville problems[END_REF] for a detailed discussion.

The paper is organized as follows. In Section 2 we recall the basic well-posedness facts related to Equation (1.2) and in Section 3 we review the construction of the discrete (elliptic) operators, the analogs of the second-order and fourth-order derivatives. In particular, the discrete biharmonic operator (DBO) is defined in (3.11). The construction of the discrete operators depends in a fundamental way on the concept of Hermitian derivative, which is defined in (3.9). The present paper aims at proving (in the general case) the convergence of the discrete approximation to the continuous solution and establishing (in the constant coefficient case) the significant (from both theoretical and practical aspects) "high-order" accuracy of the approximation. The crucial point in obtaining such estimates is that the discrete operators possess the fundamental elliptic properties (such as coercivity), in full analogy with the continuous elliptic operators.

In Section 4 we study the pure biharmonic evolution equation, and prove the "almost optimal" convergence rate (see below a more detailed explanation). We treat this case first because of its great importance as a model in many classical problems of mathematical physics, such as elasticity theory or the streamfunction formulation of the 2 -D Navier-Stokes system.

In Section 5 we introduce (Equation (5.1)) the discrete elliptic operator L A,B,h , which is the discrete analog to L A,B . The semi-discrete equation, analogous to (1.2) is given in (5.2):

∂ ∂t v(x j , t) = -L A,B,h v(x j , t) + f * (x j ), 1 ≤ j ≤ N -1, t ≥ 0.
The unknown function v(x j , t) is a grid function defined on the discrete spatial grid and on the continuous time variable. By the "convergence of the scheme" we mean the convergence (in an appropriate sense) of the functions (in time) v(x j , t), t ≥ 0 to the corresponding restrictions of the exact solution u(x j , t), t ≥ 0 as the number N of the grid points increases. The convergence of the discrete scheme to the exact solution is stated in Theorem 5.3. In Section 6, the general constant coefficient case (A(x) ≡ a, B(x) ≡ b) is studied. We show that it can be reduced to the pure biharmonic case. To measure the accuracy of the scheme, we estimate the difference between the grid function and the restriction of the exact solution in terms of a power O(h β ), where h = 1 N is the mesh size (we consider here only uniform meshes). The exponent β > 0 is the "order of accuracy". An estimate with β = 4 is an "optimal " convergence rate. Even though we believe (based on the numerical test cases) that our scheme yields the optimal convergence rate, it remains an open problem (even for the pure biharmonic case). We are able to show (Theorem 6.2) only "almost optimal" convergence rate, meaning that we can take β = 4 -ε for any ε > 0 (but with a coefficient that may blow up as ε → 0).

In Section 7 we present a collection of numerical examples that demonstrate the convergence properties claimed in the preceding sections. It is worthwhile pointing out the great improvement in accuracy achieved by using the DBO instead of the standard (five-point) bilaplacian. This is explicitly demonstrated in the fourth test case.

In fact, we go beyond the linear theory expounded here and present numerical simulations of several nonlinear evolution equations. In particular, we consider the well-studied Kuramoto-Sivashinsky equation in the seventh test case.

In the Appendix we give details of the Taylor expansions of the finite difference operators addressed in the paper. For these expansions we give the nominal order, namely disregarding boundary conditions.

THE LINEAR EVOLUTION EQUATION-BASIC FACTS

As is well-known, non-homogeneous boundary conditions are accommodated by a modification of the righthand side function f (x). We use the standard notation for Sobolev spaces; H p (Ω) is the space of functions with (distributional) derivatives that are square integrable up to order p, while H p 0 (Ω) is the closure of C ∞ 0 (0, 1) in H p .

Defining a function w(x, t) = e ct u(x, t), for some positive constant c > 0, we see that w satisfies the equation

∂ ∂t w(x, t) = -L A,B+c w(x, t).
Remark that the initial and boundary conditions of w are identical to those of u. By taking c > 0 sufficiently large, and without changing notation, we can therefore assume that L A,B is coercive in the sense that, Claim 2.1. [1, Theorem 7.5] There exists a constant C > 0 such that, for every function z(x) ∈ H 4 (Ω), satisfying the homogeneous boundary conditions

z(0) = z (0) = z(1) = z (1) = 0,
we have

(2.1) 1 0 L A,B z(x) • z(x)dx ≥ C 1 0 [z (x) 2 + z (x) 2 + z(x) 2 ]dx.
By a standard argument about parabolic equations it follows that Equation (1.2) has a unique solution u(x, t).

Assume that the initial function u 0 ∈ C ∞ (Ω). Note that in view of the homogeneous boundary conditions (1.4) all time derivatives of the solution vanish on the boundary. By differentiating the equation N ≥ 1 times with respect to t, and taking the scalar product with ∂ ∂t N u(x, t) we obtain (2.2)

1 2 d dt 1 0 ∂ ∂t N u(x, t) 2 dx = - 1 0 ∂ ∂x 2 ∂ ∂t N u(x, t) 2 dx + 1 0 A(x) ∂ ∂x ∂ ∂t N u(x, t) 2 dx - 1 0 B(x) ∂ ∂t N u(x, t) 2 dx.
The integrations-by-parts in the right-hand side are justified by the boundary conditions (1.4). Invoking the Sobolev inequality (with φ = ∂ ∂t N u(x, t))

(2.3) 1 0 |φ (x)| 2 dx ≤ ε 1 0 |φ (x)| 2 dx + C(ε) 1 0 |φ(x)| 2 dx,
and taking a sufficiently small ε > 0 we infer that

(2.4) d dt 1 0 ∂ ∂t N u(x, t) 2 dx ≤ C 1 0 ∂ ∂t N u(x, t) 2 dx,
where C > 0 depends only on A(x), B(x). Gronwall's lemma now yields (2.5)

1 0 ∂ ∂t N u(x, t) 2 dx ≤ e Ct 1 0 ∂ ∂x 4N u 0 2 dx, t ≥ 0.
The equation (1.2) enables us to replace temporal derivatives by spatial derivatives, to any order, so the Sobolev embedding theorem yields, with a constant C(N, T, u 0 ) > 0, depending only on u 0 , N, T,

(2.6) sup x∈[0,1] ∂ ∂t N u(x, t) ≤ C(N, T, u 0 ), t ∈ [0, T ].

THE BASIC SETUP for DISCRETE APPROXIMATION

We equip the interval Ω = [0, 1] with a uniform grid

x j = jh, 0 ≤ j ≤ N, h = 1 N .
The approximation is carried out by grid functions v defined on {x j , 0 ≤ j ≤ N } . The space of these grid functions is denoted by l 2 h . For their components we use either v j or v(x j ). For every smooth function G(x) we define its associated grid function

(3.1) G * j = G(x j ), 0 ≤ j ≤ N.
Clearly G * depends on the parameter h = 1 N , but we do not indicate this dependence in order to simplify the notation. The discrete l 2 h scalar product is defined by

(v, w) h = h N j=0 v j w j ,
and the corresponding norm is

(3.2) |v| 2 h = h N j=0 v 2 j .
For linear operators A :

l 2 h → l 2 h we use |A| h to denote the operator norm. The discrete sup-norm is (3.3) |v| ∞ = max 0≤j≤N {|v j |} .
The discrete homogeneous space of grid functions is defined by

l 2 h,0 = v ∈ l 2 h , v 0 = v N = 0 . Given v ∈ l 2
h,0 we introduce the basic (central) finite difference operators

(3.4) (δ x v) j = 1 2h (v j+1 -v j-1 ), 1 ≤ j ≤ N -1, (δ 2 x v) j = 1 h 2 (v j+1 -2v j + v j-1 ), 1 ≤ j ≤ N -1.
The cornerstone of our approach to finite difference operators is the introduction of the Hermitian derivative of v ∈ l 2 h,0 , that will replace δ x . It will serve not only in approximating (to fourth-order of accuracy) first-order derivatives, but also as a fundamental building block in the construction of finite difference approximations to higher-order derivatives. First, we introduce the "Simpson operator"

(3.5) (σ x v) j = 1 6 v j-1 + 2 3 v j + 1 6 v j+1 , 1 ≤ j ≤ N -1.

Note the operator relation (valid in

l 2 h,0 ) (3.6) σ x = I + h 2 6 δ 2 x ,
so that σ x is an "approximation to identity" in the following sense. Let

ψ ∈ C ∞ 0 (Ω), then (3.7) |(σ x -I)ψ * | ∞ ≤ Ch 2 ψ L ∞ (Ω) , which yields (3.8) |(σ x -I)ψ * | h ≤ Ch 2 ψ L ∞ (Ω) .
In the above estimates the constant C > 0 is independent of h, ψ. The Hermitian derivative v x is now defined by

(3.9) (σ x v x ) j = (δ x v) j , 1 ≤ j ≤ N -1.
Remark 3.1. In the definition (3.9), the values of (v x ) j , j = 0, N, need to be provided , in order to make sense of the left-hand side (for j = 1, N -1). If not otherwise specified, we shall henceforth assume that, in accordance with the boundary condition

(1.4), v x ∈ l 2 h,0 , namely (v x ) 0 = (v x ) N = 0.
In particular, the linear correspondence l 2 h,0 v → v x ∈ l 2 h,0 is well defined, but not onto, since δ x has a non-trivial kernel.

We next introduce a fourth-order replacement to the operator δ 2

x (see [START_REF] Fishelov | Recent advances in the study of a fourth-order compact scheme for the onedimensional biharmonic equation[END_REF]Equation (15)], [6, Equation (10.50

)(c)]), (3.10) ( δ 2 x v) j = 2(δ 2 x v) j -(δ x v x ) j , 1 ≤ j ≤ N -1. The biharmonic discrete operator (DBO) is given by (for v, v x ∈ l 2 h,0 ), (3.11) δ 4 x v = 12 h 2 [δ x v x -δ 2 x v].
Note that, in accordance with Remark 3.1 the operator δ 2

x is defined on grid functions v ∈ l 2 h,0 , such that also v x ∈ l 2 h,0 . The connection between the two difference operators for the second-order derivative is given by

(3.12) -δ x 2 = -δ 2 x + h 2 12 δ 4 x .
Remark 3.2. Clearly the operators δ x , δ 2 x , δ 4

x depend on h, but for notational simplicity this dependence is not explicitly indicated.

Remark 3.3. For clarity, the Taylor expansions of the operators used are collected in the Appendix in Section A.

The fact that the biharmonic discrete operator δ 4

x is positive (in particular symmetric) is proved in [6, Lemmas 10.9, 10.10]. Therefore its inverse δ 4

x -1 is also positive. Beyond just positivity, a fundamental tool (analogous to classical elliptic theory) is the coercivity property (with C > 0 independent of h) [6, Propositions 10.11,10.13],

(3.13) (δ 4 x z, z) h ≥ C(|z| 2 h + |δ 2 x z| 2 h + |δ x z x | 2 h
), valid for any grid function z ∈ l 2 h,0 such that also z x ∈ l 2 h,0 . Owing to (3.10) we can add also (3.14) (δ 4

x z, z) h ≥ C| δ 2 x z| 2 h . Notation for time-dependent grid functions. A time-dependent grid function is denoted as {v(x j , t)} N j=0 or {v j (t)} N j=0 . When there is no risk of confusion we simplify the notation , replacing {v(x j , t)} N j=0 , {u * (x j , t)} N j=0 by v(t), u * (t) respectively.

THE PURE BIHARMONIC GENERATOR

We first establish convergence estimates for the simple equation

(4.1) ∂ ∂t u = - ∂ ∂x 4 u, x ∈ Ω, t ≥ 0,
subject to initial data u(x, 0) = u 0 (x), and homogeneous boundary conditions

u(0, t) = ∂ ∂x u(0, t) = u(1, t) = ∂ ∂x u(1, t) = 0, t ≥ 0.
The initial function u 0 is assumed to be sufficiently smooth, and we take for simplicity u 0 ∈ C ∞ (Ω). The semidiscrete analog of (4.1) is

(4.2) ∂ ∂t v(x j , t) = -δ 4 x v(x j , t), 1 ≤ j ≤ N -1, t ≥ 0, subject to initial data v(x j , 0) = (u * 0 ) j = u 0 (x j ),
and homogeneous boundary conditions

v(x 0 , t) = v x (x 0 , t) = v(x N , t) = v x (x N , t) = 0, t ≥ 0.
Here δ 4

x is the discrete biharmonic operator (3.11). Clearly the grid function v(•, t) depends on h, but we refrain from indicating it explicitly, as it will always be clear from the context (compare Remark 3.2). Occasionally the simplified notation v(t) will be used instead of v(•, t), see the end of Section 3. We recall the basic optimal convergence fact [6, Theorem 10.19]:

Claim 4.1. Let f (x) be a smooth function in Ω. Let g(x) satisfy d dx 4 g(x) = f (x),
subject to homogeneous boundary conditions. Then

(4.3) |g * -(δ 4 x ) -1 f * | ∞ = O(h 4 ). 4.1.
ALMOST OPTIMAL ERROR ESTIMATE-THE PURE BIHARMONIC CASE. We note that Equation (4.1) can be rewritten as

(4.4) ∂ ∂t ∂ ∂x -4 u = -u, x ∈ Ω, t ≥ 0,
In view of Claim 4.1 (with -∂ ∂t u(•, t) playing the role of f and u(•, t) playing the role of g) we have, for any T > 0, with O(h 4 ) depending only on u 0 and T,

(4.5) ∂ ∂t δ 4 x -1 u * (x j , t) = -u * (x j , t) + O(h 4 ), 1 ≤ j ≤ N -1, 0 < t ≤ T.
Note that by definition

∂ ∂t u * (x j , t) = ∂ ∂t u * (x j , t), 1 ≤ j ≤ N -1, t ≥ 0.
Similarly, Equation (4.2) can be rewritten as

(4.6) ∂ ∂t δ 4 x -1 v(x j , t) = -v(x j , t), 1 ≤ j ≤ N -1, t ≥ 0.
The smoothness of the initial data implies the boundedness of all time derivatives of the solution, as stated in the following claim.

Lemma 4.2. For every integer N ≥ 0

(4.7) ∂ ∂t N v(•, t) h ≤ |(δ 4 x ) N u * 0 | h , t ≥ 0.
Proof. Taking the scalar product of (4.2) with v(•, t) we have, due to the positivity of δ 4

x , 1 2

∂ ∂t |v(•, t)| h = -(δ 4 x v(•, t), v(•, t)) h ≤ 0.
Thus |v(•, t)| h is a nonincreasing function of t and its value at t = 0 is |u * 0 | h . This proves the claim for N = 0. We then proceed to all N by repeated differentiation of (4.2) with respect to t, and noting that

∂ ∂t N v(t = 0) = (-δ 4 x ) N u * 0 .
The difference between the exact solution and the discrete one is given by the "error" grid function

(4.8) e(t) = u * (t) -v(t).
The following theorem gives an estimate for this error.

Theorem 4.3. For every ε > 0 and every T > 0 the error can be estimated by

(4.9) |e(t)| h ≤ Ch 4-ε , t ∈ [0, T ], h < h 0 ,
where C > 0 depends only on u 0 , T, ε.

Proof. Define the grid function (4.10)

w(t) = (δ 4 x ) -1 e(t) = (δ 4 x ) -1 u * (t) -(δ 4 x ) -1 v(t
w (t) + e(t) = O(h 4 ), t ∈ [0, T ],
with w(0) = 0. Taking the discrete scalar product of (4.11) with w and using the coercivity property (3.13), noting that e(t) = δ 4

x w(t), we get (4.12)

d dt |w(t)| 2 h + C|w(t)| 2 h ≤ |(O(h 4 ), w(t)) h | ≤ O(h 8 ) + C 2 |w(t)| 2 h . It follows from Gronwall's inequality that (4.13) |w(t)| h = O(h 4 ), t ∈ [0, T ].
Take an integer Q > 1. We define a finite difference operator

(4.14) S Q w(t) = (∆t) -1 M k=-L a k w(t + k∆t), so that (4.15) |S Q w(t) -w (t)| = O((∆t) Q ), t ∈ [L∆t, T -M ∆t].
Observe that we can find such a difference operator, to any order, since all time derivatives of w(t) are bounded in the | • | h norm in view of the bounds (2.6), Lemma 4.2 and the uniform boundedness of the operator (δ 4 x ) -1 for h < h 0 . Note also that M, L depend on Q. Plugging this into (4.11) we obtain (4. [START_REF] Fishelov | Semi-discrete time-dependent fourth-order problems on an interval: error estimate[END_REF])

S Q w(t) + e(t) = O(h 4 ) + O((∆t) Q ), t ∈ [L∆t, T -M ∆t],
and the estimate (4.13), applied to S Q w(t), yields (4.17)

e(t) = O(h 4 ) + O((∆t) Q ) + O( h 4 ∆t ), t ∈ [L∆t, T -M ∆t]. Take now ∆t = h 4 Q . The last estimate yields (4.18) e(t) = O h 4(1-1 Q ) , t ∈ [Lh 4 Q , T -M h 4 Q ].
Since Q can be taken arbitrarily large, the theorem is proved.

A DISCRETE VERSION OF THE GENERAL EVOLUTION EQUATION

Using the finite difference operators introduced in Section 3, and taking h = 1 N , we introduce the discrete operator analogous to L A,B that was defined in (1.1):

(5.1) [L A,B,h g] i = (δ 4 x g) i -A * i ( δ 2 x g) i -(A ) * i (g x ) i + B * i g i , 1 ≤ i ≤ N -1
, where f * , A * , (A ) * , B * are the grid functions corresponding, respectively, to f (x), A(x), A (x), B(x). In order to formulate a semi-discrete analog of (1.2), we introduce the unknown grid function v(x j , t), depending on the continuous time parameter. We do not indicate explicitly the obvious dependence on "h", since this will be clear from the context, as in Section 4. The semi-discrete equation is therefore

(5.2) ∂ ∂t v(x j , t) = -L A,B,h v(x j , t) + f * (x j ), 1 ≤ j ≤ N -1, t ≥ 0,
subject to initial data v(x j , 0) = (u * 0 ) j = u 0 (x j ), and homogeneous boundary conditions

v(x 0 , t) = v x (x 0 , t) = v(x N , t) = v x (x N , t) = 0, t ≥ 0.
Remark 5.1. We assume that all grid functions and their Hermitian derivatives are in l 2 h,0 . This amounts simply to extending the grid functions (whose relevant values are at the interior points {x i , 1 ≤ i ≤ N -1}) as zero at the endpoints x 0 , x N .

In the proof of the convergence theorem below, an important ingredient is the coercivity of the operator, in analogy with Claim 2.1. It is actually a generalization of (3.13) that played a fundamental role in the case of the pure biharmonic generator. Lemma 5.2. By adding to B(x) a sufficiently large constant K > 0, the following coercivity inequality holds.

There exist positive constants C, h 0 > 0 such that, for every grid function z ∈ l 2 h,0 , where also z x ∈ l 2 h,0 , we have

(5.3) (L A,B,h z, z) h ≥ C | δ 2 x z| 2 h + |z x | 2 h + |z| 2 h , 0 < h < h 0 .
Proof. We have

(L A,B,h z, z) h = (δ 4 x z, z) h -(A * δ 2 x z, z) h -((A ) * z x , z) h + (B * z, z) h ≥ (δ 4 x z, z) h -ε[| δ 2 x z| 2 h + |z x | 2 h ] - β ε |z| 2 h + (B * z, z) h
, where we have used the Cauchy-Schwarz inequality and set

β = |A * | 2 ∞ + |(A ) * | 2 ∞ .
In view of (3.13) and (3.14) we can take ε > 0 sufficiently small and then add a suitable constant K > 0 to B * so that (5.3) is satisfied.

The basic convergence result here is that "stability" implies "convergence" as follows.

Theorem 5.3. [General convergence] Let u(x, t) be the solution to (1.2) and let v(x j , t) be the solution to (5.2). Let e(t) = u * (t) -v(t) (compare (4.8)). Then, for every T > 0, we have the following convergence result:

(5.4) lim h→0 sup 0≤t≤T |e(t)| h = 0.
Proof. The coercivity property stated in Lemma 5.2 implies, exactly as in the discussion prior to Claim 2.1, that by adding a sufficiently large constant c > 0 to B(x), without changing notation, we can assume that L A,B,h satisfies (5.3). Let τ (h) = o(h) be the truncation error arising in the application of the discrete operator L A,B,h to the grid function u * (x j , t), namely

(5.5) ∂ ∂t u * (x j , t) = -L A,B,h u * (x j , t) + f * (x j ) + τ (h), 1 ≤ j ≤ N -1.
Observe that τ (h) is in fact also a function of t, such that for any ε > 0 and T > 0 there exists h 0 > 0 so that

(5.6) |τ (h)| h ≤ εh, 0 < h < h 0 , 0 ≤ t < T.
Subtracting (5.2) from (5.5) we obtain (using the simplified notation for time dependent grid functions)

(5.7) ∂ ∂t e(t) = -L A,B,h e(t) + τ (h).
Taking the scalar product (in l 2 h,0 ) of (5.7) with e(t) and using the coercivity inequality (5.3) yields

(5.8) 1 2 d dt |e(t)| 2 h ≤ - 1 2 C|e(t)| 2 h + 2 C |τ (h)| 2 h
, where the Cauchy-Schwarz unequality was used to estimate the scalar product (τ (h), e(t)) h . In view of (5.6) and Gronwall's inequality the proof is complete.

∂ ∂t u(x, t) = -L a,b u(x, t), x ∈ Ω, t ≥ 0.
The equation is supplemented with initial data u(x, 0) = u 0 (x), and homogeneous boundary conditions

u(0, t) = ∂ ∂x u(0, t) = u(1, t) = ∂ ∂x u(1, t) = 0, t ≥ 0.
The initial function u 0 is assumed to be smooth, namely u 0 ∈ C ∞ (Ω). As already observed above (see (2.1)) by taking the constant b sufficiently large the coercivity can be assumed :

(L a,b z, z) L 2 (Ω) ≥ C z 2 L 2 (Ω) + z 2 L 2 (Ω) + z 2 L 2 (Ω) .
Next we introduce the corresponding discrete operator (as in (5.1))

L a,b,h v(x j , t) = δ 4 x v(x j , t) -a δ 2 x v(x j , t) + bv(x j , t), 1 ≤ j ≤ N -1.
The coercivity in the general case (5.3) yields in this special case the following estimate (for b > 0 sufficiently large).

(6.2) (L a,b,h z, z) h ≥ C | δ 2 x z| 2 h + |z x | 2 h + |z| 2 h , 0 < h < h 0 ,
valid for every grid function z ∈ l 2 h,0 , such that also z x ∈ l 2 h,0 . Claim 2.1 and (6.2) imply in particular that the operators L a,b , L a,b,h are invertible and the norms |L -1 a,b,h | h , 0 < h < h 0 are uniformly bounded. Throughout the rest of this section we assume that b > 0 is sufficiently large, as stipulated in these claims. The semi-discrete analog of (6.1) is

(6.3) ∂ ∂t v(x j , t) = -L a,b,h v(x j , t), 1 ≤ j ≤ N -1, t ≥ 0, subject to initial data v(x j , 0) = (u * 0 ) j = u 0 (x j ), 1 ≤ j ≤ N -1. and homogeneous boundary conditions v(x 0 , t) = v x (x 0 , t) = v(x N , t) = v x (x N , t) = 0, t ≥ 0.
In analogy with (4.4) we can rewrite (1.2) as (6.4)

∂ ∂t L -1 a,b u(x, t) = L -1 a,b ∂ ∂t u(x, t) = -u(x, t), x ∈ Ω, t ≥ 0.
Similarly, the analog of (4.6) is

(6.5) ∂ ∂t L -1 a,b,h v(x j , t) = L -1 a,b,h ∂ ∂t v(x j , t) = -v(x j , t), 1 ≤ j ≤ N -1, t ≥ 0.
The crucial element in our treatment is the "optimal" estimate for the (constant coefficient) elliptic case, which extends the same one in the pure biharmonic case (Claim 4.1).

Claim 6.1. [8, Theorem 5.7] Let f (x) be a smooth function in Ω. Let g(x) satisfy

L a,b g(x) = f (x),
subject to homogeneous boundary conditions. Then

(6.6) |g * -L -1 a,b,h f * | ∞ = O(h 4 ).
We can now apply Claim 6.1 to Equation (6.4), where ∂ ∂t u(x, t) corresponds to f, while -u corresponds to g. We infer that (6.7)

L -1 a,b,h ∂ ∂t u * (x j , t) = -u * (x j , t) + O(h 4 ).
Defining the error grid function e(t) = u * (t) -v(t) (compare (4.8)) we obtain, by subtracting (6.5) from (6.7)

(6.8) ∂ ∂t L -1
a,b,h e(t) = -e(t) + O(h 4 ). The following theorem gives the "almost optimal" convergence estimate for the general equation and constitutes a generalization of Theorem 4.3. Theorem 6.2. Suppose that u 0 ∈ C ∞ (Ω). For every ε > 0 and every T > 0 the error e(t) can be estimated by (6.9)

|e(t)| h ≤ Ch 4-ε , t ∈ [0, T ], h < h 0 ,
where C > 0 depends only on u 0 , T, ε.

Proof. Define the grid function (6.10)

w(t) = L -1 a,b,h e(t) = L -1 a,b,h u * (t) -L -1
a,b,h v(t). From Equation (6.8) we get (6.11)

w (t) + e(t) = O(h 4 ), t ∈ [0, T ],
with w(0) = 0. This equation is identical to Equation (4.11), so that the proof proceeds verbatim as that of Theorem 4.3. Note that instead of using the coercivity property (3.13), we use here the coercivity property (6.2).

NUMERICAL EXAMPLES -LINEAR AND NONLINEAR

Theorem 5.3 dealt with a general convergence result for linear equations with variable coefficients while Theorem 6.2 established an almost fourth-order accuracy of the semi-discrete scheme (6.3), in the case of constant coefficients. In this section, we first corroborate the assertion of Theorem 6.2 by four numerical test cases with constant coefficients. Indeed, the numerical results yield a full fourth-order accuracy. We then consider a linear equation with variable coefficients in the fifth test case. The excellent approximation (to known analytical solutions) achieved by employing the same difference operators is indicative of the applicability of the methodology developed here to more general evolution equations. Then, with the sixth and seventh test cases, we go beyond the linear setting that has been the subject matter of the general convergence theory expounded hitherto and demonstrate the effective fourth-order accuracy on several nonlinear biharmonic equations. This includes numerical solutions of the well-known Kuramoto-Sivashinsky equation, modelling the evolution of flame fronts (as well as other physical phenomena).

Discrete time stepping: Notice that in the convergence theory of the preceding sections we have used a "semi-discrete" scheme; the time variable was continuous. Numerical implementation necessitates discrete time intervals. Semi-discrete equations considered hereafter have the general form (7.1)

∂ t v = -δ 4 x v + L(v) + H(v) + F (t) = R(v, t), where L is linear in v, v x , δ2 x v and H is nonlinear is v, v x . F (t) = [F 1 (t), . . . , F N -1 (t)
] is a time-dependent forcing function.

Two time stepping schemes are considered. The first one is the two stage Radau-IIA scheme. It is an IRK (Implicit Runge-Kutta) scheme, of third order accuracy, with A-and L-stability. Using this scheme allows to focus only on observing the accuracy in space. This is true in particular for linear problems. In the computation, due to the third order accuracy in time, the time step is selected according to a relation ∆t = Ch 4/3 . The matrix A of the Radau-IIA scheme is [12, chap. 3, pp. The stage time values are ∆t/3 and ∆t, respectively. The scheme is expressed as

(7.3)        ξ 1 = u n + ∆t a 11 R(ξ 1 ) + a 12 R(ξ 2 ) , ξ 2 = u n + ∆t a 21 R(ξ 1 ) + a 22 R(ξ 2 )
with u n+1 = ξ 2 (FSAL 1 property). The system (7.3) is solved by Newton iterations. The second time stepping scheme is the IMEX (IMplicit-EXplicit) scheme already used in [START_REF] Ben-Artzi | A Pure-Compact Scheme for the streamfunction formulation of Navier-Stokes equations[END_REF]. Here, it is used only for the Kuramoto-Sivashinski equation (test cases 7 hereafter). As a rule, we denote by u * ,n (see (3.1)) the grid values, at time level the error , using discrete l 2 h norm, for the x-derivative of the solution.

Note that v n x is the Hermitian derivative (with respect to x) of v n . However we compare to the exact derivative ∂u ∂x and not to the Hermitian derivative (u * ) n x . We shall also use the sup norm (see

(3.3)) |e n | ∞ = |v n -u * ,n | ∞ .
Recall that e n , as a grid function, depends on h. To display the convergence rate, we make it explicit by wrrting e n h . Then the convergence rate between the error e h and e h/2 is calculated by Log 2 (|e h |/|e h/2 |).

7.1. First test case: periodic biharmonic equation. We consider on the interval I = [0, 2π) the biharmonic equation

(7.6)      ∂ ∂t u = - ∂ ∂x 4 u, x ∈ I, t ≥ 0 u(x, 0) = u 0 (x) = sin(x), x ∈ I.
When periodic boundary conditions are applied, the solution is u(x, t) = sin(x) exp(-t). The semi-discrete scheme is

(7.7)            dv j (t) dt = -δ 4 x v j (t), j = 0, 1, ..., N -1, v N (t) = v 0 (t), t ≥ 0, v j (0) = (u * 0 ) j := u 0 (x j
), j = 0, ..., N. The discrete time-stepping is the Radau-IIA time scheme (7.3). The periodic boundary conditions make this test case easier to implement compared to the case of homogeneous boundary conditions. However, it enables us to focus on the spatial accuracy of the discrete biharmonic operator δ 4

x . In Table 1, the accuracy of the scheme is compared to the piecewise cubic discontinuous Galerkin method of [START_REF] Yan | Local discontinuous Galerkin methods for partial differential equations with higher order derivatives[END_REF] (without local postprocessing). Both schemes are fourth-order accurate in space. The error levels are of the same order of magnitude for the two cases. 7.2. Second test case: polynomial initial data. We consider the equation

(7.8)                ∂ ∂t u = - ∂ ∂x 4 u + a ∂ ∂x 2 u -bu + f (x, t), 0 < x < 1, t ≥ 0 u(0, t) = ∂ ∂x u(0, t) = 0, u(1, t) = ∂ ∂x u(1, t) = 0, t ≥ 0 u(x, 0) = u 0 (x), 0 ≤ x ≤ 1.
The semi-discrete analog to (7.8) is (compare (6.3)) (7.9)

           dv j (t) dt = -δ 4 x v j (t) + a δ2 x v j (t) -b v j (t) + f * j (t), j = 1, ..., N -1, t ≥ 0, v 0 (t) = 0, v N (t) = 0, v x,0 (t) = 0, v x,N (t) = 0, t ≥ 0 v j (0) = (u * 0 ) j := u 0 (x j
), j = 0, ..., N. We picked a = 1, b = -1 in (7.8) with exact solution u(x, t) = x 4 (1 -x) 4 e t . The initial condition is therefore u 0 (x) = x 4 (1 -x) 4 . The function f (x, t) is chosen so that (7.8) is satisfied. Note that in Claim 2.1 it was required that b > 0 in order to ensure coercivity, and in particular that zero is not an eigenvalue of L a,b,h . However, for any value of b ∈ R, if we know that zero is not an eigenvalue, then the discrete scheme has the "almost optimal" rate of convergence. Indeed, coercivity is restored (without changing the scheme) by multiplying the discrete equation by a suitable function e ct (see the beginning of the proof of Theorem 5.3). Several runs were performed with successively refined meshes. Due to third-order accuracy of the Radau-IIA scheme, we picked ∆t = Kh 4/3 , where the constant K is adjusted so that 8 iterations are performed for N = 8. The numerical results are displayed in Table 2 and Figure 1. 7.3. Third test case: highly oscillating solution. We consider again Equation (7.8), but now the exact solution is (7.10) u ε (x, t) = p(x) sin 1/q ε (x) sin(2πt).

The polynomial functions p(x) and q ε (x) are given by (7.11)

p(x) = 16x 2 (1 -x) 2 , q ε (x) = (x -1/2) 2 + ε, ε > 0.
The parameter ε has to be adjusted so that the function u ε (x, t) oscillates faster in the middle of the interval [0, 1].

In this example we took ε = 0.05. Numerical results using the Radau-IIA time scheme are displayed in Table 3 and Figure 2. The results, even for the very coarse grid with N = 32, show the remarkable accuracy of the scheme. 7.4. Fourth test case: linear equation with spectral damped solution. We use this example to assess the difference between the discrete operators δ 4

x , (δ 2 x ) 2 , both of which serve to approximate the fourth-order derivative. Let ϕ k (x) (resp. µ k ) be the k-th eigenfunction (resp. k-th eigenvalue) of the spectral problem (7.12) d dx

4 ϕ(x) = µϕ(x),
with homogeneous boundary conditions

(7.13) ϕ(0) = d dx ϕ(0) = 0, ϕ(1) = d dx ϕ(1) = 0.
The eigenfunctions are analytically known, see for example [START_REF] Davies | Spectral theory of differential operators[END_REF][START_REF] Fishelov | Fourth-order convergence of a compact scheme for the one-dimensional biharmonic equation[END_REF]. We consider the time dependent problem

(7.14)                          ∂ ∂t u(x, t) = - ∂ ∂x 4 u(x, t) + (µ k - 1 τ )u(x, t), 0 < x < 1, t ≥ 0, u(0, t) = 0, ∂ ∂x u(0, t) = 0, t ≥ 0 u(1, t) = 0, ∂ ∂x u(1, t) = 0, t ≥ 0 u(x, 0) = ϕ k (x), k ≥ 1, 0 ≤ x ≤ 1.
Here τ > 0 is a parameter. The exact solution is given by the damping in time of the initial function ϕ k (

, t) = ϕ k (x) exp(-t/τ ). x), (7.15) u(x
The aim here is to recover a similar behavior for the discrete solution. With the notation of Section 6 we have

(7.16) a = 0, b = -(µ k -1/τ ), f (x, t) = 0.
In the numerical example we took k = 10 and τ = 10 µ k , leading to b = -9 10 µ k . We use this example also to consider the difference between two semi-discrete schemes for (7.14). The first consists of (see (5.2)):

(7.17)                    dv j (t) dt = -δ 4 x v j (t) + (µ k - 1 τ )v j (t), j = 1, . . . , N -1, t ≥ 0, v 0 (t) = 0, v N (t) = 0, t ≥ 0 v x,0 (t) = 0, v x,N (t) = 0, t ≥ 0 v j (0) = ϕ * k,j := ϕ k (x j ).
The second scheme replaces the compact δ 4

x by the standard five-point discrete biharmonic given in [START_REF] Bjørstad | Fast numerical solution of the biharmonic Dirichlet problem on rectangles[END_REF]:

(7.18) (δ 2 x ) 2 v j = v j+2 -4v j+1 + 6v j -4v j-1 + v j-2 h 4 , j = 2, . . . , N -2.
This scheme needs to be supplemented by an appropriate meaning for (δ 2 x ) 2 v j at the near-boundary points j = 1, N -1, taking into account the Dirichlet boundary conditions as in (7.17). In [START_REF] Bjørstad | Fast numerical solution of the biharmonic Dirichlet problem on rectangles[END_REF] this is accomplished by a quadratic extrapolation based on the value at the near-boundary point; for j = 1 the imposed boundary values v 0 = v x,0 = 0 are used, along with the value v 1 , obtaining an extrapolated value at the "ghost" point j = -1. Using this extrapolation approach leads to a formula for (δ 2

x ) 2 v 1 given by A similar approximation is applied at j = N -1. The operator (δ 2 x ) 2 is second order at points 2 ≤ j ≤ N -2. The associated semi-discrete scheme is (7.20)

(7.19) (δ 2 x ) 2 v 1 = 7v 1 -4v 2 + v 3 h 4 . 0 0.2 0.4 0.6 0. 8 1 -2 -1.5 -1
                   dv j (t) dt = -(δ 2 x ) 2 v j (t) + (µ k - 1 τ )v j (t), j = 1, . . . , N -1, t > 0, v 0 (t) = 0, v N (t) = 0, t ≥ 0, v x,0 (t) = 0, v x,N (t) = 0, t ≥ 0, v j (0) = ϕ * k,j := ϕ k (x j
), j = 1, . . . , N -1, Fig. 3 displays the oscillating shape of the eigenfunction ϕ 10 . Note that ϕ 10 is antisymmetric with respect to the point x = 1 2 . The corresponding eigenvalue is µ 10 11, 410, 019.97. The parameter τ represents a damping effect in time. The selected value is τ = 10 µ10 . The grids are N = 8, 16, 32, 64 and N = 128. The final time is t f = 1/µ 10 . Fig. 4 represents the numerical convergence history for the schemes (7.17) and (7.20). In Figs. 5 and6 the results of the schemes (7.17) and (7.20) are compared at final time for N = 32 and N = 64. The superiority of the operator δ 4

x over the operator (δ 2 x ) 2 is clearly observed: with N = 32 the scheme (7.17) is more accurate than the scheme (7.20) with N = 64. 7.5. Fifth test case: linear equation with non constant A(x) and B(x). We consider the function

(7.21) u(x, t) = U (x) exp(ωt), U (x) = exp - 1 x(1 -x)
.

The graph of U (x) is displayed in Fig. 7. Note that U (0) = U (1) = U (0) = U (1) = 0. The derivatives of U (x) are given by (7.22)

U (n) (x) = p n (x) q n (x) U (x).
The polynomials p 1 (x) and q 1 (x) are (7.23)

p 1 (x) = 1 -2x, q 1 (x) = x 2 (1 -x) 2 .
For n ≥ 2, the following recursion relations hold (7.24) p n+1 = q 1 (p n q n -q n p n ) + q n p n p 1 , q n+1 = q 2 n q 1 . The polynomial q n is q n = q 2 n -1

1

. We consider the operator L A,B (1.1) applied to u(x, t), where A(x) and B(x) are selected below in (7.29) and (7.27) respectively. Recall that (7.25) L A,B u(x, t) = ∂ 4 x u(x, t) -∂ x (A(x)∂ x u(x, t)) + B(x)u(x, t). (a) Convergence rate for the scheme (7.17) at final time t f = 1/µ 10 after time-stepping. ) is approximated by (7.17). The biharmonic operator is approximated by the fourth order operator δ 4

x . In the left panel, the grid size is N = 32. In the right panel, the grid size is N = 64. In both cases, the approximate solution is remarkably close to the exact one. Using (7.24) we get (7.26) ) is approximated by (7.20). The biharmonic operator is approximated by the second order operator (δ 2 x ) 2 . In the left panel, the grid size is N = 32. In the right panel, the grid size is N = 64. In both cases, the approximate solution is far from the exact one. 1 x(1-x) ).

(∂ t + L A,B )u(x, t) = ωU (x) exp(ωt) + p 4 (x) q 4 (x) U (x) exp(ωt) -∂ x A(x) p 1 (x) q 1 (x) U (x) exp(ωt) + B(x)U (x) exp(ωt).
Given a function A(x), if we define the function B(x) to be (7.27)

B(x) = -ω - p 4 (x) q 4 (x) + A (x) p 1 (x) q 1 (x) + A(x) p 1 (x) q 1 (x) + A(x) p 1 (x) q 1 (x) 2 = -ω - p 4 (x) q 4 (x) + A (x) p 1 (x) q 1 (x) + A(x) p 2 (x) q 2 (x) .
we see readily that u(x, t) is the solution of the equation

(7.28) ∂ t u(x, t) = -L A,B u(x, t).
We take (7.29) A(x) = 1 + 0.1 sin(2πkx), and the corresponding function B(x) is determined by (7.27). We then consider Equation (7.28) with the initial condition u(x, 0) = U (x) and and the semi-discrete approximation (see (5.2)) (7.30) where the operator L A,B,h is given by (5.1). The approximate solution of (7.30) is compared to the analytical solution u(x, t) = U (x) exp(ωt). In the following we have chosen k = 10 and ω = 1. In Fig. 8 we display the solution of the semi-discrete equation at final time t = t f after time-stepping using the grids N = 32 and N = 64. This solution is compared with the exact solution u(x, t f ) = U (x) exp(ωt f ). We also record the convergence rate , adding a coarser N = 16 grid and a finer one N = 128. Note that there is a delicate interaction between the three finite difference operators δ 4

∂ t v(t) = -L A,B,h v(t),
x v, δ2

x v and the Hermitian derivative v x . 7.6. Sixth test case: a non linear example. A nonlinear modification of (7.8) is the evolution equation (7.31)

               ∂ ∂t u = - ∂ ∂x 4 u + H(u(x, t)) + f (x, t), 0 < x < 1, u(0, t) = ∂ ∂x u(0, t) = 0, u(1, t) = ∂ ∂x u(1, t) = 0, u(x, 0) = u 0 (x).
The function f (x, t) is selected so that the exact solution is u ε (x, t) with ε = 0.05 (see (7.10)). Here H(u) = 100 sin 2 u. The initial data is u 0 (x) = u ε (x, 0). The analogous semi-discrete equation is

(7.32) d dt v j (t) = -(δ 4 x v) j + H(v j (t)) + f * (x j , t), 1 ≤ j ≤ N -1.
The time discretization is performed with the Radau-IIA scheme. The numerical results are displayed in Table 4 and Figure 9.

                 ∂ ∂t u = - ∂ ∂x 4 u - ∂ ∂x 2 u -u ∂ ∂x u + f, 0 < x < 1, t > 0, u(0, t) = 0, ∂ ∂x u(0, t) = 0, u(1, t) = 0, ∂ ∂x u(1, t) = 0.
This equation has been independently derived in the context of several extended physical systems driven far from equilibrium by intrinsic instabilities, including instabilities of dissipative trapped ion modes in plasmas, instabilities in laminar flame fronts [START_REF] Sivashinsky | Nonlinear analysis of hydrodynamic instability in laminar flames I. Derivation of basic equations[END_REF], phase dynamics in reaction-diffusion systems [START_REF] Kuramoto | Persistent propagation of concentration waves in dissipative media far from thermal equilibrium[END_REF], and fluctuations in fluid films on inclines [START_REF] Sivashinsky | On irregular wavy flow of a liquid film down a vertical plane[END_REF]. The equation generically describes the dynamics near long-wave-length primary instabilities in the presence of appropriate (translational, parity and Galilean) symmetries. We present computations of three numerical examples, where the exact solutions are known and the error of the approximate solutions (and their derivatives) can be precisely measured. In the third example, there is no forcing function and the evolving solution develops steep slopes.

(1) First numerical example. We take u(x, t) as (7.34) u(x, t) = e -t x 2 (1 -x) 2 , 0 < x < 1, t > 0, so that Equation (7.33) is satisfied with f (x, t) taken accordingly. As in the previous test cases, we obtain the approximate solutions v j (t) of u(x j , t) from the knowledge of the discrete data f * j (t) = f (x j , t) on the grid 0 = x 0 < • • • < x j < x N = 1. The approximate solutions are computed at discrete time levels t n = n∆t, using a temporal discretization that is analogous to the second-order IMEX time-scheme that has been used in [START_REF] Ben-Artzi | A Pure-Compact Scheme for the streamfunction formulation of Navier-Stokes equations[END_REF][START_REF] Ben-Artzi | A High Order Compact Scheme for the Pure-Streamfunction Formulation of the Navier-Stokes Equations[END_REF].

(7.35)          v n+1/2 j -v n j ∆t/2 + 1 2 (δ 4 x v n+1/2 j + δ2 x v n+1/2 j ) = - 1 2 (δ 4 x v n j + δ2 x v n j ) -v n j v n x,j + f * ,n j , v n+1 j -v n j ∆t + 1 2 (δ 4 x v n+1 j + δ2 x v n+1 j ) = - 1 2 (δ 4 x v n j + δ2 x v n j ) + v n+1/2 j v n+1/2 x,j + f * ,n+1/2 j ,
In Table 5 we display numerical results for (7.35). We have picked ∆t = h 2 , where h = 1/N , and the final time is t f = 0.25. Observe that the scheme (7.35) achieves fourth-order accuracy for u and ∂u ∂x . In Figure 10 we display the exact solution u(x, t) (solid line) at t = 0.25 and the computed solution v (circles) for N = 32 (left) and N = 64 (center). The convergence rates are documented in the right panel. Notice the fourth-order rates both for u and ∂ ∂x u . Notice also that even with a coarse mesh N = 16 the match is excellent.

(2) Second numerical example. We consider again the Kuramoto-Sivashinsky equation (7.33), and we take the exact solution u(x, t) as (7.36) u(x, t) = e -t x 4 (1 -x) 4 , 0 < x < 1, t > 0, with a suitable f (x, t). The discrete solution v is again obtained by the scheme (7.35). In Table 6 we display numerical results. We have picked ∆t = h 2 , where h = 1/N , and the final time is t = 0.25. Observe in the right panel that fourth-order accuracy is achieved for u and ∂u ∂x . In Figure 11 we display the exact solution u(x, t) at (solid line) t = 0.25 and the computed solution v (circles)for N = 32 (left) and N = 64 (center) with no forcing term (i.e. f (x, t) = 0). Here c = -0.1, k = 0.5 11/19 and x 0 = -10. The discrete solution v is again obtained by the scheme (7.35). In Table 7 we display numerical results for (7.35). for the fourth-order scheme. We have picked ∆t = h 2 , where h = 1/N , and the final time is t = 1. Observe in the right panel that fourth-order accuracy is achieved for u and ∂u ∂x . In Figure 12 we display the exact solution u(x, t) at (solid line) t = 1 and the computed solution v (circles) for N = 121 (left) and N = 961 (center) at t = 1 . Notice that even with a coarse mesh N = 121 the match is very good. The right panel shows the fourth order convergence rate for u and u x . Appendix A. Taylor expansions

In this appendix we collect Taylor expansions of the finite difference operators introduced in the text. For u(x) a regular function, we call u * the restriction of u(x) to the grid. For the Hermitian derivative (u * ) x,j we have ). The small constants in the truncation errors observed in (A.1), (A.4) and (A.5) partly explain the accuracy observed in the numerical results.

2 u

 2 , t) = -L A,B u(x, t) + f (x) (x, t) + A (x) ∂ ∂x u(x, t) -B(x)u(x, t) + f (x), x ∈ Ω = [0, 1], t ≥ 0,

6 . 4 z -a d dx 2 z

 642 THE LINEAR EVOLUTION EQUATION WITH CONSTANT COEFFICIENTS In Equation (1.2) we now consider the case A(x) ≡ a, B(x) ≡ b, f ≡ 0, where a, b are real constants. Using the operator notation (compare (1.1)) L a,b z = d dx + bz, the equation takes the form (6.1)

Figure 1 .

 1 Figure 1. Exact and calculated solution for (7.8) with a = 1, b = -1 and u(x, t) = x 4 (1 -x) 4 e t at final time t f = 10. Left: N = 8, center: N = 16, right: the estimated convergence rate in the maximum (circles) and l 2 h (squares) norms is 4. Radau-IIA time scheme with ∆t = Kh 4/3 .There are 8 time interations with the grid N = * . The error level at final time is excellent, even on the coarse grid N = 8.

Figure 2 .

 2 Figure 2. Exact and calculated solution u(x, t) = p(x) sin(1/q ε (x)) sin(2πt) for Equation (7.8) ε = 0.05, at final time t f = 0.75. The number N of grid points is N = 32, N = 64, N = 128 and N = 256. The convergence rate is indicated at final time t f = 0.75. Left: N = 32, center: N = 256, right: the observed convergence rate is close to 4 in the max norm (circles) and the l 2 h norm (squares). Radau-IIA time scheme with ∆t = Kh 4/3 .There are 10 time interations with the grid N = 32. Notice the good error level, even on the coarse grid N = 32, with 10 time iterations.

Figure 3 .

 3 Figure 3. The eigenfunction ϕ 10 of the spectral problem (7.12) associated to the eigenvalue µ 10 .

  Convergence rate for the scheme (7.20) at final time t f = 1/µ 10 after time-stepping.

Figure 4 .δ 4 x and the grid N = 64 Figure 5 .

 4645 Figure 4. Convergence rates for the discrete approximations to the time dependent equation (7.14). (a): convergence rates for (7.17). (b): convergence rates for (7.20). For both subfigures: N = 8, N = 16, N = 32, N = 64 and N = 128. (a): fourth order convergence rates with excellent error values for δ 4x . (b): second order convergence rates with poor accuracy for (δ 2 x ) 2 .

Figure 6 .

 6 Figure 6. Approximate and exact solution at final time t f = 1/µ 10 . The (7.14) is approximated by(7.20). The biharmonic operator is approximated by the second order operator (δ 2x ) 2 . In the left panel, the grid size is N = 32. In the right panel, the grid size is N = 64. In both cases, the approximate solution is far from the exact one.

Figure 7 .

 7 Figure 7. The function U (x) = exp(-1x(1-x) ).

Figure 8 .

 8 Figure 8. Equation (7.28) approximated by (7.30) with A in (7.29) and B in (7.27). The initial time is t i = 0. Results are shown at final time t f = 3. Left panel: N = 32, the convergence in space is not achieved. Center panel: N = 64, the convergence in space is achieved. Right panel: convergence rate in the maximum (circles) and l 2 h (squares) norms) based on the values obtained for N = 16, N = 32, N = 64 and N = 128. The convergence rate in the lower part is close to 4.

7. 7 .

 7 Seventh test case : Kuramoto-Sivashinsky equation. Consider the Kuramoto-Sivashinsky equation [28, 25] (7.33)

Figure 9 .

 9 Figure 9. Nonlinear equation (7.31) with exact and calculated solution u ε (x, t) = p(x)1/q ε (x) sin(2πt), ε = 0.05, at final time t f = 0.75. The grids contains N = 32, N = 64, N = 128 and N = 256 points. Left panel: N = 32, center panel: N = 256, right panel: convergence rate in the maximum (circles) and l 2 h (squares) norms. The rate is 4. The time scheme is the Radau-IIA scheme with a time step ∆t = Kh 4/3 and 10 time iterations for N = 32.

Figure 10 .

 10 Figure 10. First KS numerical example (7.33): Exact solution (solid line) and computed solution (circles) for N = 32 (left), N = 64 (center). The convergence rate is displayed in the right panel for u (circles) and ∂u ∂x (squares).

Figure 11 .

 11 Figure 11. Second KS numerical example: Exact solution (solid line) and computed solution (circles) for N = 32 (left) and N = 64 (center) The convergence rate for the KS equation is documented in the right panel for u (circles) and u x (squares).

( 3 )

 3 Third numerical example. We consider again the Kuramoto-Sivashinsky equation(7.33). This time we take the exact solution u(x, t)[START_REF] Xu | Local discontinuous Galerkin methods for the Kuramoto-Sivashinsky equations and the Ito-type coupled KdV equations[END_REF] as (7.37) u(x, t) = c + (15/19) 11/19(-9 tanh(k(xctx 0 )) + 11 tanh 3 (k(xctx 0 )), x ∈ [-30, 30],

Figure 12 .

 12 Figure 12. Third KS numerical example: Exact solution (solid line) and computed solution (circles) for N = 121 (left) and N = 961 (center) The convergence rate for the KS equation is documented in the right panel for u (circles) and ∂u ∂x (squares).

  ).

	Subtracting Equation (4.6) from Equation (4.5) we obtain
	(4.11)

  , of the exact solution. The discrete solution at this time level is designated as v n . The computational errors are designated as follows.(7.4) |e n | h = |v nu * ,n | h the error (of the computed solution) in the discrete l 2 h norm, as in (3.2),

	(7.5)	|e n x | h = v n x -	∂u ∂x	* ,n	h

1

First Same As Last t = t n

  j + O(h6 ), and for the operator δ x (u * ) x : Since the coefficient in factor of h 4 is 1/360 in (A.2) and (A.3) there is a cancellation when evaluating the DBO operator δ 4x u * j = 12(δ x u * x,j -δ 2 x u * j )/h 2 . This gives

					2
	(A.1) x u) (A.2) (u * x ) j = (∂ x u) * j -h 4 180 (∂ 5 δ x (u * x ) j = (∂ 2 x u) * j + h 2 6 (∂ 4 x u) * j + h 4 360 (∂ 6 x u) * j -	h 6 15120	(∂ 8 x u) * j + O(h 8 ).
	On the other hand the centered operator δ 2 x u * j has the expansion
	(A.3) x u) (A.4) δ 2 x u * j = (∂ 2 x u) * j + h 2 12 (∂ 4 x u) * j + h 4 360 (∂ 6 x u) * j + h 6 20160 (∂ 8 δ 4 x u * j = (∂ 4 x u) * j -h 4 720 (∂ 8 x u) * j + O(h 6 ).
	Finally the fourth order accuracy of the modified operator δ2 x in (3.10) is seen in the expansion
	(A.5)	δ2 x u * j = (∂ 2 x u) * j +	h 4 360	(∂ 6 x u) * j + O(h 8

* * j + O(h 8 ).

Table 1. Biharmonic periodic problem (7.6) with exact solution u(x, t) = sin(x) exp(-t) on x ∈ [0, 2π) at final time is t f = ln(2). The first line indicates the errors obtained by the scheme (7.7). The second line indicates the errors obtained by a piecewise cubic discontinuous Galerkin method [START_REF] Yan | Local discontinuous Galerkin methods for partial differential equations with higher order derivatives[END_REF] (without postprocessing). The time scheme Radau-IIA is used for (7.7). The time step is ∆t = Kh 4/3 with 5 iterations for N = 10. There are 80 iterations for the grid N = 80. Note that the DG scheme has 4 times more unknowns ( 4 unknowns per cell) that the compact scheme (one unknown per point).

mesh N = 8

Rate N = 16 Rate N = 32 Rate N = 64 |e| ∞ 2.46(+0) 4.02 1.64(-1) 4.13 1.01(-2) 3.99 5.83(-4) Table 2. Equation(7.8) with a = 1, b = -1 discretized by (7.9). The exact solution is u(x, t) = x 4 (1 -x) 4 e t . The initial time t 0 = 0 and the final time is t f = 10. The time scheme is the Radau-IIA scheme. The time step is ∆t = Ch 4/3 (8 time iterations for the coarsest grid N = 8). There is 125 iterations on the finest grid N = 64.


Matthieu Garnier
email: matthieu.garnier@chu-rouen.fr

Adelya Curado

Paul Billoir

Virginie Barbay

Matthieu Demeyere

Jean-Nicolas Dacher

Imaging of Oxford/AstraZeneca® COVID-19 vaccine-induced immune thrombotic thrombocytopenia

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

and management of complications [START_REF] Hani | COVID-19 pneumonia: a review of typical CT findings and differential diagnosis[END_REF][START_REF] Ammar | Lung cavitation in COVID-19 pneumonia[END_REF][START_REF] Cellina | Acute pulmonary embolism in a patient with COVID-19 pneumonia[END_REF]. Recently, a large vaccination campaign has been set up and very rare thrombotic complications have been reported in patients who have been vaccinated. We report a patient who developed hypercoagulable state after receiving COVID-19 vaccine (Oxford-AstraZeneca ® ).

A 26-year-old woman was admitted for acute stroke 8 days after administration of COVID-19 vaccine (Oxford-AstraZeneca ® ). Right hemiplegia and aphasia (Score of 8 using the National Institute of Health Stroke Scale) occurred rapidly while the patient was hospitalized for persistent nausea and headache that had started shortly after vaccination. Initial angiography revealed a proximal left middle cerebral artery occlusion. Dual thrombo-aspiration using the direct aspiration first pass technique (ADAPT technique, i.e., without stent retriever) performed 3.5 hours after the onset of symptoms led to recanalization of left middle cerebral artery after a first pass (Fig. 1A). Final angiography showed a 2C recanalization Abbreviations: CT, Computed tomography; VITT, Vaccine-induced immune thrombotic thrombocytopenia. https://doi.org/10.1016/j.diii.2021.04.005 based on thrombolysis in cerebral infarction scale (Fig. 1B). Laboratory tests revealed major thrombocytopenia, hypofibrinogenemia and inflammation (C-reactive protein 158 mg/L). Disseminated intravascular coagulation or catastrophic antiphospholipid syndrome were initially considered. However, high titre anti-platelet factor 4 antibodies (Lifecodes PF4 IgG Assay, Immucor ® ) and positive result of functional platelet activation testing were consistent with vaccine-induced immune thrombotic thrombocytopenia [START_REF] Greinacher | Thrombotic thrombocytopenia after ChAdOx1 nCov-19 vaccination[END_REF][START_REF] Oldenburg | Diagnosis and management of vaccine-related thrombosis following AstraZeneca COVID-19 vaccination: guidance statement from the GTH[END_REF]. CT examination revealed segmental pulmonary embolism and portal thrombosis extending to the splenomesenteric trunk and ileal veins while the patient was eupneic and had no abdominal pain. On arterial phase, CT showed global arterialization of the liver parenchyma with central hypoattenuating areas, that became isoattenuating on portal venous phase. These findings consistent with transient hepatic attenuation differences due to hepatic blood flow alterations (Fig. 1C). In addition, mesenteric venous enlargement and intraluminal filling defects, bowel wall edema and decreased ileal wall enhancement suggested venous mesenteric ischemia (Fig. 1D). The patient was given corticosteroids, plasmatic exchange and anticoagulants. Magnetic resonance imaging of the brain 14 days later showed minor sylvian and border zone infarcts, with mild hemorrhagic changes. At day 15, only gripping difficulties and minor phasic troubles were remaining.

As observed in COVID-19 [START_REF] Berre | Concomitant acute aortic thrombosis and pulmonary embolism complicating COVID-19 pneumonia[END_REF], Oxford/AstraZeneca vaccine may induce hypercoagulable state that could lead to concomitant venous and arterial thromboses. In our patient, immunoassay and confirmatory functional assay for PF4 antibodies were positive and consistent with vaccine-induced immune thrombotic thrombocytopenia, a disease that resembles heparin-induced thrombocytopenia. Post-vaccination thrombosis has been recently described, but imaging findings have rarely been reported. Although very rare, occurrence of multiple thrombosis in the days following Covid-19 vaccine should urge extensive evaluation in order to approach the diagnosis. Multimodality imaging has a pivotal role in the diagnosis and treatment of such condition.

  Since the beginning of the COVID-19 pandemic, radiologists have been working at the forefront of patient care. Computed tomography (CT) is a pivotal imaging essential modality by contributing to diagnosis and differentials as well as for the follow-up

Fig. 1 .

 1 Fig. 1. 26-year-old woman was admitted for acute stroke 8 days after administration of COVID-19 vaccine. A: Digital subtracted angiography of the left internal carotid artery in anteroposterior projection shows occlusion of the M1 segment (arrow) of the left middle cerebral artery and absence of opacification of the Sylvian network. National Institute of Health Stroke Scale score was 8. B: Cerebral digital subtracted angiography after double thrombo-aspiration (ADAPT technique) shows 2C recanalization of left middle cerebral artery after a first pass based on thrombolysis in cerebral infarction scale. Immediate clinical follow up was favorable. Fifteen days later, National Institute of Health Stroke Scale score was 2. C: CT image of the liver in axial plane obtained during the arterial phase of enhancement shows global arterialization of the liver parenchyma with central and perihilar perfusion defects anomalies (arrow) due to portal trunk thrombosis. Central areas became isoattenuating on portal venous phase. D: CT image of the abdomen in the coronal plane obtained during the portal venous phase demonstrates portal (black arrow) and superior and inferior mesenteric vein (white arrow) thrombosis. Decreased ileal wall enhancement (arrowhead) consistent with venous mesenteric ischemia is also present.

Funding

This work did not receive any grant from funding agencies in the public, commercial, or not-for-profit sectors.

Human rights

The authors declare that the work described has been carried out in accordance with the Declaration of Helsinki of the World Medical Association revised in 2013 for experiments involving humans.

Informed consent and patient details

The authors declare that this report does not contain any personal information that could lead to the identification of the patients. The authors also confirm that the personal details of the patients have been removed.

Author contributions

All authors attest that they meet the current International Committee of Medical Journal Editors (ICMJE) criteria for Authorship.

Disclosure of interest

The authors declare that they have no competing interest. Matthieu Garnier a, * Adelya Curado a Paul Billoir b Virginie Barbay b Matthieu Demeyere a Jean-Nicolas Dacher a,c a Department of Radiology, CHU de Rouen, 76000

Rouen, France b Vascular Hemostasis Unit, CHU de Rouen, 76000

Rouen, France c UNIROUEN, Inserm U1096, 76000 Rouen, France


Mohammad Abobala

Keywords: Representations, 2021, hal-03265163 Refined Neutrosophic Matrices Representations ,

come

Main Discussion

Definition 1:

Let

be two vector spaces over , consider any three linear transformations . We define the corresponding full AH-refined linear transformation as follows:

. We denote it by

Theorem 2:

Let be any full AH-refined linear transformation, then is linear by classical meaning.

Proof:

we have: , , we have: , Now we compute :

, .

Definition 3:

Let be a matrix over . We say that is the matrix of iff .

Theorem 4:

Let be any full AH-refined linear transformation, then is the matrix of iff is the matrix of , is the matrix of , is the matrix of .

Proof:

Firstly, we suppose that is the matrix of , is the matrix of , is the matrix of , where are classical linear transformations.

Consider the refined neutrosophic matrices , we have for each : , , hence is the matrix of .

Conversely, we assume that is the matrix of , we must prove that is the matrix of , is the matrix of , is the matrix of .

We have , thus:

So that, This implies that is the matrix of us the matrix of , thus is the matrix of and is the matrix of , thus is the matrix of .

In order to prove that every refined neutrosophic matrix can be represented by a unique AH-refined linear transformation, we interduce the following algorithm to derive a basis of from any classical basis of .

Theorem 5: Let be a basis of over , then:

Is a basis of over .

Theorem 6: Let be a linear function, then must be a full AH-refined linear transformation.

Example 7:

Now, we give an example to clarify the structure of the refined neutrosophic basis of a refined neutrosophic vector space basing on Theorem 5.

It is well known that the is a vector space over the real field R. Consider the corresponding refined neutrosophic vector space .

The basis of is the set { }. According to Theorem 5, the related refined neutrosophic basis of is: , , .

Example 8:

Now, we illustrate an example to clarify the representation of refined linear transformations by refined neutrosophic matrices.

Let = be a refined linear transformation, we have to write as an AH-refined linear transformation.

By Theorem 6, we know that must be equal to , where .

Firstly, we have , and

This means that: so that hence .

This implies that is a zero map.

For the computing of h, we regard that: .

So that:

. By remarking that , we get:

, hence , This means that: , Thus , i.e. is a zero map.

By the previous discussion, we get that .

The matrix of , thus the refined neutrosophic matrix of is . Now, we check that M is refined neutrosophic matrix of by the following computing:

.

Example 9:

In example 8, we showed how a refined neutrosophic linear transformation can be turned into a refined neutrosophic matrix, now we do the converse. We illustrate an example to show how a refined neutrosophic matrix can be turned into a refined neutrosophic linear transformation.

Let be a refined neutrosophic matrix. M can be written in the following form:

. M can be represented by a unique refined neutrosophic linear transformation = .