# Non-negativity preservation for more general discrete mesh operators 48

In document Numerical Treatment of Linear Parabolic Problems (Pldal 49-60)

## 2.3 Discrete analogs of the qualitative properties - reliable discrete models

### 2.3.7 Non-negativity preservation for more general discrete mesh operators 48

In this section we analyze the non-negativity preservation of discrete mesh operators obtained by the discretization of the partial differential operators in a general form. We should guarantee the validity of conditions (P1)-(P3) in Theorem 2.3.15. The following statement gives a sufficient condition for this in terms of the matrices M(n) and K(n). Theorem 2.3.61 Under the following conditions

(P10) Kij(n) 0, i6=j, i = 1, ..., N, j = 1, ...,N,¯

(P20) (X(n)1 )ij =Mij(n)+θ∆t Kij(n)0, i6=j, i = 1, ..., N, j = 1, ...,N,¯ (P30) (X(n)2 )ii=Mii(n)(1−θ)∆t Kii(n)0, i= 1, ..., N,

(2.3.122)

the non-negativity assumptions (P1)-(P3) in Theorem 2.3.15 are satisfied.

Proof. According to (P10), the elements of (X(n)1 )i,j for each i6=j are non-positive.

Moreover, by using the basic estimations in Section 2.3.5, we obtain X(n)10e0 X(n)10 e0+X(n)1∂e

| {z }

≤0

=X(n)1 e= µ 1

∆tM+θK(n)

e>0.

Thus, the matrices X(n)10 are regular M-matrices and as such they are regular and their inverses are non-negative ([9]). Hence, (P1) is satisfied. Due to (P10), X(n)1∂ 0, hence (P2) is obvious. Finally, (P30) guarantees the non-negativity of X(n)2 , which, together with (P1), is sufficient for (P3). This completes the proof.

First we investigate the finite difference approximations, and then the linear/bilinear finite element discretization.

a. The non-negativity preservation property of the finite difference mesh operators

Let us consider the finite difference approximation of the differential operator in the general form (2.3.65). In the next theorem an a priori sufficient condition of the non-negativity preservation is given for the upwind finite difference operators.

Theorem 2.3.62 The upwind finite difference discrete mesh operatorL, defined by (2.3.17) with (2.3.72)-(2.3.74) possesses the discrete non-negativity preservation property if the condition

∆t 1

(1−θ) max1≤i≤NKi,i(n) (2.3.123) is satisfied.

Proof. The non-negativity preservation can be guaranteed by the Theorem 2.3.61.

The conditions (P10) and (P20) are obviously true. The assumption (2.3.123) guarantees (P30).

Remark 2.3.63 By using the coefficient functions in the continuous operatorLin (2.3.65), the condition (2.3.123) can be guaranteed by the assumption

∆t 1

(1−θ)³

2k h2min + ha

min −ainf´ = h2min

(1−θ) (2k+ahmin−ainfh2min), (2.3.124) where k = sup(x,t)∈QT{Pd

m=1km(x, t)}; and a = sup(x,t)∈QT{Pd

m=1|am(x, t)|}.9

Next we pass to the consideration of the finite difference discrete mesh operator L defined by (2.3.17) with (2.3.66)-(2.3.69) (central difference approximation). For such an approximation the non-positivity of the elements (X(n)1 )i,j for each i6=j is not automati-cally satisfied, only if the condition

hi+sign((a

m)(n)i )xm 2(km)i+0.5sign((am)(n)i )

|(am)(n)i | (2.3.125) holds. Then, by repeating the proof of Theorem 2.3.62, we arrive at the following state-ment.

Theorem 2.3.64 The finite difference discrete mesh operator L, defined by (2.3.17) and (2.3.66)-(2.3.69), under the condition (2.3.125), possesses the discrete non-negativity preservation property if the condition (2.3.123) is satisfied.

Remark 2.3.65 By using the coefficient functions in the continuous operatorLin (2.3.65), the condition (2.3.125) can be guaranteed by the assumption

hmax 2 min1≤m≤dinfQTkm(x, t)

max1≤m≤dsupQT|am(x, t)|. (2.3.126)

The condition (2.3.123) can be guaranteed by the assumption

∆t 1

(1−θ)

³2k

h2min −ainf

´. (2.3.127)

