M Ben-Artzi

J.-P Croisille

D Fishelov

A CARTESIAN COMPACT SCHEME FOR THE NAVIER-STOKES EQUATIONS IN STREAMFUNCTION FORMULATION IN IRREGULAR DOMAINS

Keywords: Biharmonic problem, embedded compact scheme, Navier-Stokes equations

we introduced an embedded Cartesian scheme for the biharmonic problem in two dimensions.

Here we extend this methodology to the 2D Navier-Stokes system. Hermite (or Birkhoff) interpolation is invoked in one and two dimensions to obtain finite difference operators. The consistency analysis of the discrete formulas for irregular grids is emphasized. Numerical results demonstrate remarkable accuracy for a series of test cases for flows in elliptical domains.

Introduction

In [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF] an embedded Cartesian finite difference scheme for the biharmonic problem in irregular domains was introduced. A regular Cartesian grid is laid out in the embedded domain. In most of the interior physical domain the regular Cartesian grid is kept unchanged. Thus, points which are far from the boundary are the center of regular elements containing the points' immediate neighbors. However, irregular elements are constructed around points which are close to the boundary. As a result, the grid is not necessarily aligned to the boundary, as in the case of body fitted grids.

In case the physical domain is not a rectangular, but rather an irregular domain, the scheme results in a Cartesian 9-point compact scheme at points "far" from the boundary. However, the Cartesian stencil is distorted around points which are close to the boundary. A crucial tool in our embedding approach consists of classifying grid points near the boundary. First, the Cartesian points are labeled as interior points if they are inside the physical domain, exterior points if they are outside the domain, and boundary points if they are on the boundary. This is called a topological coding. Points which are inside the domain are then divided into three subcategories. This is named a logical coding (see Fig. 1). A first category contains regular calculated points; their neighboring points reside on the Cartesian grid and are far enough from the boundary. A second category contains edge points; these are points which are too close to the boundary and therefore they are not included in the computational stencil, but are replaced by points on the boundary. A third category contains irregular calculated points which are close to the boundary and at least one of their neighboring points is too close to the boundary. Irregular calculated points are surrounded by a distorted stencil. Fig. 2 shows the stencil of an irregular calculated point M 0 , which includes the edge point M2 . The topological and logical coding are shown in Fig. 1 for an ellipse. We refer to [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF][START_REF] Ben-Artzi | An embedded compact scheme for biharmonic problems in irregular domains[END_REF] for more details.

Extending finite difference approximations to irregular domains is not a new idea. It originates from early works such as [START_REF] Shortley | The numerical solution of Laplace's equation[END_REF]. Important contributions in the 70's with emphasis on the capacitance matrix method include [START_REF] Buzbee | The direct solution of the biharmonic equation on rectangular regions and the Poisson equation on irregular regions[END_REF][START_REF] Pereyra | High order fast Laplace solvers for the Dirichlet problem on general regions[END_REF]. A significant renewal of interest in this approach is recently observed. For example, Poisson solvers are suggested in [START_REF] Jomaa | The Shortley-Weller embedded finite-difference method for the 3D Poisson equation with mixed boundary conditions[END_REF][START_REF] Gibou | A second-order-accurate symmetric discretization of the Poisson equation on irregular domains[END_REF]. For time dependent problems we refer to [START_REF] Abarbanel | Asymptotically stable fourth-order accurate schemes for the diffusion equation on complex shapes[END_REF][START_REF] Ditkowski | High-order embedded finite difference schemes for initial boundary value problems on time dependent irregular domains[END_REF]. Concerning the Navier-Stokes equations we refer to [START_REF] Li | A fast finite difference method for solving Navier-Stokes equations in irregular domains[END_REF][START_REF] Hou | Convergence of a finite difference scheme for the Navier-Stokes equations using vorticity boundary conditions stream function-only formulation[END_REF]. Approximating the Navier-Stokes equations in pure stream function formulation using a high order compact scheme was first introduced in [START_REF] Ben-Artzi | A pure-compact scheme for the streamfunction formulation of Navier-Stokes equations[END_REF] in the case of a square domain (see also [START_REF] Gupta | A new paradigm for solving Navier-Stokes equations: streamfunction-velocity formulation[END_REF]). Extending this kind of scheme to irregular geometries was considered in [START_REF] Pandit | On the use of compact streamfunction-velocity formulation of steady Navier-Stokes equations on geometries beyond rectangular[END_REF][START_REF] Kalita | An efficient transient Navier-Stokes solver on compact nonuniform space grids[END_REF] using a boundary fitted grid. Another approach for handling Navier-Stokes equation in irregular geometries was introduced in [START_REF] Fishelov | A new fourth-order compact sheme for the Navier-Stokes equations in irregular geometries[END_REF]. We refer to Section 4 for the relation of the scheme [START_REF] Fishelov | A new fourth-order compact sheme for the Navier-Stokes equations in irregular geometries[END_REF] to the present work.

Here the compact scheme is based on a particular interpolation polynomial of degree six, combined with a compact interpolation for the gradient. The main features of this approximation have been introduced in [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF]. As in the approach to finite differences using an interpolating polynomial presented in classic treatises [START_REF] Collatz | The Numerical Treatment of Differential Equations[END_REF][START_REF]Numerical Analysis[END_REF][START_REF] Kunz | Numerical Analysis[END_REF], the key problem is to find a suitable polynomial which interpolates the "data" on a stencil around any calculated The point M 2 is labeled as an edge point since it it too close to the boundary to be a calculated point. The point M 2 is defined as the point on the boundary which lies on the line connecting M 0 and M 2 . point x i,j . For all calculated points, this is a 9 point compact stencil. At regular points, the stencil is regular; only at irregular points (near the boundary) the stencil is distorted (see Fig. 2). Differentiating this polynomial provides finite difference operators at x i,j in terms of the data. When collecting the discrete equations for all points x i,j the data become "unknowns" to be computed as solution of the global linear system. Note that a full mathematical convergence analysis of our scheme is not available yet. However, elements of error bounds and convergence analysis are available in one dimension [START_REF] Ben-Artzi | Navier-Stokes equations in planar domains[END_REF][START_REF] Fishelov | Recent advances in the study of a fourth-order compact scheme for the onedimensional biharmonic equation[END_REF][START_REF] Ben-Artzi | Spline functions, the biharmonic operator and approximate eigenvalues[END_REF][START_REF] Ben-Artzi | Discrete fourth-order Sturm-Liouville problems[END_REF][START_REF] Ben-Artzi | Time evolution of discrete fourth-order elliptic operators[END_REF]. The interest in this approach is supported by the remarkable accuracy of the numerical results obtained so far [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF][START_REF] Fishelov | A new fourth-order compact sheme for the Navier-Stokes equations in irregular geometries[END_REF][START_REF] Ben-Artzi | An embedded compact scheme for biharmonic problems in irregular domains[END_REF].

M 1 (-h 1 , h 1 ) M 2 (0, h 2 ) M 3 (h 3 , h 3 ) M 4 (-h 4 , 0) M 5 (h 5 , 0) M 6 (-h 6 , -h 6 ) M 7 (0, -h 7 ) M 8 (h 8 , -h 8 ) M 0 (0, 0) M 2 M 3 M 5 M 8
The outline of the paper is as follows. In Section 2 we summarize the one-dimensional design of the compact scheme using a polynomial approach. In Section 3 the structure of the two dimensional polynomial, introduced in [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF] is outlined. Emphasis is given on the directional splitting between the two axis directions x and y and the two diagonal directions x + y and x -y. In Section 4, we display several facts concerning discrete operators on compact stencils, regular or not, which are used in our approximation to the Navier-Stokes equation. Finally, Section 5 displays numerical results obtained for the Navier-Stokes system in two dimensions for irregular domains.

Irregular Finite Difference operators in one dimension

In this section, we describe the basics of our finite difference scheme for the one dimensional case. Let us consider the interval I = [a, b] with the irregular grid

