DISCRETE MAXIMUM PRINCIPLE AND ADEQUATE DISCRETIZATIONS OF LINEAR PARABOLIC PROBLEMS

In this paper, we analyze the connections between the different qualitative properties of numerical solutions of linear parabolic problems with Dirichlet-type boundary condition. First we formulate the qualitative properties for the differential equations and shed light on their relations. Then we show how the well-known discretization schemes can be written in the form of a one-step iterative process. We give necessary and sufficient conditions of the main qualitative properties of these iterations. We apply the results to the finite difference and Galerkin finite element solutions of linear parabolic problems. In our main result we show that the nonnegativity preservation property is equivalent to the maximum-minimum principle and they imply the maximum norm contractivity. In one, two, and three dimensions, we list sufficient a priori conditions that ensure the required qualitative properties. Finally, we demonstrate the above results on numerical examples.

1. Introduction. Let us consider the partial differential equation Lv = f , where L denotes some linear partial differential operator, f is a given function, and v is the unknown function to be determined. A function g is also given, which prescribes the values of v on certain parts of the boundary of the solution domain.
Several problems for partial differential equations possess some characteristic qualitative properties, which are typical of the phenomenon the partial differential equation describes. The most important three of them are the maximum-minimum principle (MP), the nonnegativity preservation (NP) and the maximum norm contractivity (MNC). The MP states that the solution of a differential equation can be estimated from above and below by means of the given functions f and g. For example, when f = 0, the solution takes its extremal values on the boundary of the solution domain as well. The NP property means that the nonnegativity of f and g implies the nonnegativity of the solution v. The MNC property, which makes sense for time-dependent problems only, says that for arbitrary two initial functions the maximum norm of the difference of the solutions at every time level is not greater than the maximum norm of the difference of the initial functions.
In most cases, a partial differential equation cannot be solved in a useful form. In some particular cases, the solution can be written in the form of an infinite Fourier series but this is also useless from the practical point of view. Thus, the use of numerical methods is a necessary step to obtain an approximate solution. It is a natural requirement of an adequate numerical method that it has to possess the discrete equivalents of the qualitative properties the continuous problem satisfies.
The discrete version of the maximum-minimum principle is the so-called discrete maximum-minimum principle (DMP). The topic of validity of various DMPs arose already 30 years ago and was first investigated for elliptic problems (see e.g. [4,5,18,23]). The DMP was generally guaranteed by some geometrical conditions for the meshes. The DMP for parabolic problems was discussed in [8,10,11,20]. In [11], based on the acuteness of the tetrahedral meshes, a sufficient condition of the DMP was obtained for the Galerkin finite element solution of certain parabolic problems. In paper [8], a necessary and sufficient condition of the DMP was derived for Galerkin finite element methods and sufficient conditions were given for hybrid meshes. About the DMP, a comprehensive survey can be found in papers [2] and [3]; [17] investigates nonlinear problems.
The conditions of the discrete nonnegativity preservation (DNP) was discussed in [7,12] for linear finite elements in one, two and three dimensions, and in [6] in one dimension with the combination of the finite difference and finite element methods. The DNP is investigated for nonlinear problems in [24].
The discrete maximum norm contractivity (DMNC) was analyzed for one-dimensional parabolic problems in [13,19]. In both references the necessary and sufficient conditions were given. In the first one, the dependence on the spatial discretization was also discussed.
For one-dimensional problems, we can deduce some other remarkable qualitative properties such as the preservation of the shape and the monotonicity of the initial function, and the sign-stability (see e.g. [9,14,15,16,22]).
Albeit some sufficient conditions (and in certain cases the necessary ones as well) are known for the above listed three main qualitative properties, the relations and the implications between them have not been revealed yet. The main goal of this paper is to establish the connections between the maximum-minimum principle, the nonnegativity preservation and the maximum norm contractivity both for linear parabolic problems and their numerical solutions. In this paper, we will show the implications (D)N P ⇐⇒ (D)M P =⇒ (D)M N C, that is the maximum-minimum principle is equivalent to the nonnegativity preservation property and both imply the maximum norm contractivity. We give necessary and sufficient conditions for the adequate (qualitatively correct) discrete models. Furthermore, we also give some useful sufficient conditions. The paper is organized as follows. In §2, we discuss the qualitative properties of continuous problems. In §3, we give the linear algebraic form of the finite difference and finite element discretization of the investigated problem. We prove the existence and the uniqueness of the numerical solution. In §4, we construct a one-step iterative process, we define its qualitative properties and we formulate necessary and sufficient conditions of their preservation. In §5, we apply the results of §4 to the discretizations given in §3. In §6, some sufficient conditions are listed in order to guarantee the numerical qualitative properties a priori. In the last section, we demonstrate our results on numerical examples.
For simplicity, we denote zero matrices and zero vectors by the symbol 0, whose size is always chosen according to the context. The ordering relation for vectors and matrices is always meant elementwise.
for any arbitrary positive number τ . For some fixed number T > 0, we consider the linear partial differential operator where δ ≥ 1, ς 1 , . . . , ς d denote nonnegative integers, |ς| is defined as |ς| = ς 1 + . . . + ς d for the multi-index ς = (ς 1 , . . . , ς d ), and the coefficient functions a ς : Q T → IR are bounded on the set Q T . We define the domain of the operator L, denoted by dom L, as the space of functions v ∈ C(Q T ∪ Γ T ), for which the derivatives D ς v (0 < |ς| ≤ δ) and ∂v/∂t exist in Q T and they are bounded. It can be seen easily that Lv is bounded on Qt 1 for each v ∈ dom L and 0 < t 1 < T . Thus inf Qt 1 Lv and sup Qt 1 Lv are finite. Definition 2.1. We say that the operator L satisfies the maximum-minimum principle if for any function v ∈ dom L the inequality The following statement shows that, in case of the maximum-minimum principle, the functions from dom L are uniquely determined by their values on the boundary Γ T .
Theorem 2.2. Let t 1 ∈ (0, T ) be any arbitrary fixed number. If L satisfies the maximum-minimum principle, then for any two Proof. Let us choose an arbitrary number satisfying the condition 0 < t 0 < t 1 , and we consider the functionv = v − v . For this function the relations Lv| Qt 0 = 0 andv| Γ t 0 = 0 are valid. Thus, based on the maximum-minimum principle we havē v(x, t 0 ) = 0 for all x ∈ Ω and 0 < t 0 < t 1 . This completes the proof. Definition 2.3. The operator L is called nonnegativity preserving when the relations min Γt 1 v ≥ 0 and Lv| Qt 1 ≥ 0 imply v| Qt 1 ≥ 0 for all 0 < t 1 < T .
Definition 2.4. The operator L is called contractive in maximum norm when for all arbitrary two The main connections between the above properties of the operator L are listed in the next two theorems.
Theorem 2.5. The operator L defined in (2.1) satisfies the maximum-minimum principle if and only if it preserves the nonnegativity.
Proof. The necessity of the condition is trivial. In order to show the sufficiency, we choose an arbitrary function v ∈ dom L and apply the operator L to the function v = v − min Γt 1 v − t · min{0, inf Qt 1 Lv}. Clearly,v| Γt 1 ≥ 0. Moreover, we obtain that which implies thatv is nonnegative on Qt 1 by virtue of the nonnegativity preservation assumption. Thus the lower estimation the upper bound is proved similarly. This completes the proof.
Theorem 2.6. If the operator L defined in (2.1) is nonnegativity preserving, then it is contractive in maximum norm.
Proof. Let v and v ∈ dom L be two arbitrary functions with Lv For these functions, the estimations Lv ± | Qt 1 = 0 ≥ 0 and min Γ t 1v ± ≥ 0 are true, which implies the nonnegativity ofv ± on Qt 1 . Thus, we have This completes the proof.
In the sequel, we consider the initial-boundary value problem in the form where g : Γ T → IR, f : Q T → IR are arbitrary given functions. The operator L is defined in (2.1). (In the usual context, g| Ω×{0} is the initial and g| ∂Ω×[0,T ] is the boundary condition, respectively.) We are interested in finding a function v ∈ dom L that satisfies equalities (2.2)-(2.3). We say that the problem (2.2)-(2.3) is nonnegativity preserving / contractive in maximum norm / satisfies the maximum-minimum principle if the operator L in the equation (2.2) possesses the same properties. According to Theorem 2.5 and Theorem 2.6, the maximum-minimum principle is valid for the problem (2.2)-(2.3) if and only if it is nonnegativity preserving; in this case the problem is contractive in maximum norm as well. If the maximum-minimum principle is valid for the problem, then its solution (when it exists) is unique (cf. Theorem 2.2).
In this paper we investigate the parabolic problem where C : Ω → IR is a known bounded function with the property 0 < C min ≤ C ≡ C(x) ≤ C max , the known bounded function κ : Ω → IR has continuous first derivatives and fulfills the property 0 < κ min ≤ κ ≡ κ(x) ≤ κ max , the function g : Γ T → IR is continuous on Γ T , the function f : Q T → IR is bounded in Q T , ∇ denotes the usual nabla operator, and the solution v is sought in C 2,1 (Q T ∪ Γ T ).
Remark 2.7. The problem (2.4)-(2.5) is generally applied for the description of heat conduction processes, where v denotes the temperature, C is the product of the specific heat and the density, κ is the thermal conductivity and f gives the density of heat sources. The variables x and t play the role of the space and time variables, respectively. Problem (2.4)-(2.5) is suitable to describe diffusion and transport phenomena as well.
We show that the problem (2.4)-(2.5) is nonnegativity preserving, contractive in maximum norm and satisfies the maximum-minimum principle. Because (2.4)-(2.5) can be written in the form of the problem (2.2)-(2.3) dividing both sides of the equation by the positive function C (that is L ≡ ∂/∂t − (1/C)∇(κ∇)), we have to show only the validity of the nonnegativity preservation (see Theorem 2.5 and Theorem 2.6). The proof will be based on the next theorem.
Theorem 2.8. For any t 1 ∈ (0, T ) and x ∈ Ω, the solution v of the problem (2.4)-(2.5) satisfies the inequality sup λ>0 e λt 1 min min Proof. For any arbitrary number λ > 0 we define the functionv(x, t) = v(x, t)e −λt , where v stands for the solution of (2.4). It can be seen easily thatv in Q T . Asv is continuous onQ t 1 , it takes its maximum either on the boundary Γ t 1 or in Ω × (0, t 1 ] at some point (x 0 , t 0 ). In the first case we trivially have In the second case, the relations hold for i = 1, . . . , d, where the last two relations imply that ∇(κ∇v)(x 0 , t 0 ) ≤ 0. This relation, combined with the second one in (2.9) and with (2.7), results in Thus, in general case, using the upper bounds (2.8) and (2.10) we obtain the estimationv Multiplying both sides by e λt 1 and taking into account that the relation is true for all positive numbers λ > 0, we obtain the inequality on the right-hand side of (2.6). The lower bound can be proved similarly.
Remark 2.9. The proof of the above theorem is based on the proof of Theorem 2.1 in [20]. Let us notice, however, that we did not confine ourselves to second order operators, and leaving the zeroth order derivative out of the expression of L, we arrived at a stronger estimation.
Theorem 2.10. The problem (2.4)-(2.5) preserves the nonnegativity, and consequently it is contractive in maximum norm and satisfies the maximum-minimum principle.
Proof. The proof is based on the previous theorem. Let t 1 ∈ (0, T ) be an arbitrary number, and let us suppose that 3. Numerical Models of Linear Parabolic Problems. The two most widely used numerical methods for solving the problem (2.4)-(2.5) are the finite difference and the Galerkin finite element methods. Finite difference methods are typically applied when the solution domain Ω is relatively simple (an interval, a rectangle or a block), while finite element methods can be used also with more complicated geometrical structures (e.g. polyhedrons). Both methods start with the discretization of the computational space Ω ⊂ IR d in order to get a system of ordinary differential equations. This is the so-called semi-discretization. Then the system of ordinary differential equations is solved by some time-integration method, such as the wellknown θ-method. In this paper we consider only the 3D case; the 1D and 2D cases can be derived in the same manner by omitting the corresponding terms.

