ANALYSIS OF A FOURTH ORDER SCHEME AND APPLICATION OF LOCAL DEFECT CORRECTION METHOD

Abstract

This paper provides a new application similar to the Local Defect Correction (LDC) technique to solve Poisson problem -u"(x) = f(x) with Dirichlet boundary conditions. The exact solution is supposed to have high activity in some region of the domain. LDC is combined with a fourth order compact scheme which is recently developed in Abbas (Num. Meth. Partial differential equations, 2013). Numerical tests illustrate the interest of this application.

1. Introduction

The Local Defect Correction (LDC) method was introduced byW. Hackbusch [8] for solving elliptic boundary value problem. LDC is a domain decomposition technique in which the local domain fully overlaps the global one. It is a generic iterative algorithm of multiscale type for the resolution of discrete problem. Two or more grids of calculation can be used. We refer to [6] for applications to problems in combustion and numerical simulations of the flow and heat transfer in a glass tank. An analysis of the LDC technique in combination with finite difference discretizations is presented in [4]. In [7] the method is extended to include adaptivity, multilevel refinement, domain decomposition and regredding. The LDC method is combined with finite volume discretizations in [5] and finite elements discretizations in [1,2].

Let’s briefly outline the basic version of the LDC technique. Consider the elliptic boundary value problem

where L is a linear elliptic differential operator, ƒ and g are the source term and Dirichlet boundary condition, respectively. To discretize (1), we first choose a global coarse grid (grid spacing H), which we denote by ΩH. Let LH be a discrete operator approximating the continuous operator L. An initial approximation i = 0, on ΩH is obtained by solving the system

In (2), the vector ƒH contains the source term ƒ and the Dirichlet boundary condition g as well. Suppose that LH is invertible. Moreover, suppose that the exact solution u(x, y) of (1) has a high activity region in some (small) part of the domain. We select a subdomain Ωl ⊂ Ω such that the high activity region of u is contained in Ωl. The subdomain Ωl is discretized by a local fine grid (grid spac-ing h). We denote it by . The fine grid is built such that This means that the coarse grid points that lie in the region of refinement are also points of . In order to formulate a discrete problem on , we should define artificial boundary conditions on Ґ, the interface between Ωl and Ω \ Ωl, (Figure 1). Actually, we use an interpolation operator Ph,H to obtain artificial boundary conditions on Ґ. The operator Ph,H gives the values of the fine grid points on Ґh by interpolating the values of the coarse grid points on the interface ҐH.

FIGURE 1.Discretization in one dimension. The big points represent the coarse grid ΩH, the small points represent the fine domain .

On ∂Ωl \ Ґ, we obtain boundary conditions using the Dirichlet boundary condi-tions g in (1)b. On the fine grid we consider the following discrete problem, at iteration i = 0,

where

In (3), i represents an index of iteration and is a square matrix representing the dependence of the fine grid solution on the artificial boundary conditions on Ґ. In (3), the matrix (assumed to be invertible) is a discrete approximation of the operator L on the subdomain Ωl. In the right term in (3), contains the source term ƒ and the Dirichlet boundary condition g on ∂Ωl \ Ґ given in (1). Notice that Ґ ∩∂Ω can be empty as in the figure 1. If we know the values of the defect dH = LH(u|ΩH )−ƒH, we can use it to find more accurate approximation on the coarse grid. This can be obtained by replacing dH in the right hand side of (2). However, we do not know the exact solution of the exact continuous problem then we cannot calculate dH. We will use calculated on the fine grid to approximate dH. Indeed, for all i, index of iteration, the function is the function on the coarse grid defined by:

For each index of iteration i we define by

Actually, provides an estimate of the local discrestization error at all points of . We observed numerically it is better to use the approximation (6) on a proper subset of only. In particular, the points of the coarse grid near ∂Ωl should be excluded. This leads to introduce a “Safety region” denoted by (Figure 2). The estimation of the local error of discretization of the coarse grid is placed in the right hand side of the equations corresponding to coarse grid points belonging to only. Then we apply the correction step on the coarse