(2.1) a = x 0 < x 1 < ... < x N -1 < x N = b.
We denote by h j the spacing between neighboring points, thus

h j = x j -x j-1 , j = 1, ..., N and h = [h 1 , . . . , h N ]. x j x j-1 x j+1 u j-1 u j-1 u j u j u j+1 u j+1
Figure 3. One dimensional Hermite interpolation (2.4). The unique polynomial p(x) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + a 4 x 4 + a 5 x 5 with the six data u j-1 , u j-1 , u j , u j , u j+1 , u j+1 is given in (2.5).

Furthermore, and h are defined by (see [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF]),

     h = max j h j , (i) = min j h j+1 h j , h j h j+1 . (ii) (2.2)
The design of finite difference operators is based on a particular interpolation procedure of the data. For clarity we consider the following two steps in the analysis.

• Interpolation of the data u j , u j : We first consider the case where each interpolating point x j carries the two interpolating data u j and u j . The values u j and u j are supposed to be known data. The primary data u j stands for some approximation of u(x j ) and the value u j approximates the derivative u (x j ), j = 1, . . . , N -1. In approximation theory [START_REF] Powell | Approximation theory and methods[END_REF], these data are considered to be given, either as exact data or as measurements. To each point x j , j = 1, . . . , N -1, we associate a polynomial p(x)

(2.3) p(x) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + a 4 x 4 + a 5 x 5 .
Assuming the knowledge of the six data u j , u j , u j+1 , u j+1 , u j-1 , u j-1 (see Fig. 3), p(x) is the solution of (2.4)

p(x j ) = u j , p(x j-1 ) = u j-1 , p(x j+1 ) = u j+1 , p (x j ) = u j , p (x j-1 ) = u j-1 , p (x j+1 ) = u j+1 .
The coefficients a α in (2.3) are expressed in terms of u j-1 , u j , u j+1 , u j-1 , u j , u j+1 by (2.5)

a 0 = u j , a 1 = u j (a α ) j = A 1,α (h j , h j+1 )u j-1 + (-1) α A 1,α (h j+1 , h j )u j+1 + A 2,α (h j , h j+1 )u j + B 1,α (h j , h j+1 )u j-1 -(-1) α B 1,α (h j+1 , h j )u j+1 + B 2,α (h j , h j+1 )u j ,
The rational fractions

A 1,α (h 1 , h 2 ), A 2,α (h 1 , h 2 ), B 1,α (h 1 , h 2 ) and B 2,α (h 1 , h 2 ), 2α ≤ 5 are (2.6)          A 1,2 (h 1 , h 2 ) = h 2 2 (5 h 1 + 3 h 2 ) h 2 1 (h 1 + h 2 ) 3 , A 2,2 (h 1 , h 2 ) = - 3 h 2 1 -4 h 2 h 1 + 3 h 2 2 h 2 1 h 2 2 , B 1,2 (h 1 , h 2 ) = h 2 2 h 1 (h 1 + h 2 ) 2 , B 2,2 (h 1 , h 2 ) = -2 h 1 -h 2 h 2 h 1 , (2.7)
         A 1,3 (h 1 , h 2 ) = -2 h 2 h 2 h 1 -h 2 2 + 5 h 1 2 h 1 3 (h 1 + h 2 ) 3 , A 2,3 (h 1 , h 2 ) = 2 (h 1 -h 2 ) h 1 2 -3 h 2 h 1 + h 2 2 h 1 3 h 2 3 , B 1,3 (h 1 , h 2 ) = - h 2 (2 h 1 -h 2 ) h 1 2 (h 1 + h 2 ) 2 , B 2,3 (h 1 , h 2 ) = h 1 2 -4 h 2 h 1 + h 2 2 h 1 , (2.8)
         A 1,4 (h 1 , h 2 ) = -5 h 2 h 1 -4 h 2 2 + 5 h 1 2 h 1 3 (h 1 + h 2 ) 3 , A 2,4 (h 1 , h 2 ) = 4 h 1 2 -7 h 2 h 1 + 4 h 2 2 h 1 3 h 2 3 , B 1,4 (h 1 , h 2 ) = -2 h 2 + h 1 h 1 2 (h 1 + h 2 ) 2 , B 2,4 (h 1 , h 2 ) = 2 h 1 -h 2 h 1 2 h 2 2 .
(2.9)

       A 1,5 (h 1 , h 2 ) = 2 2 h 1 + h 2 h 1 3 (h 1 + h 2 ) 3 , A 2,5 (h 1 , h 2 ) = 2 h 1 -h 2 h 1 3 h 2 3 , B 1,5 (h 1 , h 2 ) = 1 h 1 2 (h 1 + h 2 ) 2 , B 2,5 (h 1 , h 2 ) = 1 h 1 2 h 2 2 . Denoting (2.10) ū = [u 1 , . . . , u N -1 ], [u 1 , . . . , u N -1 ] ,
the discrete operator (D α x ū) j d α dx α u(x j ) (also denoted ∂ α x u(x j )) at x j is related to a α by (2.11) (D α x ū) j = α! (a α ) j α = 2, 3, 4, 5. We should denote the coefficients a α by a α (ū, h) and the discrete derivative (D α x ū) j by D α x (ū, h) j , but instead we adopt the notation (2.11). The operator (D α x ū) j approximates ∂ α x u(x j ) in terms of the six data (u k , u k ) j-1≤k≤j+1 and of the two mesh spacing h j , h j+1 . The truncation error analysis for D α x ū is given in Section 6.

• Approximating the derivative u j by Hermitian derivatives u x,j The discrete operators above depend on the gridfunction ū in (2.10). In particular, the vector of the derivatives [u 1 , ..., u N -1 ] is assumed to be known. However, in our context below, the data u j is not available 1 . Thus, u j must be interpolated from the primal data u j . Here, this interpolation is obtained by mean of an interpolation polynomial called q(x). This polynomial is a priori different from p(x) in (2.3). We define the polynomial q(x) associated with x j as (2.12)

q(x) = b 0 + b 1 x + b 2 x 2 + b 3 x 3 + b 4 x 4 .
The gridfunction u x = [u x,1 , . . . , u x,N -1 ] is calculated implicitly as follows. Suppose that u x,j-1 and u x,j+1 are given. Then the vector b

= [b 0 , b 1 , b 2 , b 3 , b 4 ] T ∈ R 5
is the solution of the Hermite interpolation problem (2.13) q(x j-1 ) = u j-1 , q(x j ) = u j , q(x j+1 ) = u j+1 , q (x j-1 ) = u x,j-1 , q (x j+1 ) = u x,j+1 .

1 In some numerical methods, as the the finite element Discontinuous Galerkin approach, u is kept as a primal unknown.

The vector u x = [u x,j ] 1≤j≤N -1 is then defined from the five data u j-1 , u j , u j+1 , u x,j-1 , u x,j+1 by (2.14) u x,j = q (x j ).

Relations (2.14) for j = 1, ..., N -1 form a linear system with unknown u x , which is given for a general grid by (2.15) α 1,j u x,j + u x,j + α 2,j u x,j = β 1,j u j-1 + β 2,j u j + β 3,j u j+1 .

The five coefficients α 1,j , α 2,j , β 1,j , β 2,j and β 3,j are (2.16)

         α 1,j = h 2 j+1 (h j+1 + h j ) 2 , α 2,j = h 2 j (h j+1 + h j ) 2 , β 2,j = 2(h j+1 -h j ) h j h j+1 , β 1,j = - 2h 2 j+1 (2h j + h j+1 ) h j (h j + h j+1 ) 3 , β 3,j = 2h 2 j (2h j+1 + h j ) h j+1 (h j + h j+1 ) 3 .
Relation (2.15) is called the Hermitian closure for the derivative.