Spatial Semi-Discretization with the Finite Difference Method.
In the case of the finite difference method, the function v is approximated at the points of a rectangular mesh defined on the rectangular domain Ω. Let us denote the interior mesh points by P 1 , . . . , P N ∈ Ω, and the points on the boundary by P N +1 , . . . , P N +N ∂ ∈ ∂Ω, respectively. For the sake of simplicity, we also use the notation P i+x (P i−x ) for the grid point adjoint to P i in positive (negative) x-direction. We also defineN = N +N ∂ . After denoting the semi-discrete approximation of v(x, t) at a grid point x = P i by v i (t), we can approximate (2.4) at each inner point P i of Ω as of the function v at the point P i and at the time instant t, and has the form The distances between the points P i+x , P i and P i−x , P i are denoted by h i+x and h i−x , respectively. The values κ i−x/2 and κ i+x/2 denote the approximate values of the material parameter κ on the segments [P i−x , P i ] and [P i , P i+x ] (typically the midpoint values), C i denotes the approximate value of C at the point P i (typically C i = C(P i ) for continuous functions). The terms v iyy (t) and v izz (t) are defined similarly.
Hence the semi-discretization (3.1) and the boundary condition (2.5) result in a Cauchy problem for the system of ordinary differential equations and M and K are sparse N ×N matrices. For the entries of M we have and the entries of K can be computed using equations (3.1) and (3.2). The nonzero 3.2. Semi-Discretization with the Galerkin Finite Element Method. In the case of the Galerkin finite element method we cover the domain Ω with the following meshes T h (h is a discretization parameter): the 1D mesh consists of intervals, the 2D mesh is a so-called hybrid mesh, which is a combination of triangles and rectangles, and the 3D mesh is a tetrahedron or a block mesh. Figure 3.1 shows a 2D hybrid mesh. As before, P 1 , . . . , P N denote the interior nodes of T h , and P N +1 , . . . , PN the boundary ones. Let φ 1 , . . . , φN be basis functions defined as follows: each φ i is required to be continuous, piecewise linear (over intervals, triangular or tetrahedral elements) or multi-linear (over rectangular or block elements), such (Multi-linearity means that it is bilinear in case of rectangles, and trilinear in case of blocks.) It is obvious that such basis functions have the properties We search for the semi-discrete solution of (2.4)-(2.5) in the form Using the weak formulation of the problem, we arrive at a Cauchy problem again that has the form (3.1) with (see [8]). The above defined matrices M and K are called mass and stiffness matrices, respectively.