It is worth characterizing the sharpness of the estimates (2.3.124) and (2.3.127) on the heat conduction operator, for which we know the exact bounds. For this operatorkst = 1, ast = 0 and ainf = 0, hence, both the estimations (2.3.124) and (2.3.127) result in the bound (2.3.109), which is the exact bound for all N = 1,2, . . ..

b. The non-negativity preservation property of the linear finite element mesh operators

In the sequel we analyze the discrete mesh operator obtained by linear finite element space and theθ-method time dicretization. First we investigate the one dimensional case, and then the higher dimensional case will be considered.

b1. Non-negativity preservation property in 1D case

In this part we give the conditions of the non-negativity preservation for the linear finite element discretization on uniform mesh of the operator

L≡

∂t−

∂x(k(x, t)

∂x)−a(x, t)

∂x −a0(x, t), (2.3.128)

9We recall thatainf= infQTa0, according to the proof of Theorem 2.2.12.

We assume again that the bounds k k(x, t) k > 0; a ≥ |a(x, t)| and ainf a0(x, t)0 are finite.

The non-negativity preservation can be guaranteed again by the use of Theorem 2.3.61.

Due to (2.3.78), for the operatorLin (2.3.128) the entries of the stiffness matrix are defined from the relation

(P10) is satisfied. Due to the positivity of k?, for sufficiently small h, this condition can always be satisfied. Then, (P20) can be satisfied by the condition

∆t 1

Hence, for (P30) the condition reads as

∆t 2h2

(1−θ) (6k?+ 3ha?2h2ainf). (2.3.135) Hence, we can summarize our results as follows.

Theorem 2.3.66 The discrete mesh operator L, obtained by linear finite element dis-cretization of the operator L in (2.3.128), possesses the discrete non-negativity preserva-tion property if, according to (2.3.132), the space discretization parameterh is sufficiently small, and the time discretization parameter ∆t satisfies both the lower and upper bounds, given in (2.3.133) and (2.3.135), respectively.

Remark 2.3.67 We note that for a special case, namely, for the heat conduction operator L in (2.3.81), the bounds (2.3.133) and (2.3.135) turn into the exact bound (2.3.120), which is valid for all N.

Remark 2.3.68 Naturally, the lower bound should not exceed the upper bound. This results in certain restrictions for the possible choice of the parameter θ. Namely, if we introduce the notation µ=µ(k) :=k?/k? for the oscillation of the function k(x, t), then, under the condition

θ≥θ0(µ) = µ

2 +µ, (2.3.136)

for sufficiently small h we can always find suitable values for ∆t. Since µ≥1, therefore θ0 1/3 and θ0 = 1/3 only for a(x, t) =a0(x, t) = 0 and k(x, t) =constant. Since θ0(µ) in (2.3.136) tends to one monotonically as µ tends to infinity, therefore, for any linear finite element discretization with sufficiently small h, there always exists suitable θ0, such that for any θ 0,1], under the conditions (2.3.133) and (2.3.135), the discrete mesh operator is non-negativity preserving. It is worth mentioning that the Crank-Nicolson scheme (θ = 0.5) belongs to this interval only when µ 2, which corresponds to the condition k? 2k?.

b2. Non-negativity preservation property in higher dimensions

To give the exact condition for the non-negativity preservation property of the discrete mesh operator, obtained by the finite element method, for dimensions d 2 is a diffi-cult task, even the simplest case L = ∆d, where ∆d denotes the d-dimensional Laplace operator. (The problem is the inversion of the block tridiagonal matrices.) The prob-lem turns into more complex task when we consider the discretization of the operator L in the general form (2.3.65). In this part we analyze this problem for a special case:

am(x, t) = a0(x, t) = 0 and km(x, t) = constant. (Without loss of generality we as-sume that this constant equals one.) We give sufficient condition for the non-negativity preservation property.

Hence, we consider the operator L≡

∂t Xd

m=1

2

∂x2m =

∂t d. (2.3.137)

First we assume that d= 2 and Ω is a polygonal domain in IR2 with a boundary ∂Ω, T >0. Let Ω be covered by a hybrid mesh Th (see Figure 2.3.3), where h stands for the discretization parameter. (We note that in extremal cases the hybrid mesh covers two special types of meshes - a triangular one and a rectangular one.) LetP1, ..., PN denote the interior nodes, and PN+1, . . . , PN¯ the boundary ones in Th. We also define N := ¯N−N.