In [3, Lemma 4.2], it is shown that (2.17) max 1≤j≤N -1 |u (x j ) -u x,j | ≤ C h 4 max x∈[a,b] |∂ 5 x u(x)|,
where is given in (2.2) ii . In the case of a regular grid h = h j = h j+1 , j = 1, ..., N -1, (2.15) becomes the standard Hermitian relation given by (2.18)

   1 6 u x,j-1 + 2 3 u x,j + 1 6 u x,j+1 = u j+1 -u j-1 2h , 1 ≤ j ≤ N -1, u x,0 = u (a), u x,N = u (b).
Assuming periodicity, the truncation error on a uniform grid is given by 2 .

(

2.19) (u * ) x,j -(∂ x u) * j = - h 4 180 (∂ 5 x u) * j + O(h 6 ).
Replacing the exact derivatives u j (= ∂ x u j ) in (2.11) by the approximate values u x,j gives the finite difference operator denoted by D α x u:

(2.20) ( D α x u) j = α! (a α ) j ([u], [u x ]) α = 2, 3, 4, 5. The gridfunction u x = [u x,1 , . . . , u x,N -1 ] is not [u 1 , . . . , u N -1 ]
but it is rather its approximation by (2.15). In addition this approximation is non-local (see Section 6). This implies that the truncation error at x j depends on the full vector h and not only on h j , h j+1 as in (6.5). Consider the particular case of the two operators D 2

x and D 4

x , which satisfy

(2.21) ( D 2 x u) j ∂ 2 x u(x j ), ( D 4 x u) j ∂ 4
x u(x j ). We have Proposition 2.1. Let u(x) be a smooth function. Assuming that the grid (2.1) is such that

(2.22) c h ≤ h j ≤ h, c > 0,
The truncation errors for D 2 x and D 4

x

(2.23) t2 (u) = ( D 2 x u * ) -(∂ 2 x u) * , t4 (u) = ( D 4 x u * ) -(∂ 4 x u) * satisfy the estimates (2.24) | t2 (u)| ≤ Ch 3 , (i) | t4 (u)| ≤ Ch, (ii)
where C is a constant depending only on u(x).

2 Refer to Section 6 for more comments on truncation error analysis.

Proof. The truncation error t2 (u) is decomposed as

(2.25) t2 (u) = D 2 x u * -(∂ 2 x u) * = D 2 x u * -D 2 x u * (I) + D 2 x u * -(∂ 2 x u) * (II)
.

The term (I) is the truncation error t 2 (u) in (6.5) a . It satisfies the estimate

(2.26) |(II)| ≤ C(u)h 4 .
Invoking (2.5), (2.6), (2.17), (2.20) and assuming that c ≤ , it is easy to see that (II) satisfies

(2.27)

|(II)| ≤ 2B 1,2 (h j , h j+1 )|u x,j-1 -∂ x u j-1 | + 2B 1,2 (h j+1 , h j )|u x,j+1 -∂ x u j+1 | + 2B 2,2 (h j , h j+1 )|u x,j -∂ x u j | ≤ Ch 3 ,
where C is a constant depending only on u. A similar argument is used for (2.24) ii .

In the case of a regular grid, (i.e with h = h j , 1 ≤ j ≤ N -1), the operators D 2 x u and D 4 x u are (2.28)

   D 2 x u j = δ2 x u j = 2δ 2 x u j -δ x u x,j , D 4 x u j = δ 4 x u j = 12 h 2 δ x u x,j -δ 2 x u j .
The nominal truncation error (i.e. assuming periodicity) of D 2 x and of of D 4

x are

(2.29)

       ( D 2 x u * ) j -(∂ 2 x u) * j = h 4 360 (∂ 6 x u) * j + O(h 6 ), ( D 4 x u * ) j -(∂ 4 x u) * j = - h 4 720 (∂ 8 x u) * j + O(h 6 ).
In [START_REF] Ben-Artzi | Convergence of a compact scheme for the pure streamfunction formulation of the unsteady Navier-Stokes system[END_REF][START_REF] Ben-Artzi | Navier-Stokes equations in planar domains[END_REF], the operators δ2

x and δ 4 x were introduced using the polynomial q(x) only, and not the polynomial p(x). The reason is that when using u x,j instead of ∂ x u(x j ) in p(x), then p(x) and q(x) are identical. Here for clarity of the presentation we choose to introduce separately the two polynomials p(x) and q(x). Another Hermitian closure for the derivative, different from (2.14), could be selected to define the approximate values u x,j (see for example Equation (3.37) in [START_REF] Ben-Artzi | A high order compact scheme for the pure-streamfunction formulation of the Navier-Stokes Equations[END_REF]). In fact, in the two dimensional setting introduced in [3], the polynomials p and q are distinct, as will be seen below in Section 3. Remark 2.2. We refer to [START_REF] Kuntzmann | Méthodes numériques. Interpolation -Dérivées[END_REF][START_REF] Prenter | Splines and variational methods[END_REF][START_REF] Powell | Approximation theory and methods[END_REF] for a general presentation of Hermite interpolation.

Remark 2.3. The linear system (2.15) is distinct from the one employed to define the cubic splines derivatives, [2, (2.1.14) p.12]. However, the two systems are close. Both coincide with (2.18) in the case of a regular grid. See also [START_REF] Ben-Artzi | Spline functions, the biharmonic operator and approximate eigenvalues[END_REF].

Finite difference scheme on two dimensional irregular grids by interpolation

In two dimensions we consider a Cartesian grid in a rectangular domain embedding an irregular domain Ω. The definition of the finite difference operators on Cartesian grids can be obtained in two ways. The first is by expanding finite difference formulas in Taylor series (see [START_REF] Lele | Compact finite-difference schemes with spectral-like resolution[END_REF]). In the second approach a local interpolating polynomial is defined around each computational point which is the center point of an element (regular or irregular) inside the physical domain. Successive differentiation of this polynomial provides discrete differential operators. This is done in [START_REF] Stephenson | Single cell discretizations of order two and four for biharmonic problems[END_REF] for the Biharmonic operator. The latter has been extended in [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF] to define a high order approximation to the biharmonic operator on the stencil shown in Fig. 4. The polynomial P (x, y), defined in [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF], is the solution of a specific Birkhoff interpolation problem. 3 The polynomial P (x, y) ∈ P 19 = Span{l k , k = 1, . . . 19} where the Second, the 10 data corresponding to directional derivatives are marked with arrows.