Fully Discretized Model.
To get a fully discrete numerical scheme, we choose a time-step ∆t > 0. We denote the approximation to v(n∆t) by v n , and we set f n = f (n∆t) (n = 0, 1, ...). The time-discretization of (3.1) with the θ-method (θ ∈ [0, 1] is a given parameter) can be written in the form of the systems of linear algebraic equations where n = 0, 1, . . . and θ ∈ [0, 1] is a given parameter. Clearly, (3.8) can be rewritten as In the sequel, the matrices M + θ∆tK and M − (1 − θ)∆tK will be denoted by A and B, respectively. We shall use the following partitions of the matrices and vectors: Similar partition is used for the matrices M and K. Then, iteration (3.9) can be written as The relation (3.10) defines a one-step iterative process with respect to the unknown vector u n+1 .
Remark 3.1. For the finite difference method, the choice θ = 0 results in an explicit scheme, due to the diagonality of M. This is the so-called explicit Euler method. The choice θ = 1 gives the so-called implicit Euler method and θ = 1/2 is the Crank-Nicolson method. However, the schemes obtained by the Galerkin finite element method are always implicit for any choice of θ.
In order to guarantee the existence and the uniqueness of u n+1 in (3.10), we have to show that A 0 is a nonsingular matrix.
Theorem 3.2. The matrix A 0 ∈ IR N ×N is nonsingular for both the finite difference and the Galerkin finite element methods.
Proof. The case of the finite difference method: Because of the relation A 0 = M 0 + θ∆tK 0 and the equalities in (3.6), the off-diagonal elements of A 0 are nonpositive and the diagonal ones are positive. Moreover, considering (3.6) and (3.7) we have Thus, the estimate is true (i = 1, . . . , N ), which shows that A 0 is a strictly diagonally dominant matrix, hence it is nonsingular. The case of the Galerkin finite element method: The elements of A 0 are The right-hand side defines a scalar product on the vector space span{φ 1 , . . . , φ N }. Thus A 0 is a Gram-matrix of the linearly independent functions φ 1 , . . . , φ N , hence it is nonsingular and also positive definite. Remark 3.3. It is important to notice that in the case of the finite difference method, A 0 is a so-called M -matrix (see [1, p. 137]). This can be seen from the sign-structure and the strict diagonal dominance of A 0 . That is, in addition to the nonsingularity, A 0 has a nonnegative inverse too. In the finite element method the nonnegativity of A −1 0 is not satisfied automatically. 4. Qualitative Properties of One-Step Iterative Models. As we mentioned in the Introduction, it is a natural requirement for the numerical solution that it has to possess some basic qualitative properties. The numerical solution can be obtained by the iteration (3.10). Hence, the qualitative properties of the numerical solution will be defined as the qualitative properties of such iteration processes. We