FIGURE 2.The small points consist of points of Ωh, the big ones are the points of the safety region .

As (7) contains estimations of the local error of the discretization on the coarse grid, we expect to be a more accurate approximation than Then we obtain better boundary conditions on Ґ and a solution on the fine local grid by (3) with i = 1. By performing the iterations on the index i, we obtain the following algorithm:

Algorithm 1.1. • Initialization

-Solve the problem(2) on the coarse grid ΩH with i = 0. We obtain the vector

-Solve the problem(3) on the fine local grid with i = 0. We obtain the vector

• Iteration i = 1, 2...

- Compute by (5).

- Compute by (6).

- Re-solve the problem (7) on the coarse grid ΩH. We obtain the vector

- Solve the problem (3) on the local fine grid We obtain the vector

In practice, one iteration is enough to obtain a satisafactory approximation on the composite grid ΩH,h which is the grid formed by the union of the coarse grid and the fine grid, The algorithm 1.1 is an elementary version of the LDC method. Several generalizations are possible:

• Use of several fine grids. Notice that the local problems are independent each other and can be solved simultaneously. Recursive refinement where a local fine grid can be a coarse grid for another local fine grid.

• Use of different discretizations. We refer to [6] for a detailed analysis of the behavior of the convergence for the Poisson problem solved with classic five points finite difference scheme.

 

2. Main results

2.1. Analysis of a fourth order scheme. We recall and analyse a new scheme called hermitian Box-scheme (HB-scheme), for the Poisson problem in one dimension based on the previous work [12]. This scheme appears to be fourth order accurate in practice. It has been successfully extended to two dimensions in [9] and three dimensions in [10]. Let’s recall the principle of this scheme. Consider the one-dimensional Poisson problem on the interval Ω = (a, b) with length L = b − a,

Equation (8) is recast in mixed form:

We lay out on Ω a regular grid xj = a+jh, 0 ≤ j ≤ N with stepsize h = L/N, (Figure 3). The unknowns are denoted by uj ≈ u(xj) and ux,j ≈ u′(xj). The vectors U,Ux ∈ ℝN−1 stand for the unknowns at internal points,

FIGURE 3.Grid in dimension 1. The two boundary points x0, xN (denoted by “◦”). The points xi, i = 1, · · · ,N − 1 are the interior points (denoted by “•”).

