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

(P1* ^{0}*)

*K*

_{ij}^{(n)}

*≤*0,

*i6=j, i*= 1, ..., N, j = 1, ...,

*N,*¯

(P2* ^{0}*) (X

^{(n)}

_{1})

*=*

_{ij}*M*

_{ij}^{(n)}+

*θ∆t K*

_{ij}^{(n)}

*≤*0,

*i6=j, i*= 1, ..., N, j = 1, ...,

*N,*¯ (P3

*) (X*

^{0}^{(n)}

_{2})

*=*

_{ii}*M*

_{ii}^{(n)}

*−*(1

*−θ)∆t K*

_{ii}^{(n)}

*≥*0,

*i*= 1, ..., N,

(2.3.122)

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

Proof. According to (P1* ^{0}*), 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)}_{10}e_{0} *≥*X^{(n)}_{10} e_{0}+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 (P1* ^{0}*), X

^{(n)}

_{1∂}

*≤*0, hence (P2) is obvious. Finally, (P3

*) guarantees the non-negativity of X*

^{0}^{(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*−θ) max*_{1≤i≤N}*K*_{i,i}^{(n)} (2.3.123)
*is satisfied.*

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

The conditions (P1* ^{0}*) and (P2

*) are obviously true. The assumption (2.3.123) guarantees (P3*

^{0}*).*

^{0}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^{∗}*h*^{2}_{min} + _{h}^{a}^{∗}

min *−a*_{inf}´ = *h*^{2}_{min}

(1*−θ) (2k** ^{∗}*+

*a*

^{∗}*h*

_{min}

*−a*

_{inf}

*h*

^{2}

_{min})

*,*(2.3.124)

*where*

*k*

*= sup*

^{∗}_{(x,t)∈Q}

_{T}*{*P

_{d}*m=1**k** _{m}*(x, t)}; and

*a*

*= sup*

^{∗}_{(x,t)∈Q}

_{T}*{*P

_{d}*m=1**|a** _{m}*(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

*h*_{i+sign((a}

*m*)^{(n)}* _{i}* )x

*m*

*≤*2(k

*)*

_{m}*i+0.5sign((a*

*m*)

^{(n)}

*)*

_{i}*|(a** _{m}*)

^{(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*

*h**max**≤* 2 min1≤m≤dinf*Q**T**k**m*(x, t)

max_{1≤m≤d}sup_{Q}_{T}*|a** _{m}*(x, t)|

*.*(2.3.126)

*The condition* (2.3.123) *can be guaranteed by the assumption*

∆t*≤* 1

(1*−θ)*

³2k^{∗}

*h*^{2}_{min} *−a*_{inf}

´. (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 operator*k**st* = 1,
*a** _{st}* = 0 and

*a*

_{inf}= 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* *−a*0(x, t), (2.3.128)

9We recall that*a*inf= inf*Q**T**a*0, 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 *a*_{inf}*≤*
*a*_{0}(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 operator*L*in (2.3.128) the entries of the stiffness matrix are defined
from the relation

(P1* ^{0}*) is satisfied. Due to the positivity of

*k*

*, for sufficiently small*

_{?}*h, this condition can*always be satisfied. Then, (P2

*) can be satisfied by the condition*

^{0}∆t *≥* 1

Hence, for (P3* ^{0}*) the condition reads as

∆t *≤* 2h^{2}

(1*−θ) (6k** ^{?}*+ 3ha

^{?}*−*2h

^{2}

*a*

*)*

_{inf}*.*(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 parameter*h* *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) =a*_{0}(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 ∆

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

_{d}*L*in the general form (2.3.65). In this part we analyze this problem for a special case:

*a** _{m}*(x, t) =

*a*

_{0}(x, t) = 0 and

*k*

*(x, t) =*

_{m}*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* *−*
X*d*

*m=1*

*∂*^{2}

*∂x*^{2}* _{m}* =

*∂*

*∂t* *−*∆_{d}*.* (2.3.137)

First we assume that *d*= 2 and Ω is a polygonal domain in IR^{2} with a boundary *∂Ω,*
*T >*0. Let Ω be covered by a hybrid mesh *T** _{h}* (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.) Let

*P*

_{1}

*, ..., P*

*denote the interior nodes, and*

_{N}*P*

_{N+1}*, . . . , P*

*N*¯ the boundary ones in

*T*

*. We also define*

_{h}*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

*φ*

*(P*

_{i}*) =*

_{j}*δ*

_{ij}*, i, j*= 1, . . . ,

*N*¯, where

*δ*

*is the Kronecker symbol. For these basis functions we have*

_{ij}*φ*_{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 *T** _{h}* consists of
triangles

*T*

*i*and rectangles

*R*

*j*, 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**∈T**h*

*σ*_{T}*, λ*^{4}* _{min}* = min

*T**∈T**h*

meas_{2}*T, λ*^{4}* _{max}*= max

*T**∈T**h*

meas_{2}*T,* (2.3.140)
*where* meas_{2}*T* *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 min

^{2}

*{a, b} −*max

^{2}

*{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∈T**h*

*µ*_{R}*, λ*^{#}* _{min}* = min

*R∈T**h*

meas_{2}*R, λ*^{#}* _{max}*= max

*R∈T**h*

meas_{2}*R,* (2.3.142)
*where* meas_{2}*R* *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* *T*_{h}*is of compact type if both* *σ* *≥*0 *and*
*µ≥*0. We say that the hybrid mesh *T*_{h}*is of strictly compact type if* *σ >*0 *and* *µ >*0.

In fact, any nonzero entry*K** _{ij}* of the stiffness matrixK, and any nonzero entry

*M*

*of the mass matrixM, presents a sum of several contributions calculated over triangles and rectangles forming the intersection of the supports of basis functions*

_{ij}*φ*

*and*

_{i}*φ*

*. In what follows, we find what these contributions are equal to.*

_{j}The contributions to the mass matrix M over the triangle *T*:
*M*_{ij}*|** _{T}* = meas

_{2}T

12 (i*6=j*), *M*_{ii}*|** _{T}* = meas

_{2}T

6 *.* (2.3.143)

The contribution to the stiffness matrix K over*T* (with the vertices denoted by e.g., *P** _{i}*,

*P*

*, and*

_{j}*P*

*) is equal to*

_{k}*K*_{ij}*|** _{T}* =

*−*1

2ctg *α** _{ij}* (i

*6=j*), (2.3.144)

where *α**ij* is the interior angle opposite the edge *P**i**P**j*. If the triangle *T* is non-obtuse,
then, obviously, *K*_{ij}*|** _{T}* is non-positive. Further,

*K*_{ii}*|** _{T}* =

*l*

^{2}

_{i}4 meas_{2}*T,* (2.3.145)

where *l** _{i}* is the length of the edge of the triangle

*T*opposite the vertex

*P*

*.*

_{i}Now we pass to the investigation of the contributions from the rectangular elements, i.e.,
form the rectangular element *R* with the edges *a** _{R}* and

*b*

*. We can easily show that*

_{R}*M*_{ij}*|*_{R}*∈*
If the rectangle *R* is non-narrow, then *K*_{ij}*|** _{R}* is non-positive. Also,

*K*_{ii}*|** _{R}* =

*a*

^{2}

*+*

_{R}*b*

^{2}

_{R}3 meas_{2}*R.* (2.3.148)

In the following we formulate three lemmas which give sufficient conditions for the
requirements (P1* ^{0}*)-(P3

*) in Theorem 2.3.61, respectively.*

^{0}Lemma 2.3.73 *Let the hybrid mesh* *T*_{h}*be of compact type, then* *K*_{ij}*≤*0 *for* *i6=j, i* =

because the values*K**ij**|**R* and*K**ij**|**T* are non-positive for any non-narrow rectangle and any
non-obtuse triangle, respectively.

Proof. We denote suppφ_{i}*∩*suppφ* _{j}* by

*S, then*

In the formulation of the lemma we have used the notation *l** _{max}* for the length of the
longest edge in

*T*and

*γ*_{max}^{#} = max
Proof. We have the following lower estimations:

*M**ii**−*(1*−θ)∆tK**ii* =

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{ ^{µ}

3λ^{#}*max**,* ^{σ}

*λ*^{4}*max**}* (2.3.153)

*and*

∆t*≤* 1

9 (1*−θ) max{*^{γ}^{max}^{#}

3λ^{#}_{min}*,* ^{γ}^{max}^{4}

4*λ*^{4}_{min}*}* (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*≥* *h*^{2}

3θ (2.3.156)

*and*

∆t *≤* *h*^{2}

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 Ω*⊂*IR* ^{d}* is a polytopic domain with a boundary

*∂Ω.*

^{10}Let Ω be covered by a simplical mesh

*T*

*. For this case we are also able to define the elements of the local mass and stiffness matrices [20].*

_{h}The contributions to the mass matrix M over the simplex *T*:
*M*_{ij}*|** _{T}* = 1

(d+ 1)(d+ 2)meas_{d}T, (i*6=j),* *M*_{ii}*|** _{T}* = 2

(d+ 1)(d+ 2)meas_{d}T. (2.3.158)
The contribution to the stiffness matrix K over the simplex*T* is equal to

*K*_{ij}*|** _{T}* =

*−*(meas

_{d−1}S

_{i})(meas

_{d−1}S

_{j})

*d*^{2}meas_{d}T cos*γ*_{ij}*,* (i*6=j*), K_{ii}*|** _{T}* = (meas

_{d−1}S

_{i})

^{2}

*d*^{2}meas_{d}T *.*(2.3.159)
Here*T* is a simplex with vertices*P*_{1}*, . . . , P** _{d+1}*,

*S*

*is the (d*

_{i}*−*1)-dimensional face opposite to the vertex

*P*

*i*, and cos

*γ*

*ij*is the cosine of the interior angle between faces

*S*

*i*and

*S*

*j*. In the following we formulate those conditions which guarantee the requirements in The-orem 2.3.61.

10Mostly we assume that Ω is a*d-dimensional rectangle, which is also called as orthotope, *
hyperrec-tangle, or*d-dimensional box.*

(P1* ^{0}*) When the condition
Hence, a sufficient condition of (P2

*) is that each term in the above sum is non-positive, i.e.,*

^{0}A sufficient condition of (P3* ^{0}*) is that each term in the above sum is non-negative,
i.e.,
Let us introduce the notations:

*S**?* = min*T**⊆T**h*(measd−1Si), *S** ^{?}* = max

*T*

*⊆T*

*h*(measd−1Si), (2.3.165)

*T*

*= min*

_{?}

_{T}

_{⊆T}*(meas*

_{h}_{d}T

_{i}),

*T*

*= max*

^{?}

_{T}

_{⊆T}*(meas*

_{h}_{d}T

_{i}), (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* *T*_{h}*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*
*d*^{2}
*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* *T**h* *the number*
*µ** _{d}*(T

*) = 1*

_{h}*γ*_{?}^{(d)}
µ*S*^{?}

*S*_{?}

¶_{2} µ
*T*^{?}*T*_{?}

¶_{2}

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

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

1*−θ* (2.3.170)

is true. Let us notice that (2.3.170) can be written as
*θ≥θ*_{0}^{(d)} := 2*µ** _{d}*(T

*)*

_{h}2 +*µ** _{d}*(T

*)*

_{h}*.*(2.3.171)

Remark 2.3.81 *We note that on the uniform mesh* *T*_{h}*we have* *S** ^{?}* =

*S*

_{?}*and*

*T*

*=*

^{?}*T*

_{?}*.*

*Hence, for this case*

*µ** _{d}*(T

*) = 1*

_{h}*γ*

*?*

^{(d)}

and *θ*^{(d)}_{0} = 1

2γ*?*^{(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 dom*L.) Applying some approximation to the operator* *L, we obtain the*
discrete mesh operator *L, and hence we get the equation*

*Lν* =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}*=

*r*g

*(q,M*

_{stab}^{(n)}

_{0}

*,*K

^{(n)}

_{0})ν

^{n−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. When* Ω *is unbounded, then* M^{(n)}_{0} *and* K^{(n)}_{0} *are infinite matrices.*

The function *r*g* _{stab}*(q,M

^{(n)}

_{0}

*,*K

^{(n)}

_{0}) is called the

*stability function of the method*and it is derived in the following way: for the numerical solution of the semidiscretized problem

M^{(n)}_{0} *y** ^{0}*(t) = 1

*h*^{2}K^{(n)}_{0} *y(t);* M^{(n)}_{0} *y(0) =y*_{0}*,* (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 method^{12} corresponds to the choice *θ* = 0.5, i.e, for the finite
difference approximation its stability function reads as

*r** _{CN}*(qK

_{0}) = (I

*−*0.5qK

_{0})

*(I+ 0.5qK*

^{−1}_{0}), (2.4.6) where, for simplicity, we omit the notation of the dependence of the matrix K

_{0}on

*n.*

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 condition*q≤*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 and*h.*

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.