PO-iterations.
Based on the iteration (3.10), we will investigate the general one-step iteration u n g n + h n , n = 0, 1, . . . Here h n ≡ [h n 1 , . . . , h n N ] ∈ IR N and all matrices and vectors have the same dimension as in the previous section. We emphasize that now they are assumed to be arbitrary with the only restriction that the matrix A 0 is nonsingular. Iteration (4.2) is called partitioned one-step iteration or shortly PO-iteration.

.2. A PO-iteration satisfies the DNP if and only if the inequalities
The statement follows directly from the following equivalent form of (4.1) is valid, whereû n+1 andũ n+1 are computed from (4.1) by setting u n =û n and u n =ũ n .

Theorem 4.4. A PO-iteration satisfies the DMNC if and only if
We apply (4.3) with the two different vectors u n =û n and u n =ũ n . Hence we obtainû The relation (4.4) is valid for allû n andũ n vectors if and only if A −1 0 B 0 ∞ ≤ 1. Proof. Choosing the vectors h n = 0, u n = e 0 , g n = g n+1 = e ∂ , (4.5) results in u n+1 = e 0 . That is Ae = Be in view of (4.1), which completes the proof. Proof. The DNP ensures the relations −A −1

Definition 4.5. A PO-iteration is said to satisfy the discrete maximum-minimum principle (DMP) if the relation
This shows that A −1 0 B 0 ∞ ≤ 1, that is the DMNC is valid. Remark 4.9. Because the DMP implies both the DNP and condition Ae = Be, the DMP implies the DMNC.
by virtue of the nonnegativity of the vector A −1 0 Be. This completes the proof. 4.3. Qualitatively adequate one-step iterations. We have seen in the previous section that the properties are necessary conditions of the DMP for PO-iterations. In the next theorem, we prove that they are sufficient conditions as well. In order to make the expressions much simpler, we introduce some notations. We define the values v n min = min{v n }, v n max = max{v n },