M 1 (-h 1 , h 1 ) M 2 (0, h 2 ) M 3 (h 3 , h 3 ) M 4 (-h 4 , 0) M 0 (0, 0) M 5 (h 5 , 0) M 6 (-h 6 , -h 6 ) M 7 (0, -h 7 ) M 8 (
M 1 (-h, h) M 2 (0, h) M 3 (h, h) M 4 (-h, 0) M 0 (0, 0) M 5 (h, 0) M 6 (-h, -h) M 7 (0, -h) M 8 (h, -h)
                       l 1 (x, y) = 1, l 2 (x, y) = x, l 3 (x, y) = x 2 , l 4 (x, y) = x 3 , l 5 (x, y) = x 4 , l 6 (x, y) = x 5 , l 7 (x, y) = y, l 8 (x, y) = y 2 , l 9 (x, y) = y 3 , l 10 (x, y) = y 4 , l 11 (x, y) = y 5 , l 12 (x, y) = xy, l 13 (x, y) = xy(x + y), l 14 (x, y) = xy(x -y), l 15 (x, y) = xy(x + y) 2 , l 16 (x, y) = xy(x -y) 2 , l 17 (x, y) = xy(x + y) 3 , l 18 (x, y) = xy(x -y) 3 , l 19 (x, y) = x 2 y 2 (x 2 + y 2 ).
The polynomial P (x, y) is equivalently expressed as (3.2) P (x, y) = P (0, 0) + P 1 (x) + P 2 (y) + xyQ(x, y), where Consider ψ(x, y) a given function and ψ * i,j its restriction to the grid. The quintic polynomials P (0, 0) + P 1 (x) and P (0, 0)+P 2 (y) solve interpolation problems of the form (2.4). The polynomial P (0, 0)+P 1 (x) ∈ Span{1, x, x 2 , x 3 , x4 , x 5 } is given by the 6 horizontally aligned data (see Fig. 4):

(3.3)                  (a) P 1 (x) = a 2 x + a 3 x 2 + a 4 x 3 + a 5 x 4 + a 6 x 5 , (b)
(3.4)    Γ 4 (ψ) = ψ * (M 4 ), Γ 5 (ψ) = ψ * (M 0 ), Γ 6 (ψ) = ψ * (M 5 ), Γ 13 (ψ) = -∂ x ψ * (M 4 ), Γ 14 (ψ) = ∂ x ψ * (M 0 ), Γ 16 (ψ) = ∂ x ψ * (M 5 ).
Similarly, the polynomial P (0, 0) + P 2 (y) ∈ Span{1, y, y 2 , y 3 , y 4 , y 5 } is determined by the 6 vertically aligned data

(3.5)    Γ 8 (ψ) = ψ * (M 7 ), Γ 5 (ψ) = ψ * (M 0 ), Γ 2 (ψ) = ψ * (M 2 ), Γ 18 (ψ) = -∂ y ψ * (M 7 ), Γ 15 (ψ) = ∂ y ψ * (M 0 ), Γ 11 (ψ) = ∂ y ψ * (M 2 ).
Consider now the 8 diagonally aligned data. They are denoted as:

(3.6)          Γ 3 (ψ) = ψ * (M 3 ), Γ 7 (ψ) = ψ * (M 6 ), Γ 9 (ψ) = ψ * (M 8 ), Γ 1 (ψ) = ψ * (M 1 ), Γ 12 (ψ) = (∂ x + ∂ y )ψ * (M 3 ), Γ 17 (ψ) = (-∂ x -∂ y )ψ * (M 6 ), Γ 19 (ψ) = (∂ x -∂ y )ψ * (M 8 ), Γ 10 (ψ) = (-∂ x + ∂ y )ψ * (M 1 ).
Proposition 3.1. For any values α 1 , . . . , α 8 , there exists a unique polynomial Q(x, y) of the form (3.3) c satisfying the 8 relations

(3.7) Γ k (Q) = α k ,
where k ∈ I Q = {3, 7, 9, 1, 12, 17, 19, 10}

Proof. Since the problem is linear, it is sufficient to prove that α k = 0, k ∈ I Q , implies that Q(x, y) ≡ 0. This is easily obtained by analyzing Q(x, y) along the two diagonal directions x = y, x = -y: it is found that Q is a quartic polynomial along these two directions with 5 zero data. This implies Q(x, y) ≡ 0 by considering the coefficient of xy(x 2 + y 2 ). As a corollary, there exists a unique solution P (x, y) which is a solution of the 19 × 19 linear system

(3.8) Γ k (P ) = Γ k (ψ * ), 1 ≤ k ≤ 19.
where the 19 data Γ k (ψ * ) are defined in (3.4), (3.5) and (3.6). Denoting

(3.9) Γ(ψ * ) = [Γ k (ψ * )] 1≤k≤19 :
the polynomial P (x, y) is expressed as, (3.10)

P (x, y) = 19 k=1 a k Γ(ψ * ) l k (x, y).
The expression of P (x, y) in terms of Γ(ψ * ) ∈ R 19 and the lengths h ∈ R 8 is obtained by symbolic calculation 4 . The decomposition (3.2) is applied as follows. First, the polynomials P 1 (x) and P 2 (y)) are calculated. They are of the form (2.3), in the directions x and y respectively. The coefficients of P (x, 0) (resp. P (0, y)) are obtained by (2.5) in terms of h 4 , h 5 (resp. h 2 , h 7 ). Second, the polynomial Q(x, y) is calculated by solving a 8 × 8 system deduced from Prop. 3.1. Third, the 8 data Γ l (Q) in (3.6) are expressed in terms of the 19 data Γ k (P ), k = 1, . . . , 19. This is easily obtained using (3.2).

Discrete operators on compact stencils for the Navier-Stokes equations

Let Ω ⊂ R 2 be a regular bounded domain. The Navier-Stokes system in pure stream function form and with viscosity ν > 0 is (4.1)

               ∂ t ∆ψ(t, x) + ∇ ⊥ ψ • ∇(∆) ψ(t, x) -ν∆ 2 ψ(t, x) = f (t, x), ψ(t, x) = g 1 (t, x), ∂ψ ∂n (t, x) = g 2 (t, x), x = (x 1 , x 2 ) ∈ ∂Ω. with x = (x 1 , x 2 ) ∈ Ω, t > 0. The velocity v(t, x) is related to ψ by (4.2) v(t, x) = ∇ ⊥ ψ(t, x) = (-∂ 2 ψ, ∂ 1 ψ).
Three spatial operators must be approximated: the Biharmonic operator ∆ 2 ψ, the Laplacian operator ∆ψ and the convective term

C(ψ) = ∇ ⊥ ψ • ∇(∆) ψ.
The approximations, described hereafter, are all based on the polynomial P . They are called ∆ 2 h , ∆ h and C h respectively. The numerical scheme used hereafter is second order accurate. It is expressed as

(4.3)        ∆ h ψ n+1/2 -∆ h ψ n ∆t/2 = -C h (ψ n ) + ν 2 ∆ 2 h ψ n + ∆ 2 h ψ n+1/2 + (f (t n+1/4 , .)) * , ∆ h ψ n+1 -∆ h ψ n ∆t = -C h (ψ n+1/2 ) + ν 2 ∆ 2 h ψ n + ∆ 2 h ψ n+1 + (f (t n+1/2 , .)) * .
Refer to [START_REF] Ben-Artzi | A pure-compact scheme for the streamfunction formulation of Navier-Stokes equations[END_REF][START_REF] Ben-Artzi | Navier-Stokes equations in planar domains[END_REF] for the derivation of (4.3).

4.1. The Biharmonic operator. In this section, we briefly review various finite difference Biharmonic operators used in regular stencils [START_REF] Ben-Artzi | A pure-compact scheme for the streamfunction formulation of Navier-Stokes equations[END_REF][START_REF] Ben-Artzi | Navier-Stokes equations in planar domains[END_REF] or irregular stencils [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF][START_REF] Ben-Artzi | An embedded compact scheme for biharmonic problems in irregular domains[END_REF][START_REF] Fishelov | A new fourth-order compact sheme for the Navier-Stokes equations in irregular geometries[END_REF]. The 9 point stencil in regular form in shown in Fig. 5 and in irregular form in Fig. 4. The Biharmonic operator is ,y). The basic discrete Biharmonic is the one of Stephenson [START_REF] Stephenson | Single cell discretizations of order two and four for biharmonic problems[END_REF]. It operates on the 9 point regular stencil. It is expressed as

(4.4) ∆ 2 ψ(x, y) = ∂ 4 x ψ(x, y) + ∂ 4 y ψ(x, y) + 2∂ 2 x ∂ 2 y ψ(x
(4.5) ∆ 2 h ψ i,j = δ 4
x ψ i,j + δ 4 y ψ i,j + 2δ 2 x δ 2 y ψ i,j . The Hermitian derivative ψ x,i,j (resp. ψ y,i,j ) is involved in δ 4

x , (resp. δ 4 y ), (see Section 6). Both δ 4

x and δ 4 y are fourth order. The operator ∆ 2 h is second order with truncation error

(4.6) ∆ 2 h (ψ) * i,j -(∆ 2 ψ) * i,j = h 2 6 ∂ 4 x ∂ 2 y + ∂ 2 x ∂ 4 y )ψ * + O(h 4 ).
The second order accuracy in (4.6) is due to the approximation of the mixed term (∂

2 x ∂ 2 y ψ) * i,j by δ 2 x δ 2 y ψ * i,j . A modified version ∆2 h of ∆ 2 h was introduced in [7, 17]: (4.7) ∆2 h ψ i,j = ∆ 2 h ψ i,j - h 2 6 (δ 2 x δ 4 y + δ 2 y δ 4 x )ψ i,j .
The perturbation is O(h 2 ). It is designed to cancel the leading truncation term in (4.6). The operator ∆2 h ψ i,j is fourth order. The truncation error is (4.8)