In analogy to the original box-scheme [13], the HB-scheme is deduced from the integration of the two equations (9)a,b on a box Kj =]xj-1, xj+1[ of length 2h.

Suppose given the averaged values of the source term ƒ(x) on the box Kj,

In practice, (11) can be approximated using Simpson formula

The conservation equation (9)a becomes

Second, the equation (9)b is integrated on the box Kj. This yields

Approximating the integral in the left-hand side of (14) by Simpson formula suggests the following fourth-order hermitian approximation

Equations (15) require approximations of the derivatives on the boundaries. Here, we consider the following third-order approximations

Relations (16) are obtained using Taylor expansions. In summary, the equations (13), (15), (16) with Dirichlet boundary conditions translate to the following HB-scheme: Find u = (ui)0≤i≤N and v = (vi)0≤i≤N solution of:

This HB-scheme is found numerically to be fourth order accurate as it is shown in the numerical tables.

Lemma 2.1. Suppose that Π0ƒj is approximated by Simpson formula:

and (uj)j∈ℤ, (ux,j)j∈ℤare two sequences verifying

Then (uj)j∈ℤ verifies

Proof. Denote by Ej, j ∈ ℤ, the j-th equation of (19)a. Using the sum

and (19)b, we find that the left term of (21) verifies

The right term of (21) verifies

where the equality (20) for j ∈ ℤ.

We refer to [11] for a detailed proof of the following lemma concerning the maximum principles of the standard three point Laplacian scheme.

Lemma 2.2. Let u = (ui)0≤i≤N. Suppose that then

Similarly, if then

The scheme (17) verifies the following discrete maximum and minimum prin-ciples.

Proposition 2.1 (Principles of maximum and minimum). Let u = (ui)0≤i≤N be the solution of HB-scheme (17). Suppose that ƒi ≤ 0, ∀ 1 ≤ i ≤ N − 1, then

Similarly, if ƒi ≥ 0, ∀ 1 ≤ i ≤ N − 1, u verifies

This means that u attains its maxima and minima on the boundary points x0 and xN.

Proof. Suppose that ƒi ≤ 0, ∀ 1 ≤ i ≤ N − 1 then

By (20),

If N is even, using the maximum principle of Lemma 2.2 we obtain

In the same manner if N is odd, we obtain

It results for all N,

Then to prove (25) it is enough to prove that

Suppose that

We have

This gives ux,0 ≤ ux,2 and ux,1 ≤ ux,3. The HB-scheme (17) verifies

Moreover,

Therefore,

But we have u2 ≤ u0, then

therefore

where (32). Suppose that

The equations (17) give

In addition, using (34), ux,N−2 ≤ ux,N, therefore

as uN−1 ≥ uN−3, and uN ≥ uN−2, we get

We deduce

where

Finally,

where (25). In the same manner we can prove the minimum principle (26).

Corollary 2.1 (existence and uniqueness of solution). The HB-scheme (17) has unique solution

Proof. We have (17) a linear system with 2N equations and 2N unknowns. To prove the existence and uniqueness of the solution, it is enough to prove that Π0ƒi = 0, 1 ≤ i ≤ N − 1, implies u = v = 0. Let u = (ui)0≤i≤N and v = (vi)0≤i≤N such that

Using the discrete maximum principle of proposition 2.1, we obtain

then

which proves that the matrix associated to the HB-scheme (17) has nul kernel. As it is square matrix then it is invertible, where the existence and uniqueness of the solution of HB-scheme.

2.2. Matrix form of hermitian Box-scheme. We summarize the finite-difference and matrix notations used in the sequel.

• The tridiagonal matrix is

• The Simpson matrix Ps ∈ is

where I is the identity matrix of order N − 1.

• The antisymmetric matrix K ∈ given by

• Denoting (ei)1≤i≤N−1 the canonical basis of ℝN−1, the matrices F1, F2 ∈ are defined by

Let us turn now to the matrix form of the scheme. We use the notation UL = u0, UR = uN, ux,L = ux,0, ux,R = ux,N for the boundary values.

Proposition 2.2. For all linear approximation of the derivatives on the bound-ary Ux,L, Ux,R in terms of U, Ux and the values on the boundary UL, UR, there exist matrices such that

Proof. Any linear approximation of ux,0 in terms of u = (ui)0≤i≤N, ux = (ux,i)>)0≤i≤N and u0, uN can be written in general form as:

with

Similarly,

with Multiplying (53) by e1 and (54) by eN−1 we get

Adding the equations (55) together

Using the fact that and we can write (56) as

Finally to complete the proof, let be the matrices

We claim that (16) translates into matrix form as [9]

where the matrices A, B, C are given by

and the vector of derivative is given in terms of U and the boundary conditions by

where the matrices D and ε are

Using (59) and (61), Ux can be eliminated which gives the expression of the HB-scheme in the sole unknown U as

where F = [Π0ƒ1, · · · , Π0ƒN-1]T (Π0ƒj is given in (11)). The matrices are

If needed, the gradient is recovered as a postprocessing by (61).

2.3. Two-grid refinement algorithm. Here, a multiscale technique is com-bined with the hermitian Box-scheme. We consider the problem

The domain Ω is discretized with stepsize H = 1/N such that the discrete coarse points are The coarse grid is ΩH. The unknowns vectors are UH and such that To avoid any confusion between variables, we perform a change of notations for the sequel

The hermitian Box-scheme corresponding to the problem (65) has the following matrix form

with F = F = [Π0ƒ1, Π0ƒ2, Π0ƒN-1]T. By comparing with (2), let

By resolving (67) we obtain the solution vector i = 0, on Ωhc and the global vector solution