Let φ1, . . . , φN¯ be basis functions defined as follows: each φi is required to be contin-uous piecewise linear (over triangular elements) and bilinear (over rectangular elements) such that φi(Pj) =δij, i, j = 1, . . . ,N¯, whereδij is the Kronecker symbol. For these basis functions we have

φi 0, i= 1, ...,N,¯ (2.3.138)

N¯

X

i=1

φi 1 in Ω, (2.3.139)

i.e., the assumptions (2.3.75) and (2.3.80) are satisfied. As before, we will apply Theorem 2.3.61 to prove the non-negativity preservation property of such a linear/bilinear finite element discrete mesh operator.

Figure 2.3.3: An example of a hybrid mesh.

To do this, we give some preliminaries. The hybrid unstructured mesh Th consists of triangles Ti and rectangles Rj, which together cover the solution domain Ω (see Figure 2.3.3).

Definition 2.3.69 LetT be a triangle with interior anglesα1, α2, andα3. Let the number σT be defined asσT = min{ctg α1, ctgα2, ctgα3}. We say that the triangle is non-obtuse if σT 0. We say that the triangle is of acute type if σT > 0. We also introduce the following parameters:

σ = min

T∈Th

σT, λ4min = min

T∈Th

meas2T, λ4max= max

T∈Th

meas2T, (2.3.140) where meas2T denotes the area of the triangle T.

Definition 2.3.70 Let R be a rectangle with the length of edges a and b. Let us define the number

µR= 2 min2{a, b} −max2{a, b}

ab . (2.3.141)

We say that the rectangle is non-narrow if µR 0. We call the rectangle strictly non-narrow type if µR>0. We also introduce the following parameters:

µ= min

R∈Th

µR, λ#min = min

R∈Th

meas2R, λ#max= max

R∈Th

meas2R, (2.3.142) where meas2R denotes the area of the rectangle R.

Remark 2.3.71 A rectangle is non-narrow if its longest edge is not greater than 2 times the shortest one.

Definition 2.3.72 We say that the hybrid mesh Th is of compact type if both σ 0 and µ≥0. We say that the hybrid mesh Th is of strictly compact type if σ >0 and µ >0.

In fact, any nonzero entryKij of the stiffness matrixK, and any nonzero entryMij of the mass matrixM, presents a sum of several contributions calculated over triangles and rectangles forming the intersection of the supports of basis functions φi and φj. In what follows, we find what these contributions are equal to.

The contributions to the mass matrix M over the triangle T: Mij|T = meas2T

12 (i6=j), Mii|T = meas2T

6 . (2.3.143)

The contribution to the stiffness matrix K overT (with the vertices denoted by e.g., Pi, Pj, and Pk) is equal to

Kij|T =1

2ctg αij (i6=j), (2.3.144)

where αij is the interior angle opposite the edge PiPj. If the triangle T is non-obtuse, then, obviously, Kij|T is non-positive. Further,

Kii|T = l2i

4 meas2T, (2.3.145)

where li is the length of the edge of the triangle T opposite the vertex Pi.

Now we pass to the investigation of the contributions from the rectangular elements, i.e., form the rectangular element R with the edges aR and bR. We can easily show that

Mij|R If the rectangle R is non-narrow, then Kij|R is non-positive. Also,

Kii|R = a2R+b2R

3 meas2R. (2.3.148)

In the following we formulate three lemmas which give sufficient conditions for the requirements (P10)-(P30) in Theorem 2.3.61, respectively.

Lemma 2.3.73 Let the hybrid mesh Th be of compact type, then Kij 0 for i6=j, i =

because the valuesKij|R andKij|T are non-positive for any non-narrow rectangle and any non-obtuse triangle, respectively.

Proof. We denote suppφisuppφj by S, then

In the formulation of the lemma we have used the notation lmax for the length of the longest edge in T and

γmax# = max Proof. We have the following lower estimations:

Mii(1−θ)∆tKii =

This completes the proof.

From Lemmas 2.3.73–2.3.75 it follows immediately

Theorem 2.3.76 The linear / bilinear finite element heat conduction discrete mesh oper-ator on a hybrid mesh of strictly compact type is non-negativity preserving if the conditions

∆t 1

6θ min{ µ

#max, σ

λ4max} (2.3.153)

and

∆t 1

9 (1−θ) max{γmax#

#min, γmax4

4λ4min} (2.3.154)