Theorem 4.11. A PO-iteration satisfies the DMP if and only if it satisfies the conditions (P 1) − (P 3).
Proof. The necessity of the condition follows directly from Theorem 4.2 and Lemmas 4.6 and 4.7. To show the sufficiency, we first prove the inequality on the right-hand side in (4.5). It follows from (P 3) that Bv n max = Av n max . Using (P 1), we Regrouping the above inequality, we get Hence, for the i-th coordinate of both sides we obtain where we applied property (P 2) and Lemma 4.10. Finally, expressing u n+1 i we obtain the required inequality. The inequality on the left-hand side of (4.5) can be proved similarly. This completes the proof. 5. Qualitative Properties of the Numerical Solutions of Parabolic Problems. We apply the results, which were formulated for arbitrary PO-iterations, of the previous section for the fully discretized numerical model (3.9). Based on Theorem 3.2, the finite difference and Galerkin finite element methods can be written in the form of a PO-iteration with the choice h n = ∆tf (n,θ) . This is why the qualitative properties (DMP, DNP and DMNC) of the numerical methods can be defined in the same manner like for the PO-iterations.  Proof. The first part of the statement follows from Theorem 4.2, Lemma 4.6 and Theorem 5.1. The second one is a consequence of Lemma 4.8.

Corollary 5.3. A finite difference or a Galerkin finite element method for the problem (2.4)-(2.5) is qualitatively adequate if and only if it satisfies properties
, that is if the method preserves the nonnegativity. Remark 5.4. In view of the results of [13], there exist parameter choices when the finite difference or the finite element method is contractive in maximum norm, but it does not preserves the nonnegativity.
Remark 5.5. Let us notice that our results can be applied not only for the finite difference and the Galerkin finite element method but also for any numerical method that can be written in the form of a PO-iteration.
6. Sufficient a Priori Conditions of the Qualitative Properties. The necessary and sufficient conditions in the previous sections can not be checked easily, without computing the elements of the matrices. More precisely, the theorems do not state anything explicitly about the choice of the mesh and the choice of the parameters θ and ∆t (which define the elements of the matrices M and K) in order to guarantee the qualitative properties. Our aim is to find conditions that can be checked a priori and imply the DNP / DMP. Then, these conditions, according to Theorem 5.2, imply the qualitative adequateness of the discrete models. In this section, we will give some a priori conditions. In order to present a complete work for researchers involved in scientific computing, we also added two corresponding results from the literature. Proof. We have to show that under the assumptions of the theorem the properties (P 1) − (P 2) hold. For the finite difference method A −1 0 ≥ 0 and conditions (S1) and (S2) are fulfilled automatically. This implies property (P 2). Condition (S3) implies that the main diagonal of B is nonnegative. The off-diagonal of B is also nonnegative because the off-diagonal of K is non-positive. Hence, by virtue of A −1 0 ≥ 0, we obtain property (P 1).
In the case of the finite element method, relations (S1) and (S3) yield B ≥ 0. Condition A ∂ ≤ 0 follows from (S2). A 0 is a Gram-matrix, and because of (S2), the off-diagonal elements of A 0 are non-positive. This implies that A 0 is an M -matrix (see [1, p. 134]), thus A −1 0 ≥ 0. These facts yield (P 1) and (P 2), thus the DMP. In the expressions of this section, for the sake of simplicity, fractions with zero denominators are understood as infinity.
The condition guarantees the validity of the above relations, where h min is the minimal spatial step size in the mesh and κ max = sup x∈Ω {κ(x)}. In 2D, similarly, we obtain the sufficient condition The above results can be summarized as follows. We remark that for 1D problems on uniform meshes of step-size h, with constant material parameters (C 0 , κ 0 ), the necessary and sufficient condition of the DMP is 6.2. Conditions for the finite element method. For finite element methods, condition (S1) is generally guaranteed by some geometrical requirements for the mesh. Moreover, conditions (S2) and (S3) can be achieved by some lower and upper bounds for the time-step.