∆2 h (ψ) * i,j -(∆ 2 ψ) * i,j = -h 4 1 720 (∂ 8 x ψ) * + (∂ 8 y ψ) * + 1 72 ∂ 4 x ∂ 4 y ψ * - 1 180 (∂ 2 x ∂ 6 y ψ) * + (∂ 6 x ∂ 2 y ψ) * ) + O(h 6 ).
k=1 a k l k (x, y) is associated to the point M 0 with global coordinates (x i , y j ) and local coordinates (0, 0). The approximate Biharmonic operator at M 0 is defined in terms of the data

Γ(ψ * ) = [Γ l (ψ * )] 1≤l≤19 via (4.9) ∆ 2 h ψ i,j = 19 k=1
a k (Γ(ψ * ))∆ 2 l k (0, 0) = 24 a 5 (Γ(ψ * )) + a 10 (Γ(ψ * )) + 16 a 15 (Γ(ψ * )) -a 16 (Γ(ψ * )) .

It has been proved in Theorem 3.11 of [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF] that if (u, ∂ x u, ∂ y u) are considered to be known exactly, then in analogy with (6.5) c , ∆ 2 h ψ i,j in (4.9) satisfies (4.10)

|∆ 2 h (ψ * ) i,j -(∆ 2 ψ) * i,j | ≤ Ch 2 .
In fact, this estimate is optimal, since on a regular one dimensional grid, the operator ∆ 2 h coincides with D 4 x u j in (2.11), which is second order accurate (see (6.5)(c)).

In Section 5 we approximate the first order derivatives (∂ x u, ∂ y u) via the Hermitian derivative (2.15). The latter are invoked in the approximation (4.9) of the Biharmonic operator. At regular points, the compact stencil is regular (see Fig. 5), i.e. h k = h, 1 ≤ k ≤ 8. It has been proved in [7, pp. 227-228] that the operator (4.9) coincides with (4.7), and therefore is fourth order accurate. At irregular points, which are near boundary, the accuracy drops to a lower value.

Remark 4.1. An alternative approximate Biharmonic operator may be derived from the identity [START_REF] Ditkowski | [END_REF] (4.11)

∆ 2 ψ = 2 3

∂ 4 ψ ∂x 4 + ∂ 4 ψ ∂y 4 + 2 3 ∂ 4 ψ ∂η 4 + ∂ 4 ψ ∂ξ 4 .
In (4.11) (η, ξ) are the diagonal variables defined by η = (x + y)/ √ 2, ξ = (y -x)/ √ 2. A discrete counterpart of (4.11) was introduced in [START_REF] Fishelov | A new fourth-order compact sheme for the Navier-Stokes equations in irregular geometries[END_REF] as

(4.12) ∆ 2 h u i,j = 2 3 (δ 4 x u i,j + δ 4 y u i,j ) + 2 3 (δ 4 η u i,j + δ 4 ξ u i,j
). If all the first-order derivatives with respect to x, y, ξ, η are approximated via the Hermitian derivative, then the scheme is fourth-order accurate. If the first-order derivatives with respect to x, y are computed via the Hermitian derivative and the first-order derivatives with respect to ξ, η are computed from the first-order derivatives with respect to x, y via the chain rule, then the scheme is second order accurate. In [START_REF] Fishelov | A new fourth-order compact sheme for the Navier-Stokes equations in irregular geometries[END_REF], the approximate Biharmonic operator (4.12) was used only at near boundary points, whereas (4.7) was used at interior points. This resulted in an observed fourth order accuracy of the numerical results. This gives

(4.14) ∆ h ψ i,j = 2 a 2 (Γ(ψ * )) + a 8 (Γ(ψ * )) .
The discrete Laplacian is easily proved to be fourth order. On a regular grid, ∆ h satisfies

(4.15) ∆ h = δ2 x + δ2 y .
In what follows, the convective term C(ψ) is approximated by C h (ψ) defined by

(4.16) C h (ψ) = -ψ y (∆ • ∂ x )P (0, 0) + ψ x (∆ • ∂ y )P (0, 0) .
In terms of the coefficients a k (Γ(ψ * )) this gives the formula: (4.17)

C h (ψ) = -ψ y 6a 4 (Γ(ψ * )) + 2 {a 13 (Γ(ψ * ))) -a 14 (Γ(ψ * )} + ψ x 6a 9 (Γ(ψ * )) + 2a 13 {Γ(ψ * )) + a 14 (Γ(ψ * )} ,
where a 4 , a 9 , a 13 , a 14 are given in (3.3). In the case of a regular grid, (4.17) is second order accurate. For this reason, the global order in space of the semi discrete scheme (4.3) is second order. This accuracy will be observed numerically below in Section 5. Remark 4.2. We have chosen to approximate the convective term by (4.16) mainly for simplicity. We refer to [START_REF] Ben-Artzi | A high order compact scheme for the pure-streamfunction formulation of the Navier-Stokes Equations[END_REF], formula (3.37) for a fourth order accurate convective term in the regular case. Extending the approach in [START_REF] Ben-Artzi | A high order compact scheme for the pure-streamfunction formulation of the Navier-Stokes Equations[END_REF] to the present embedding context will be the topic of a forthcoming work. See also [START_REF] Fishelov | A new fourth-order compact sheme for the Navier-Stokes equations in irregular geometries[END_REF] for an alternative approach.

Numerical results for Navier-Stokes equation in pure streamfunction

In order to assess the accuracy of our scheme for the system (4.1) three numerical test cases are presented. In section 5.1 we refer to Test cases 1 and 2, where the Navier-Stokes problem (4.1) is approximated in an ellipse. In Test cases 1 and 2 the exact solutions are chosen to be polynomials in space and exponential in time. We inspect the accuracy of the numerical solution by measuring the errors in ψ, ∂ x ψ and ∂ y ψ. Test case 2 was also considered in [START_REF] Li | A fast finite difference method for solving Navier-Stokes equations in irregular domains[END_REF], where a Cartesian embedding scheme for the Navier-Stokes equations was applied.

In section 5.2 we approximate the solution of the Navier-Stokes equation in a square. The physical domain (the square) is rotated with some angle θ and then embedded in a bigger square, which is aligned to the x, y grid. We analyze numerically the way the angle of rotation θ of the physical domain affects the accuracy of the computation. For this purpose we consider an analytic solution of the Navier-Stokes equations in the embedded square. We observe the behaviour of the numerical solution as a function of the rotation angle θ. From the numerical results one may observe that the angle of rotation has limited impact on the quality of the numerical results. x 2 0.5 2 + y 2 0.25 2 = 1. This ellipse is equivalently defined as the level set zero of the function

(5.2) ϕ(x, y) = 4x 2 + 16y 2 -1.
In Test case 1, the forcing function f (t, x, y) is chosen to be the right-hand side of Equation (4.1), for which the exact solution is (5.3) ψ e (t, x, y) = (x 4 + y 4 ) 2 exp(t).

The spatial part (x 4 + y 4 ) 2 is a polynomial of order 8. This polynomial is not included in the space P 19 in (3.1). The viscosity is chosen as ν = 1/1000. The following errors are reported in Table 2 at the final time T f = 0.5:

(5.4) Table 2. Case 1: Time dependent Navier-Stokes equation with exact solution ψ(t, x, y) = (x 4 + y 4 ) 2 exp(t). The domain is the ellipse x 2 /0.

   max i,j |(ψ e ) * i,j (T f ) -ψ i,j (T f )|, max i,j |(∂ x ψ e ) * i,j (T f ) -ψ x,i,j (T f )|, max i,j |(∂ y ψ e ) * i,j (T f ) -ψ y,i,j (T f )|.
5 2 + y 2 /0.25 2 = 1. The embedding square is [-1, 1] × [-1, 1]
, which corresponds to a coarser grid of the ellipse. The final time is T f = 0.5 and the viscosity is ν = 0.001. The asymptotic order of convergence is close to 2 for ψ and grad ψ. This is due to the convective term (4.16) which is approximated to order 2. Note that the errors are very small in magnitude.

a finer grid. The detailed embedding grid is shown in Fig. 6 for the two cases. In comparing the two cases, we would like to see whether the accuracy deteriorates in cases where the irregular domain is embedded in a larger computational square. According to Table 1, there are 921 calculated points in the finer grid (which corresponds to embedding the domain in [-0.6, 0.6] × [-0.6, 0.6]) and 313 in the coarser grid (which corresponds to embedding the domain in [-1, 1] × [-1, 1]). As expected, a better accuracy is observed with the finer grid. However, in both cases the errors are small and comparable in magnitude and the convergence rates are satisfactory (around second order and above). See Tables 2,3 and Fig. 7.

Remark 5.1. The size of the computational domain is usually adjusted to the physical domain as much as possible to reduce the computational effort. However, with our present coding, a direct solver is used to solve the global linear systems. The points exterior to the physical domain as well as the edge points are dummy points associated to the identity matrix in global linear systems. Therefore, the computational effort is in fact smaller with the large square [-1, 1] × [-1, 1] than with the small square [-0.6, 0.6] × [-0.6, 0.6].

5.1.2. Test case 2. We consider again the same ellipse (5.1) embedded in the square [-0.6, 0.6] × [-0.6, 0.6], which corresponds to the fine grid in Sec. 5.1.1. The exact solution is where r(x, y) = x 2 + 4y 2 , α > 0 and K are parameters to be picked. The corresponding velocity is