are fulfilled.

Remark 2.3.77 For pure rectangular meshes we have the weaker lower bound for ∆t in the form

∆t λ#max

3θ µ, (2.3.155)

because in (2.3.151) we can apply a weaker estimation [44]. For a square mesh with the step-size h the sufficient condition of the non-negativity preservation is

∆t h2

3θ (2.3.156)

and

∆t h2

6(1−θ), (2.3.157)

respectively. This shows that in this case the non-negativity preservation can be guaranteed only (with our sufficient condition) for methods with θ 2/3, i.e., the Crank-Nicolson scheme is not included.

Now we consider the operator L, defined by (2.3.137), in higher dimensions i.e., for d≥ 3. We assume again that ΩIRd is a polytopic domain with a boundary ∂Ω.10 Let Ω be covered by a simplical meshTh. For this case we are also able to define the elements of the local mass and stiffness matrices [20].

The contributions to the mass matrix M over the simplex T: Mij|T = 1

(d+ 1)(d+ 2)measdT, (i6=j), Mii|T = 2

(d+ 1)(d+ 2)measdT. (2.3.158) The contribution to the stiffness matrix K over the simplexT is equal to

Kij|T =(measd−1Si)(measd−1Sj)

d2measdT cosγij, (i6=j), Kii|T = (measd−1Si)2

d2measdT .(2.3.159) HereT is a simplex with verticesP1, . . . , Pd+1,Si is the (d1)-dimensional face opposite to the vertex Pi, and cosγij is the cosine of the interior angle between faces Si and Sj. In the following we formulate those conditions which guarantee the requirements in The-orem 2.3.61.

10Mostly we assume that Ω is ad-dimensional rectangle, which is also called as orthotope, hyperrec-tangle, ord-dimensional box.

(P10) When the condition Hence, a sufficient condition of (P20) is that each term in the above sum is non-positive, i.e.,

A sufficient condition of (P30) is that each term in the above sum is non-negative, i.e., Let us introduce the notations:

S? = minT⊆Th(measd−1Si), S? = maxT⊆Th(measd−1Si), (2.3.165) T? = minT⊆Th(measdTi), T? = maxT⊆Th(measdTi), (2.3.166) γ?(d) = min cosγij. (2.3.167) We can summarize our results as follows.

Theorem 2.3.78 Let us assume that the simplical mesh Th is of the strictly acute type, i.e., the geometrical condition (2.3.160) is satisfied. Then, for ∆t chosen in accordance with the upper and lower bounds (2.3.162)and (2.3.164), the linear finite element discrete mesh operator L is non-negativity preserving.

Corollary 2.3.79 Under the condition d2 the linear finite element discrete mesh operator on the strictly acute simplical mesh is non-negativity preserving.

In practice, we apply the condition (2.3.168). Therefore, it is worth analyzing when it is applicable.

Definition 2.3.80 For an acute simplical mesh Th the number µd(Th) = 1

γ?(d) µS?

S?

2 µ T? T?

2

(2.3.169) is called the uniformity number of the partition.

Clearly,µd(Th)>1 and the sufficient condition (2.3.168) can be applied only when µd(Th) 2θ

1−θ (2.3.170)

is true. Let us notice that (2.3.170) can be written as θ≥θ0(d) := 2µd(Th)

2 +µd(Th). (2.3.171)

Remark 2.3.81 We note that on the uniform mesh Th we have S? = S? and T? = T?. Hence, for this case

µd(Th) = 1 γ?(d)

and θ(d)0 = 1

?(d)+ 1. (2.3.172) In the 2D case this implies the following. Since min(maxγij) = π/3, therefore γ?(2) = 0.5.

Therefore, on the uniform mesh, according to (2.3.172), we get θ0(2) = 0.5, which means that the Crank-Nicolson scheme is applicable. However, as one can see, for the other cases θ(2)0 > 0.5, i.e., the Crank-Nicolson scheme is not included. In the 3D case for the uniform partition we get γ?(3) = cosγij = 1/3. Therefore, on the uniform mesh, according to (2.3.172), we get θ(3) = 0.6. This means that the Crank-Nicolson scheme is not applicable. Finally we remark that with the increase of the uniformity number of the partition, the interval of the admissible values θ is getting shorter and it is approaching the only possible choice θ= 1.