Finite element method in 1D.
Computing the elements of the stiffness matrix, we can check that the off-diagonal elements are non-positive, thus condition (S1) is satisfied. Computing additionally the elements of the mass matrix we obtain that (S2) and (S3) are valid if Thus we obtain the next theorem.
. Let us notice that unlike in the case of finite difference methods, the above theorem yields a lower bound for the time-step. This is a general phenomenon for finite element methods. This lower bound exists not only in sufficient conditions but in necessary and sufficient conditions, too. This is shown by the next theorem obtained in [6] for the DNP.

The Galerkin finite element solution of the one-dimensional problem (2.4)-(2.5) with piecewise linear elements on uniform grid and with constant material parameters satisfies the DMP if and only if the relations
hold.

Finite element method in 2D.
In 2D we apply the hybrid mesh in the spatial discretization. On triangular elements linear and on rectangular elements bilinear basis functions are used.
Let us consider a hybrid mesh T h . Let us define the number for each triangle T of the mesh, where α 1 , α 2 and α 3 denote the angles of the triangle. Let us define σ = min T ∈T h σ T . Furthermore, let us introduce the value for each rectangle R, where a and b denote the edges of the rectangle. Let us define µ = min R∈T h µ R . The hybrid mesh T h is called to be of strictly compact type if σ > 0 and µ > 0. The condition σ > 0 means that each angle of the triangles in the mesh is of acute type, and µ > 0 means that the longer edges of the rectangles are not greater than √ 2 times the shorter ones. The main result of [8] is formulated as follows.
Theorem 6.6. The Galerkin finite element solution of the two-dimensional problem (2.4)-(2.5) using linear/bilinear basis functions on a hybrid mesh of strictly compact type satisfies the DMP if The symbols area(T ) and area(R) denote the area of the triangle T and the rectangle R, respectively. l max is the length of the longest edge in T . Thus, for strictly compact meshes, the off-diagonal elements of the stiffness matrix are non-positive, hence condition (S1) is fulfilled. Relation (6.1) implies the conditions (S2) and (S3).
Remark 6.7. For purely rectangular meshes, we have the weaker lower bound for ∆t in the form For a square mesh with step size h and with constant material parameters a sufficient condition of the DMP is (In case of θ = 1 the upper bound disappears.) This shows that in this case the DMP can be guaranteed only (with our sufficient condition) for methods with θ ≥ 2/3.