The domain Ωl has been choosen such that Ωl = (a1, b1), with a1, b1 two fixed scalars between a and b, a ≤ a1 < b1 ≤ b. The choice of a1 and b1 is adapted such that for all N there exist i0 and j0 between 1 and N − 1 such that

that is to say a1 and b1 are two points of the coarse grid Ωhc. In order that has the same number of points as ΩH, we choose a1 and b1 such that Suppose that we know the values ga1 = u(a1) and gb1 = u(b1) where u(x) is the exact solution. Hence, on the domain Ωl we can write

with ul = u|Ωl and ƒl = ƒ|Ωl. The domain Ωl is discretized with stepsize such that the discrete fine points are with Notice the equality between N and Nƒ when The fine grid is The unknown vectors are such that and 1 ≤ i ≤ N − 1. The interface between Ωl and Ω is ∂Ωl = {a1, b1}. Based on the matrix form (63), the hermtian Box-scheme corresponding to (71) is

with Fl = [Π0ƒ1, Π0ƒ2, ..., Π0ƒN-1]T ∈ ℝN-1. The problem is that we do not know the values ga1 and gb1 on ∂Ωl in (71). To approach these values, we replace it by artificial boundary conditions

with i0 and j0 are the indices of the points a1 and b1 in the coarse grid defined in (70). Here, in particular, the artificial boundary ∂Ωl = {a1, b1} is composed of two points of the coarse grid Ωhc then there is no need to interpolate therefore the matrix is zero. By comparing (72) and (3), let

By solving (72) with i = 0 we get

then we obtain the approximation

Following (5) the vector verifies

The defect vector is calculated by

Finally, the algorithm is the following:

Algorithm 2.1. •Initialization

- Solve the problem (67) on the coarse grid Ωhc. The vector solution is denoted by

- Solve the problem (75) on the local grid The vector solution is denoted by

• Iteration i = 1, 2, ...

- Compute the vector by (76).

- Compute the defect vector by (77).

- Solve on the coarse grid

with ε > 0. The solution vector is denoted by

- Solve the problem (75) on the fine grid. The vector solution is denoted by

2.4. Numerical results. In this part, we display some numerical results prov-ing the interest of the technique presented in the last section. The calculation is performed using Matlab. The CPU time is calculated by the functions tic and toc. Notice that the CPU time can be extremely reduced using FFT or another programming language [3]. The algorithm used so far is the algorithm 2.1 after only one iteration. We intend to study in another paper the behavior of the convergence of this algorithm. In the numerical tables, uex, ux,ex are the exact solution and derivative, u, ux are the computed solution and derivative on the uniform grid and are the computed solution and derivative on the composite domain. We use the discrete norm L∞ to estimate the error:

The convergence rate is calculated using the formula

where eN is the error obtained on a mesh of size N.

TABLE 1.Error and convergence rate for Test 1 on uniform grid

TABLE 2.Error and convergence rate for Test 1 on composite grid

Test 1: We consider the Gaussian function This function has a high activity region around The results obtained on the composite grid are compared with those computed by the HB-scheme on the uniform grid with fine stepsize hƒ. Observe in tables 1 and 2 the same accuracy obtained by the two approaches with less points on the composite grid.

Test 2: Suppose that the exact solution is on Ω = (0, 1). This function has a deep gradient around We choose The numerical results are given in tables 3 and 4.

 

3. Conclusion

This paper provides a new application of the multiscale technique combined with a new high order scheme called HB-scheme. The high order accuracy of this scheme is particularly interesting in applications where physical fields have to be calculated as an outcome of a potential solution of a Poisson problem (electromagnetism, gravitation). The main features of this work are the high order of precision on the boundaries and the calculation of the derivatives. This application illustrates the gain in memory and CPU time using this technique as it is shown in the numerical tables. For instance, we have obtained for Test 1 the same precision using 3072 points instead of 4096. This gain is expected to be more interesting in higher dimensions. A generalization of this application to higher dimensions using cubic spline interpolations, based on the previous results in [9,10] will be a part of future work.

TABLE 3.Error and convergence rate for Test 2 on uniform grid

TABLE 4.Error and convergence rate of Test 2 on composite grid