(5.5) ψ e (t, x, y) = K r(x, y) α -( 1 4 ) α cos(t), if ϕ(x, y) ≤ 1 0, otherwise , Grid N × N N ite e ∞ Rate ( e
(5.6) v 1 (t, x, y) = -∂ y ψ e (t, x, y), v 2 (t, x, y) = ∂ x ψ e (t, x, y).
The vorticity is (5.7) ω e (t, x, y) = ∆ψ e (t, x, y).

In [START_REF] Li | A fast finite difference method for solving Navier-Stokes equations in irregular domains[END_REF] the particular case α = 1, K = 1 is considered. The shape of the exact and approximated solution are shown in Fig. 8 (left panel). In Table 4 and in Fig. 9 (left panel) the numerical results are reported using our scheme. Note that in this case the spatial part of the function ψ e belongs to Span(x 4 , y 4 , x 2 y 2 ) ⊂ P 19 . So that there is no truncation error in the spatial direction x, y when using the scheme (4.3) and therefore there is no error in space. The errors originate from the temporal discretization. Therefore, fourth order accuracy is expected and observed in Table 4.

A more difficult case is obtained when we chose the exact solution of (5.5) with parameters α = 2 and K = 20. The shape of the exact and approximated solutions are shown in Fig. 8 (right panel). The error levels for ψ, ∂ x ψ and ∂ y ψ are reported in Table 5. The convergence slopes are displayed in Fig. 9. The observed accuracy is around 4 for ψ(t, x, y) grad ψ(t, x, y) for the case α = 1, K = 1 and around 2 for the case α = 2, K = 20. The asymptotic convergence rate is close to 4 for ψ and grad ψ. In this case, the truncation error in space vanishes since the exact solution belongs to the space P 19 . The error is only due to the time integration scheme.

Grid Second, a grid 30 × 30 is used. In each case, the maximum errors for ψ(T f , x, y) and ∂ x ψ(T f , x, y) are evaluated at final time T f = 0.25 using the scheme (4.3) Several positions of Ω are shown in Fig. 10. The maximum errors for ψ(T f , x i , y j ) and ∂ x ψ(T f , x i , y j ) with the grid 20 × 20 in function of the position angle are reported in Fig. 11. The error for ∂ y ψ(T f , x i , y j ) is similar to the one obtained for ∂ x ψ(T f , x i , y j ). In Fig. 11, one sees that the error level slightly fluctuates in function of the angle within the interval [10 -4.6 , 10 -3.4 ] for ψ and [10 -3.3 , 10 -2.2 ] for ∂ x ψ. In Fig. 12 the error levels are reported with a grid 30 × 30. As expected, the error levels are smaller with the grid 30 × 30 than with the 20 × 20 grid. The behavior of the error appears to be piecewise linear in all cases. The jumps in the errors are attributed to the change of labelling of near boundary points when the physical domain rotates. Despite these jumps in the errors, the magnitude of the error varies by a factor of 10 only and remains small. Overall, the errors behave satisfactorily.

Appendix: Truncation error analysis for Finite differences

In this section, we present a truncation error analysis for finite differences on a bounded interval. Suppose given a regular function u(x) on [a, b]. The operator u(x) → u * maps the function u(x) on the gridfunction obtained by evaluating u(x) at gridpoints. We write (6.1)

u * = [u(x 1 ), . . . , u(x N -1 )],
meaning that the evaluation holds at all scales h > 0. The notion of gridfunction is convenient to express the truncation error of a finite difference operator. For example the truncation error u → t(u) of the operator δ 2

x is the gridfunction obtained by (6.2)

t(u) = δ 2 x (u * ) -(∂ 2 x u) * . The components of t(u) are (6.3) t(u) j = (δ 2 x u * ) j -(∂ 2 x u) * j = h 2 12 (∂ 4 x u) * j + O(h 4 ).
In (6.3) u(x) is assumed to be periodic on [a, b]. This corresponds to a nominal order of 2. For any regular function u(x), we denote ū * = [u * , (∂ x u) * ] where u * , (∂ x u) * are the gridfunctions associated to u(x) and ∂ x u(x) respectively. Consider for example the discrete operator D α x in (2.11). The truncation error u → t α (u) is (6. The leading term of the truncation error of D α x is found5 in the case of a general irregular grid (2.1) for α = 1, 2, 3 and 4 as (6.5)

                     t 2 (u) j = (D 2 x ū * ) j -(∂ 2 x u) * j = - 1 360 h 2 j h 2 j+1 ∂ 6 x u j + O(h 5 ), (a) t 3 (u) j = (D 3 x ū * ) j -(∂ 3 x u) * j = 1 60 (h 2 j h j+1 -h j h 2 j+1 )∂ 6 x u j + O(h 4 ), (b) t 4 (u) j = (D 4 x ū * ) j -(∂ 4 x u) * j = - 1 30 (h 2 j + h 2 j+1 -4h j h j+1 )∂ 6 x u j + O(h 3 ), (c) t 5 (u) j = (D 5 x ū * ) j -(∂ 5 x u) * j = 1 3 (h j+1 -h j )∂ 6 x u j + O(h 2 ) (d).
The leading terms in the right-hand sides of (6.5) are associated solely with the three gridpoints {x j-1 , x j , x j+1 } and the two mesh spacing h j , h j+1 . In the case where the boundary conditions are taken in account, the truncation error will in general depend on the full set of values ∂

x u(x k ), 1 ≤ k ≤ N -1. In this case, only an estimate of the truncation error is available in general. Even in the case of a regular grid with h = h j , the leading term of the truncation error also depends on the full set of values ∂ [START_REF] Ben-Artzi | A high order compact scheme for the pure-streamfunction formulation of the Navier-Stokes Equations[END_REF] x u(x k ), 1 ≤ k ≤ N -1. As an example we have for the Hermitian derivative on a regular grid the following Proposition 6.1 (Pointwise truncation error of the Hermitian derivative). Let u → u x be the Hermitian gradient defined in (2.18) with boundary data u x,0 = 0, u x,N = 0. Then the truncation error is expressed as the multi-point α i j (∂ (5) x u) * (x j ) + O(h 6 ).

where (6.7) Using the fact that -δ 2 x Z k = λ k Z k and that (Z k ) 1≤k≤N -1 is an orthonormal basis of R N -1 gives the gridfunction component t(u) i is

α i j = N -1 k=1 Z i k Z j k 1 -λ k /
t(u) i = - h 4 180 N -1 k=1 (1 -λ k /6) -1 Z i k N -1 j=1 Z j k (∂ 5 x u)(x j ) + O(h 6 ) = - h 4 180 N -1 j=1 α i j (∂ 5 x u) * j + O(h 6 ),
where the coefficients α i j are given by (6.7). Therefore the leading term in the truncation error involves all the values (∂ 5

x u) * j as given in (6.6).

Conclusion

The numerical approach to the Navier-Stokes equation in pure streamfunction introduced in [START_REF] Ben-Artzi | A pure-compact scheme for the streamfunction formulation of Navier-Stokes equations[END_REF] is extended to irregular domains using embedded finite differences. The methodology presented uses high order interpolating polynomials. The numerical results presented are promising. The accuracy obtained so far is remarkable and is in line with the one observed in the purely Cartesian case.

The approach needs to be extended to more complex problems, with higher Reynolds numbers and more irregular geometries. In addition higher order time schemes with A-and L-stable must be considered. Bifurcation phenomena for flows in elliptical domains will be explored in the near future. Finally performing a stability analysis in the spirit of [START_REF] Carpenter | The stability of numerical boundary treatments for compact high-order finite difference schemes[END_REF][START_REF] Svärd | Review of summation-by-parts schemes for initial-boundary-value problems[END_REF][START_REF] Wang | Convergence of summation-by-parts finite difference methods for the wave equation[END_REF] is also a perspective.

Figure 1 .

 1 Figure 1. Ellipse embedded in a square supporting a Cartesian grid. Left panel: Topological coding of the points: the points are exterior (black small circles), boundary (triangles) or internal (clear small circles). Right panel: the logical coding of the different kinds of interior points are subdivided into Regular calculated points (small circles), Irregular calculated points (large circle), Edge points (open squares).

Figure 2 .

 2 Figure 2. The figure shows the stencil of the irregular calculated point M 0 as the dash line. It consists of the set {M 1 , M 2 , M 3 , M 5 , M 8 , M 7 , M 6 , M 4 }. The points M 1 , M 4 , M 6 and M 7 belong to the Cartesian grid. The points M 2 , M 3 , M 5 and M 8 belong to the boundary of the domain. The point M 2 is labeled as an edge point since it it too close to the boundary to be a calculated point. The point M 2 is defined as the point on the boundary which lies on the line connecting M 0 and M 2 .

Remark 2 . 1 .

 21 The truncation results in Prop. 2.1 are in agreement with Equations (4.27) and (4.35) of[START_REF] Fishelov | A new fourth-order compact sheme for the Navier-Stokes equations in irregular geometries[END_REF].

Figure 5 .

 5 Figure 5. Same as for Fig. 4, but in the case of a regular stencil with step size h. homogeneous polynomials l k (x) are

4. 2 .

 2 Discrete Laplacian and convective term. Consider first the Laplacian. The discrete Laplacian used in the sequel is the operator ∆ h defined by (4.13) ∆ h ψ i,j = 19 k=1 a k (Γ(ψ * ))∆l k (0, 0).

5. 1 .

 1 Navier-Stokes equation in an ellipse. 5.1.1. Test case 1. In Test case 1 and Test case 2, we consider the Navier-Stokes equation in an irregular domain Ω, which is defined by the ellipse (5.1)

Fig. 7

 7 Fig.7put later displays the least-square convergence rates corresponding to this test case. Two embedding squares are considered. In the first case the embedding square is [-1, 1]×[-1, 1] and in the second it is [-0.6, 0.6]×[-0.6, 0.6]. The number of points in each direction x, y are chosen to be the same for both cases. Therefore, the square [-1, 1] × [-1, 1] corresponds to a coarser grid of the ellipse, and the square [-0.6, 0.6] × [-0.6, 0.6] corresponds to

Figure 6 .

 6 Figure 6. Cases 1 and 2: The ellipse x 2 /0.5 2 + y 2 /0.25 2 = 1. Left panel: coarser grid with embedding square is [-1, 1]×[-1, 1]. Right panel: finer grid with embedding square is [-0.6, 0.6]× [-0.6, 0.6]. The three types of points inside the domain are: Regular calculated points (small circles), Irregular calculated points (large circle) and Edge points (open squares).

Table 3 .

 3 Case 1: Time dependent Navier-Stokes equation with exact solution ψ(t, x, y) = (x 4 + y 4 ) 2 exp(t). The domain is the ellipse x 2 /0.5 2 + y 2 /0.25 2 = 1. The embedding square [-0.6, 0.6] × [-0.6, 0.6], which corresponds to a fine grid of the ellipse. The final time is T f = 0.5 and the viscosity is ν = 0.001.

Figure 7 .

 7 Figure 7. Case 1: Navier Stokes equations in the ellipse x 2 /0.5 2 + y 2 /0.25 2 = 1 with exact solution ψ(t, x, y) = (x 4 + y 4 ) 2 exp(t). Regression line for the errors on ψ and ∂ y ψ. Left panel: coarse grid with embedding square [-1, 1] × [-1, 1]. The averaged rate of convergence for |(∂ x ψ) * i,j (T f ) -ψ x,i,j (T f )|) is 3.41 (not shown). Right panel: fine grid with embedding square [-0.6, 0.6]×[-0.6, 0.6]. The rate of convergence for |(∂ x ψ) * i,j (T f )-ψ x,i,j (T f )|) is 3.47 (not shown). The grids are N = 20, 40, 80, 160. The accuracy for ψ and grad ψ ranges between 2 and 4. In addition, the error levels are very good, due to the polynomial regularity of the exact solution.

Figure 8 .Table 4 .

 84 Figure 8. Case 2: Ellipse (5.1) embedded in the square [-0.6, 0.6] × [-0.6, 0.6]. The Cartesian embedding grid is 40 × 40 and there are 80 time iterations in (4.3). The exact and approximated solutions ψ(T f = 0.5, x, y) are shown. Left panel: parameters α = 1, K = 1 in (5.5). Right panel: parameters α = 2, K = 20 in (5.5).

Case 2 : 5 . 2 .

 252 Time dependent Navier-Stokes equation with ν = 0.001 in the ellipse x 2 /0.5 2 + y 2 /0.25 2 = 1 with exact solution (5.5) with the parameter α = 2 and K = 20. The final time is T f = 0.5. The embedding square is [-0.6, 0.6] × [-0.6, 0.6]. The asymptotic order of convergence is 2. The error levels for ψ, ψ x and ψ y are very good. Invariance under grid rotation: a numerical study. In this test case we consider the domain Ω = [-0.5, 0.5] × [-0.5, 0.5] embedded in the computational square [-1, 1] × [-1, 1]. The Navier Stokes equations (4.1) with ν = 1/1000 are solved in Ω with a forcing term f (t, x, y) such that the radial, polynomial in space solution is ψ e (t, x, y) = (x 2 +y 2 ) 4 e -t . The computation is repeated for a series of rotated positions of [-0.5, 0.5]×[-0.5, 0.5] in [-1, 1] × [-1, 1]. This permits to evaluate the effect of the position of the calculated points on the global accuracy. In particular, the labelling of near-boundary points changes when the domain changes of position.We let rotate Ω from the position θ = 0 to the position θ = π/4 with an angle step ∆θ = π/360 (0.5 deg). For each position θ k = k∆θ, k = 0, . . . 90 the same computation is reproduced. First, a grid 20 × 20 is considered.