Finite element method in 3D.
In 3D we investigate two different meshes. The first one is the rectangular and the second one is the tetrahedral mesh. Let us start with the first one, where trilinear basis functions are applied. With simple but tedious calculations we can compute the elements of the mass and stiffness matrices, and we can notice that the condition (S1) can be valid only for uniform meshes. Moreover, we can observe that condition (S2) cannot be guaranteed, because there are also positive off-diagonal elements in the matrix A. Thus A 0 is not an Mmatrix in this case. We have to guarantee the nonnegative matrix inverse by means of other tools. In work [12], the so-called Lorenz-criterion ( [21]) was applied.

Theorem 6.8. The Galerkin finite element solution of the three-dimensional problem (2.4)-(2.5) using trilinear basis functions on a uniform rectangular mesh with constant material parameters satisfies the DMP if
For tetrahedral meshes the DMP was discussed in paper [11]. We will recall this result. Let us consider a tetrahedral mesh T h of the solution domain Ω, and define the number σ T for each tetrahedron T of T h as σ T = min{cos β 1 , . . . , cos β 6 }, where β 1 , . . . , β 6 denote the angles made by any two faces of the tetrahedron T . If each angle is less than π/2, then σ T > 0. Let us set σ = min T ∈T h σ T . If σ > 0, then the mesh is called strictly acute type.
Theorem 6.9. The Galerkin finite element solution of the three-dimensional problem (2.4

)-(2.5) with constant material parameters using linear basis functions on a strictly acute type tetrahedral mesh satisfies the DMP if the condition
holds, where p min and p max denote the minimal and maximal perpendicular length in the mesh.
7. Numerical Examples. In this section, we demonstrate the validity of our results on some numerical examples in one, two and three dimensions. For the sake of simplicity, we suppose that the problem (2.4)-(2.5) describes a heat conduction process (see Remark 2.7). The material parameters are set to be constant one. Furthermore, we suppose that there are no heat sources or sinks present in the computational domain.

Examples in 1D.
In the first example, the 1D heat equation is solved with the Crank-Nicolson (θ = 1/2) finite difference method on the interval [0,1]. We are interested in the approximation of the temperature at the time-level t = 1. We apply an equidistant mesh with the spatial step-size h = 1/30. The temperatures at the left and right boundaries are given by the functions µ 1 (t) = 1 and µ 2 (t) = 0, respectively. A nonnegative approximation of a continuous and nonnegative initial function can be seen in Figure 7.1. Using this initial grid-function, the approximations of the temperature at the time-level t = 1 can be seen in Figure 7.2. The figures on the left-hand and right-hand sides were obtained using the time-steps ∆t = 1/11 and ∆t = 1/31, respectively. These approximations show that neither the nonnegativity preservation nor the maximum-minimum principle is valid for the above chosen timesteps. In order to get a qualitatively correct approximation, ∆t has to satisfy the condition (c.f. equation (6.3)). This bound shows that ∆t was chosen to be too large in the previous numerical example. Naturally, the above upper bound is the necessary and One can think that choosing the time-step to be sufficiently small we can obtain a qualitatively adequate solution. The next example shows that this is generally not the case. Namely, as we have shown, for finite element methods not only upper but lower bounds do exist. Let us solve the 1D heat equation with the Galerkin finite element method on the interval [0,1] using the Crank-Nicolson (θ = 1/2) scheme in the time integration. Let us suppose that the initial function is a nonnegative continuous function and it takes on values only from the interval [0, 1]. The discretization of the initial function is carried out on an equidistant mesh with h = 1/10 (Figure 7.4). The temperatures at the left and right boundaries are given by the functions µ 1 (t) = 1 and µ 2 (t) = 0, respectively. The numerical solutions calculated at the time-levels t = 1/1000 and t = 1 with the time steps ∆t = 1/100000 and ∆t = 1/9 can be seen, respectively, on the left-hand and right-hand sides of Figure 7.5. Both numerical solutions are qualitatively incorrect. Namely, the maximum-minimum principle is broken because the solutions have nonnegative values and values that are greater than one. In order to obtain qualitatively adequate solutions we can apply Theorem 6.4, which states that the maximum-minimum principle is valid providing that ∆t is chosen according to the relation 1/300 = 0.0033 ≤ ∆t ≤ 0.0067 = 2/300. In Figure 7.6, the numerical solution is calculated at the time-level t = 1 with the time-step ∆t = 0.005, which falls into the above interval. As it was expected, the solution has values only from the interval [0, 1].
In the third example we solve again the 1D heat equation on the interval [0,1]. The temperature at the right-hand side is equal to one (µ 2 (t) = 1), while the temperature In the case of the numerical solutions shown in Figure 7.7, the maximum-minimum principle is broken. The time-step was chosen to be too large (c.f. (6.3)). We can observe that the numerical solution will be periodic after the transient period (t > 3) (the boundary condition is periodic and the effect of the initial function disappears). The approximations shown in Figure 7.7 appear regularly in every third time-steps. Thus the numerical example demonstrates that it is possible that the qualitatively bad behavior of the numerical solution does not disappear while increasing the time-level at which the numerical solution is computed.