### 2.4 The Crank-Nicolson scheme to the heat equation

In what follows we apply our results to solving the homogeneous heat equation with some prescribed initial condition. (We assume that the boundary conditions, if they exist, are included into domL.) Applying some approximation to the operator L, we obtain the discrete mesh operator L, and hence we get the equation

=0. (2.4.1)

This yields the following algebraic problem: knowing the values of the mesh function ν at the points of the discrete parabolic boundary, using (2.4.1), we seek its values at the interior mesh points. When L is a two-level discrete mesh operator, then the above algebraic problem leads to the one-step iterative method of the form

νn=rgstab(q,M(n)0 ,K(n)0n−1, (2.4.2) for n= 1,2, . . . with given ν0, which is defined from the initial condition.

Remark 2.4.1 For simplicity, for bounded Ω, we assume that the boundary conditions are homogeneous. Whenis unbounded, then M(n)0 and K(n)0 are infinite matrices.

The function rgstab(q,M(n)0 ,K(n)0 ) is called thestability function of the methodand it is derived in the following way: for the numerical solution of the semidiscretized problem

M(n)0 y0(t) = 1

h2K(n)0 y(t); M(n)0 y(0) =y0, (2.4.3) we apply some one-step numerical integration method, which is based on the rational approximation r(z) of the exponential function exp(z).11 When

r(z) =rθ(z) = 1 + (1−θ)z

1−θz , (2.4.4)

with θ [0,1], then the corresponding stability function is e

rθ(q,M(n)0 ,K(n)0 ) =rθ µ

q³

M(n)0 ´−1 K(n)0

= (M(n)0 −θK(n)0 )−1(M(n)0 +(1−θ)K(n)0 ) (2.4.5) for n = 1,2, . . .. Let us notice that (2.4.2) and (2.4.5) are equivalent to (2.4.1) with the choice L as in (2.3.17).

The Crank-Nicolson method12 corresponds to the choice θ = 0.5, i.e, for the finite difference approximation its stability function reads as

rCN(qK0) = (I0.5qK0)−1(I+ 0.5qK0), (2.4.6) where, for simplicity, we omit the notation of the dependence of the matrix K0 onn.

In what follows we consider this method in more details. As it is known, due to its higher time accuracy, it is of a special interest. It is well known that this method is absolute stable, therefore it is convergent for any choice of the discretization parameters.

In this part we analyze the usage of this method for the simplest heat equation operator (2.3.81) in 1D.

According to Table 2.3.2, the conditionq≤1 is a necessary and sufficient condition of the discrete non-negativity preservation, for any number of space partition. Due to Theorem 2.3.10, this condition also guarantees the discrete maximum norm contractivity property.

However, this is not a necessary condition for it. As it is shown in [80] and [69], the exact condition in this case is q≤1.5.

Corollary 2.4.2 According to the above observation, the implication DNP DM NC in Figure 2.3.1 cannot be reversed. We note that the necessary condition of the discrete non-negativity preservation for all numbers of partition (q≤2(2−√

2)) is more restrictive than the above exact condition of the discrete maximum norm contractivity.

Hence, we can summarize the results as follows.

11We attract attention, that, inspite of the previous notation in (2.3.67) and later, here the matrixK0

does not depend on ∆t andh.

12John Crank (1915 - 2006) originally worked in industry on the modelling and numerical solution of diffusion in polymers. In 1943, working with Phyllis Nicolson (1917-1968) on finite difference methods for the time dependent heat equation, he proposed the Crank-Nicolson method which has been incorporated universally in the solving of time-dependent problems since then. Their first result on this method was published in [29].

Theorem 2.4.3 The finite difference Crank-Nicolson mesh operator, which corresponds to the one-dimensional continuous differential operator L in (2.3.81), has the following qualitative properties:

It satisfies the discrete strong maximum principle (and hence all other qualitative properties) for any number of the uniform space partition if and only if the condition q∈(0,1] holds.

It satisfies the discrete strong maximum principle (and hence all other qualitative properties) for a sufficiently large number of the uniform space partition if and only if the condition q∈(0,2(2−√

2)] holds.

For the values q (0,1.5] it is contractive in the maximum norm.

In the sequel, we analyze the qualitative behaviour of the operator “after the death”, i.e., in the case q >1.5.

In document Numerical Treatment of Linear Parabolic Problems (Pldal 49-60)