Figure 9 .

 9 Figure 9. Case 2: Navier Stokes equations in the ellipse x 2 /0.5 2 + y 2 /0.25 2 = 1 with exact solution (5.5). The embedding square is [-0.6, 0.6] × [-0.6, 0.6]. Left panel: Parameters α = 1 and K = 1 in (5.5). The rate of convergence for |(∂ x ψ) * i,j (T f ) -ψ x,i,j (T f )|) is 3.94 (not shown). Right panel: Parameters α = 2 and K = 20 in (5.5). The rate of convergence for |(∂ x ψ) * i,j (T f )ψ x,i,j (T f )|) is 2.59 (not shown).

4 Figure 10 .

 410 Figure 10. Effect of the position of the domain in the embedding square on the accuracy. The square [-0.5, 0.5] embedded in the computational square [-1, 1] × [-1, 1] is shown at positions θ = 0, θ = π/16, θ = π/8, θ = 3π/16 and θ = π/4. The type of near boundary points (Regular calculated, Irregular calculated or Edge) varies according to their position with respect to the boundary. See Fig. 1.

Figure 11 .

 11 Figure 11. Case 3: Maximum error for the Navier-Stokes equation in the square with size [-0.5, 0.5] × [-0.5, 0.5] embedded in [-1, 1] × [-1, 1]. The grid is 20 × 20. The exact solution is ψ(t, x, y) = (x 2 + y 2 ) 4 exp(-t). The final time T f = 0.25 is reached in 20 iterations. The same computation is repeated 91 times on the domain [-0.5, 0.5] × [-0.5, 0.5] rotated with an angle θ = k π 360 for k = 0, ..., 90. Left panel: maximum error for ψ(t, x, y) at time T f = 0.25 in function of the position of the domain. Right panel: maximum error for ψ x (t, x, y) at time T f in function of the position of the domain. The accuracy for ψ y is similar to the one for ψ x .

Figure 12 .

 12 Figure 12. Case 3: Maximum error for the Navier-Stokes equation in the square [-0.5, 0.5] × [-0.5, 0.5] embedded in [-1, 1] × [-1, 1]. The grid is 30 × 30. The exact solution is ψ(t, x, y) = (x 2 + y 2 ) 4 exp(-t). The final time T f = 0.25 is reached in 45 iterations. The same computation is repeated 91 times on the domain [-0.5, 0.5] × [-0.5, 0.5] rotated with an angle θ = k π 360 for k = 0, ..., 90. Left panel: maximum error for ψ(t, x, y) at time T f = 0.25 in function of the position of the domain. Right panel: maximum error for ψ x (t, x, y) at time T f in function of the position of the domain. The accuracy for ψ y is similar to the one for ψ x .

  h 8 , -h 8 ) Figure 4. Case of an irregular stencil: the 8 neighbor nodes of M 0 are the points M 1 , M 2 , M 3 , M 4 , M 5 , M 6 , M 7 and M 8 with h = [h 1 , h 2 , h 3 , h 4 , h 5 , h 6 , h 7 , h

[START_REF] Ben-Artzi | Time evolution of discrete fourth-order elliptic operators[END_REF]

]. The 19 data determining uniquely the polynomial P (x, y) in (3.2) are of two kinds. First, the 9 data for u are marked with circles.

  P 2 (y) = a 7 y + a 8 y 2 + a 9 y 3 + a 10 y 4 + a 11 y 5 ,

	(c)

Q(x, y) = a 12 + a 13 (x + y) + a 14 (x -y) + a 15 (x + y) 2 + a 16 (x -y) 2 + a 17 (x + y) 3 + a 18 (x -y) 3 + a 19 xy(x 2 + y 2 ).

Table 1 .

 1 Cases

	grid 61 × 61	Embed. sq. [-1, 1] × [-1, 1] Embed. sq. [-0.6, 0.6] × [-0.6, 0.6]
	Total number of points	61 × 61 = 3721	61 × 61 = 3721
	Number of regular calculated points	233	779
	Number of irregular calculated points	80	142
	Total number of calculated points	313	921
	Number of edge points	38	52
	Number of boundary points	6	26
	1 and 2: in the case of a computational grid 61 × 61, the two embedding squares
	[-1, 1] × [-1, 1] and [-0.6, 0.6] × [-0.6, 0.6] are compared. The embedding square [-1, 1] × [-1, 1]
	corresponds to a coarse grid with 313 vertices inside the ellipse (5.1). The embedding square
	[-0.6, 0.6] × [-0.6, 0.6] correspond to a fine grid with 779 vertices inside the ellipse (5.1).

Table 5 .

 5

	N × N N ite	e ∞	Rate	(e x ) ∞	Rate	(e y ) ∞	Rate
	20 × 20	20	1.3111(-2)		9.9770(-2)		3.1718(-1)
	40 × 40	80	1.7879(-3) 2.87 1.7643(-2) 2.50 3.5885(-2) 3.15
	80 × 80	320 4.3797(-4) 2.03 4.5494(-3) 1.996 6.9431(-3) 2.37
	160 × 160	1280 1.0803(-4) 2.02 1.1213(-3) 2.021 1.7036(-3) 2.03

  This gives thatu x = (I + h 2 6 δ 2 x ) -1 • δ x u. Let x ∈ [0, 1] → u(x)be a regular function. The truncation error u → t(u) is (6.10)t(u) = (u * ) x -(∂ x u) *

											6	,
	and
	(6.8)	        	Z i k = λ k = 4 sin 2 kπ 2 N sin N	ikπ N , 1 ≤ k ≤ N -1. , 1 ≤ i, k ≤ N -1,
	Proof. The relation (2.18) is rewritten
	(6.9)							(I +	h 2 6	δ 2 x )u x = δ x u.
		I +	h 2 6	δ 2 x	-1 δ x u * -(I +	h 2 6	δ 2 x )(∂ x u) *
		= I +	h 2 6	δ 2				-	h 2 6	δ 2 x (∂ x u) *	.
											(I)	(II)
	The term (I) is expressed as
	(6.11) x u) The term (II) is (I) = h 2 6 (∂ 3 x u) * + h 4 120 (∂ 5
	(6.13)	(II) =	h 2 6	δ 2 x (∂ x u) * =	h 2 6	(∂ 3 x u) * +	h 2 12	(∂ 5 x u) * +	h 4 360	w ,
	where w satisfies
	(6.14) w L Therefore, the truncation error is
	(6.15)	t = (I +	h 2 6	δ 2 x ) -1 -	h 4 180	(∂ 5 x u)

x

-1 (δ x u * -(∂ x u) * ) * + v, where v satisfies (6.12) v L ∞ ≤ Ch 6 (∂ 7 x u) * ) L ∞ . ∞ ≤ C ((∂ 7 x u) * ) L ∞ . * + O(h 6 ) .

In a Birkhoff problem, the set of data carried by each point in the stencil is dependent on the point. We refer to[START_REF] Lorentz | Multivariate Birkhoff Interpolation[END_REF] for a theory of multivariate Birkhoff interpolation.

performed in MAPLE

with MAPLE

Consider now extending (4.7) to the irregular stencil shown in Fig. 4. The polynomial P (x, y) = Matania Ben-Artzi: Institute of Mathematics, The Hebrew University, Jerusalem 91904, Israel E-mail address: mbartzi@math.huji.ac.il Jean-Pierre Croisille: Department of Mathematics, IECL, UMR 7502, Univ. de Lorraine, Metz 57045, France E-mail address: jean-pierre.croisille@univ-lorraine.fr Dalia Fishelov: Afeka -Tel-Aviv Academic College of Engineering, 218 Bnei-Efraim St., Tel-Aviv 69107, Israel E-mail address: daliaf@afeka.ac.il