Examples in 2D.
We solve the two-dimensional heat equation on the unit square [0, 1]×[0, 1] with homogeneous boundary condition. We apply the finite element method with bilinear elements on a rectangular mesh with ∆x = 1/10 and ∆y = 1/12. In the time-discretization the θ-method is used with θ = 0.9 and with a fixed time-step ∆t. A nonnegative discretization of a nonnegative and continuous initial function can be seen on the left-hand side of Figure 7 remark that, with ∆t = 0.0005, the numerical solution has negative components at the first 11 time-levels, while with ∆t = 0.5, the solution has negative components at every time-level. Choosing the time-step to be ∆t = 0.01, we obtain an adequate approximation (right-hand side of Figure 7.8).
We can calculate the necessary and sufficient condition of the DNP using Theorem 4.2, which results in the condition 0.0050 ≤ ∆t ≤ 0.0155. This shows that the chosen time-steps were too small or too large. A sufficient condition obtained applying the results of Theorem 6.6 is 0.0066 ≤ ∆t ≤ 0.0136.
In the second example, we apply the finite difference method with θ = 0.9 for the initial approximation seen on the left-hand side of Figure 7 numerical approximation, it is enough to choose ∆t to be not greater than 0.017 (Theorem 6.2).

Example in 3D.
In order to give a 3D example, we solve the homogeneous heat equation on the unit cube. We apply the finite difference method with the Crank-Nicolson time-integration. The spatial discretization is performed with an equidistant mesh, where ∆x = ∆y = ∆z = 1/10. The approximation of a continuous nonnegative initial function is zero in every grid point excepting 27 grid points in the middle of the region, where the temperature is approximated by 50. The approximations of the temperature at the point (5/10, 4/10, 4/10) using the time-steps ∆t = 0.05 and ∆t = 0.01 can be seen in Figure 7 temperatures, which contradicts to the nonnegativity preservation property. Theorem 6.2 gives a sufficient condition of the nonnegativity, namely ∆t must be not greater than 1/300.

Conclusions.
In this paper, we have shad light on the connections between the main qualitative properties of parabolic partial differential equations and their numerical solution methods. We have found that the maximum-minimum principle, the nonnegativity preservation and the maximum norm contractivity are equivalent properties for the continuous problem, but this is not the case when we consider the finite difference or Galerkin finite element solutions of the problem. In this case, the maximum-minimum principle is equivalent to the nonnegativity preservation, and the maximum norm contractivity yields a weaker condition as this property is implied by the properties mentioned earlier. Thus, to achieve a qualitatively adequate method, we have to apply a nonnegativity preserving method, which can be guaranteed by the sufficient conditions of Theorem 6.1 and by its special cases listed in §6.