- Email: [email protected]
Maximum norm contractivity in the numerical solution of the one-dimensional heat equation
Applied Numerical Mathematics 31 (1999) 451–462
Maximum norm contractivity in the numerical solution of the one-dimensional heat equation Róbert Horváth ∗ University of Sopron, Institute of Mathematics, Ady E. u. 5, Sopron, H-9400, Hungary
Abstract In this paper we consider the one-dimensional heat conduction equation (Friedmann, 1964). To the numerical solution of the problem we apply the so-called (σ, θ )-method (Faragó, 1996; Thomée, 1990) which unites a few numerical methods. With the choice σ = 0 we get the finite difference θ -method and the choice σ = 16 results in the finite element method with linear elements. The most important question is the choice of the suitable mesh-parameters. The basic condition arises from the condition of the convergence (Faragó, 1996; Samarskii, 1977; Thomée, 1990). Further conditions can be obtained aiming at preserving some qualitative properties of the continuous problem. Some of them are the following: non-negativity, convexity, concavity, shape preservation and contractivity in some norms (Dekker and Verwer, 1984). Now we shall study the maximum norm contractivity. There are some results in the literature for the parameter choices which guarantee this property (Kraaijevanger, 1992; Samarskii, 1977; Thomée, 1990). However these papers specialize only on the finite difference methods and give sufficient conditions. We determine the necessary and sufficient conditions related to the (σ, θ )-method. We close the paper with numerical examples. 1999 Elsevier Science B.V. and IMACS. All rights reserved. Keywords: Numerical solution; Qualitative properties; Heat equation; Maximum norm contractivity
1. Introduction Let us consider the following parabolic partial differential equation [8], the so-called one-dimensional heat conduction equation ∂u ∂ 2 u = 2, ∂t ∂ξ u(t, 0) = u(t, 1) = 0,
t > 0,
u(0, ξ ) = u0 (ξ ),
ξ ∈ [0, 1],
ξ ∈ (0, 1),
t > 0,
∗ E-mail: [email protected].
0168-9274/99/$20.00 1999 Elsevier Science B.V. and IMACS. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 9 9 ) 0 0 0 0 7 - 0
(1)
452
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
in the domain Ω = (0, ∞) × (0, 1), where u0 is a suitable smooth function for which there exists a unique solution on Ω. The function u : Ω → R is called the solution of the problem (1) if u ∈ C 1,2 (Ω) ∩ C(Ω) and it satisfies the equations in (1). This problem has some characteristic qualitative properties as the physical model [3,11,14]. Applying the well-known maximum principle [14] we can show that the function t 7→ maxx∈[0,1] |u(t, x)| is monotonically decreasing, that is the absolute value of the temperature does not increase in time. This property is called maximum norm contractivity. It can be shown that the function t 7→ u(t, x) tends to 0 for all x ∈ [0, 1] if t → ∞. (If u0 is a concave or a convex function then the convergence is monotone with respect to x.) It is said that the solution composes itself. Moreover, we say that the problem (1) preserves the non-negativity if u0 > 0 implies that u > 0 in Ω. Among the above properties we shall deal with the maximum norm contractivity. For a numerical solution to the problem (1) we define the following mesh for fixed τ ∈ R+ and h = 1/(n + 1) (n ∈ N) [13,15]:
Ωτ,h := (tj , xi ) ∈ Ω | xi = ih (i = 0, . . . , n + 1), tj = j τ (j ∈ N) . Denoting the approximation to the exact solution at the j th time level by y j ∈ Rn and using the so-called (σ, θ)-method [7] we define the following one-step iteration process: y j +1 − y j 1 1 (2) = −θ 2 Qy j +1 − (1 − θ) 2 Qy j , j = 0, 1, . . . . τ h h Here Q = tridiag[−1, 2, −1] ∈ Rn×n and M = I − σ Q ∈ Rn×n are matrices (the matrix I denotes the unit matrix), θ ∈ [0, 1] and σ ∈ [0, 14 ) are arbitrary parameters and the vector y 0 is a suitable approximation of the initial function u0 . This method unites a few numerical methods. If σ = 0 then we get the finite difference method and in the case of σ = 16 the method (2) results in the finite element method with linear elements. First of all we reformulate our problem to a problem of linear algebra. Let us introduce the notations q = τ/ h2 , z = θq − σ . If z = 0 then the method (2) is called explicit. Let us denote the matrix I − qQ by X. If z 6= 0 then we call the method implicit and introduce the notations s = q − z, x = 2 + 1/z, X 1 = Q + (x − 2)I and X 2 = I − sQ. We notice that s > 0 as s = q(1 − θ) + σ and |x| > 2 because z > − 14 . In this case X denotes the matrix (1/z)X−1 1 X 2 . The matrix X is well-defined because the matrix X 1 is invertible for all values |x| > 2 [12]. Thus the iteration (2) can be rewritten in the form M
y j +1 = Xy j ,
j = 0, 1, . . . .
(3)
It is an interesting question to examine the behaviour of the numerical solution defined by the iteration (3). We seek conditions for both the method and the mesh under which the qualitative properties of the continuous problem are preserved in the numerical solution [3–7,9,10]. We require the iteration (3) to be composing, too. This means that the vector-sequence y j (j = 0, 1, . . .) has to tend to the zero vector for all initial vectors y 0 , that is the matrix X must be a so-called convergent matrix. In this case we call the method (2) weak convergent method. We get the following condition for the weak convergence. Lemma 1 [7]. Assume that σ ∈ [0, 14 ) and θ ∈ [0, 0.5). Then the (σ, θ)-method is weak convergent for fixed n if and only if the condition 1 1 − 4σ cos2 (π/(2(n + 1))) 2(1 − 2θ) cos2 (π/(2(n + 1))) is satisfied. Moreover, it is weak convergent for all n if and only if the condition 1 − 4σ q6 2(1 − 2θ) q6
(4)
(5)
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
453
holds. For the parameters σ ∈ [0, 14 ) and θ ∈ [0.5, 1] the method is always weak convergent. Naturally in this paper we consider schemes satisfying (4) and we analyze the properties of the discrete maximum norm contractivity. Definition 2. We call the numerical process (2) contractive in the maximum norm if the inequalities
j +1
y
∞
6 y j ∞
are fulfilled for each j ∈ N. We note that the maximum norm contractivity is not required for neither the stability nor the convergence of the method, it is desirable only from qualitative point of view. Obviously the method (3) is contractive in the maximum norm if and only if the maximum norm of the matrix X is not greater than 1, that is if the condition kXk∞ 6 1 is fulfilled. Our aim is to give conditions for the values of the mesh-parameters by a fixed (σ, θ)-method under which the numerical solution is contractive in maximum norm. 2. The condition of the contractivity In this section we give the condition of the maximum norm contractivity both by fixed n and for all values n. Let us consider the case of z = 0. It follows from the form of the matrix X, that if n = 1 or n = 2, then the necessary and sufficient condition of the contractivity is q 6 1 or q 6 23 , respectively. If n > 3 then the condition is q 6 12 . Now let us pass over the implicit methods. Let G ∈ Rn×n be a symmetric, one-pair matrix defined as follows: (G)i,j = (Γ )i,j if x > 2 and (G)i,j = (−1)i+j −1 (Γ )i,j if x < −2, where sh(iα) sh((n + 1 − j )α) |x| γi,j , if i 6 j , (Γ )i,j = γi,j = , α = arch . (6) γj,i , if i > j , sh α sh((n + 1)α) 2 For our purpose we rewrite the matrix X to an explicit form. The lemma in question can be stated as follows. Lemma 3. Assume that z 6= 0. For the matrix X defined by (1/z)X−1 1 X 2 we have the relation
1 q X= G − sI . z z
(7)
Proof. The matrix G is the inverse of the matrix Q + (x − 2)I [12]. Multiplying the matrix X by the matrix z · (Q + (x − 2)I ) from the left we obtain the matrix I − sQ. 2 To find the maximum norm of the matrix X we have to examine the structure of the matrix and to sum the absolute values of the elements in each row. Using the notation ν = [(n + 1)/2] (the integer part of (n + 1)/2) we can state the following: Lemma 4. For the diagonal elements of the matrix G are valid the following relations: γi,i = γn+1−i,n+1−i ,
i = 1, . . . , ν,
γ1,1 < γ2,2 < · · · < γν,ν .
(8)
454
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
Proof. To prove this lemma it is sufficient to analyze the continuous function ω 7→ sh(ωα) sh((n + 1 − ω)α) on the interval [1, n]. 2 Introducing the notation Gi =
Pn
j =1 |(G)i,j |
=
Pn
j =1 (Γ )i,j
we have
Lemma 5. The sum of the modulus of the elements in the ith row of the matrix G can be expressed as follows: 1 Gi = (9) (1 − γ1,i − γi,n ), i = 1, . . . , n. |x| − 2 Proof. Using the notation bz = 1/(|x| − 2) we have the identity
T
Q + (|x| − 2)I · bz, . . . , bz
T
= (1, . . . , 1)T + bz, 0, . . . , 0, bz .
Multiplying (10) by the matrix Γ we obtain the statement.
(10)
2
Because of Lemma 4 for the sign of z and for the sign distribution of diagonal elements of the matrix X the possible cases are the following: Case I. z > 0 and (X)i,i > 0, i = 1, . . . , n. Case II. z > 0 and (X)i,i < 0, i = 1, . . . , n. Case III. z > 0 and there exists k ∈ {1, . . . , ν −1} such that (X)1,1 , . . . , (X)k,k < 0 and (X)k+1,k+1 , . . . , (X)ν,ν > 0. Case IV. z < 0 and (X)i,i < 0, i = 1, . . . , n. Case V. z < 0 and (X)i,i > 0, i = 1, . . . , n. Case VI. z < 0 and there exists l ∈ {1, . . . , ν − 1} such that (X)1,1 , . . . , (X)l,l > 0 and (X)l+1,l+1 , . . . , (X)ν,ν < 0. For the different cases we examine the condition of the maximum norm contractivity. Our results can be summarized as follows. Lemma 6. Assume that the method (2) is weak convergent. Then in the case V the maximum norm contractivity is not fulfilled. In cases I, III, IV and VI the contractivity is valid without further restrictions for the parameters. Finally, in the case II we get the inequality
sh((n + 1 − ν)α) sh(να) sh((n + 1 − ν)α) + sh(να) + 2(|x|−2) q s sh α 1− + 61 z sh((n + 1)α) z by fixed n and the condition 8σ (θ − 1) + 2 − θ +
p
(8σ (θ − 1) + 2 − θ)2 − 16(1 − θ)2 σ (4σ − 1) 8(1 − θ)2 for all n ∈ N which guarantee the contractivity. In case of θ = 1 no restrictions arise. q6
(11)
(12)
Proof. Because of the simple but long calculations of the proof, the steps are not detailed here. For the detailed proof we refer to Appendix. 2 It is remains to examine the six cases above, by which conditions for the parameters are realized. Comparing the bounds, we can give a necessary and sufficient condition of the maximum norm
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
455
contractivity. We are not able to give the conditions by fixed n because most of our inequalities are implicit. To see some numerical results by fixed n, we refer to Section 3. We can formulate the most important result of our paper as follows: Theorem. Let θ ∈ (0, 1) and σ ∈ [0, 14 ) be fixed. The method (2) is both weak convergent and contractive in maximum norm for all values of n if and only if θ > 2σ and the parameter q satisfies the following condition: σ 8σ (θ − 1) + 2 − θ + 6q 6 θ
p
(8σ (1 − θ) − 2 + θ)2 − 16(1 − θ)2 σ (4σ − 1) . 8(1 − θ)2
(13)
In case of θ = 1 the condition is q > σ and in the case θ = 0 (consequently σ = 0) the condition is 0 6 q 6 12 . Proof. The proof is based on Lemma 6. First let us consider the case z < 0. According to Lemma 6 the weak convergent method (2) is contractive if and only if there exists a negative element in the main diagonal of the matrix X (in cases IV and VI). From this condition we get the inequality (X)ν,ν < 0. This is fulfilled for all n if it is fulfilled for n = 1. It yields the condition q > (1 − 2σ )/(2(1 − θ)). We have to compare this bound with the bounds z < 0 and (5). From the first comparison we get the necessary condition θ < 2σ . On the other hand, in this case 1 − 2σ 1 − 4σ > 2(1 − θ) 2(1 − 2θ) so the conditions 1 − 2σ q> 2(1 − θ)
and
q6
1 − 4σ 2(1 − 2θ)
cannot be fulfilled simultaneously. Thus in the case z < 0 there is no parameter choice σ, θ, q which yields both weak convergence and maximum norm contractivity for all n. Now let us examine the case z > 0. Because of the conditions (5) and z > 0 we get the necessary condition θ > 2σ . In cases I and III the contractivity is fulfilled automatically. These cases are realized if and only if (X)ν,ν > 0. This is fulfilled for all n if it is fulfilled for n = 1. It yields the condition q 6 (1 − 2σ )/(2(1 − θ)). Comparing this condition with the condition θ > 2σ we obtain the inequality σ/θ < q 6 (1 − 2σ )/(2(1 − θ)). (If θ = 1 then q > σ .) If q > (1 − 2σ )/(2(1 − θ)) then there exists a natural number n0 that for all n < n0 the elements in the main diagonal of the matrix X are negative. But in this case we got the following condition for q in Lemma 6: q6
8σ (θ − 1) + 2 − θ +
p
(8σ (1 − θ) − 2 + θ)2 − 16(1 − θ)2 σ (4σ − 1) . 8(1 − θ)2
Thus in the case of z > 0, the conditions of the contractivity are θ > 2σ and σ 8σ (θ − 1) + 2 − θ +
p
(8σ (1 − θ) − 2 + θ)2 − 16(1 − θ)2 σ (4σ − 1) . 8(1 − θ)2
(If θ = 1 then q > σ .) Let us take into consideration the explicit methods as well. In this case we have got the condition q 6 12 . Thus the theorem is proved. 2
456
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
To analyze the conditions by fixed σ we introduce the following functions: lσ (θ) = cσ (θ) =
1 − 4σ 2(1 − 2θ)
θ<
1 2
pσ (θ) =
,
8σ (θ − 1) + 2 − θ +
p
1 − 2σ 2(1 − θ)
(θ < 1),
eσ (θ) =
σ θ
(8σ (1 − θ) − 2 + θ)2 − 16(1 − θ)2 σ (4σ − 1) 8(1 − θ)2
(θ > 0), (θ < 1).
They denote the barrier functions of the convergence. By the other values of θ let us define the values of the functions as ∞. Then the following statements are true: (a) Let σ 6= 0. Then lσ (2σ ) = cσ (2σ ) = pσ (2σ ) = eσ (2σ ) = 12 , if θ > 2σ then lσ (θ) > cσ (θ) > pσ (θ) > eσ (θ) and if θ 6 2σ then lσ (θ) 6 cσ (θ) 6 pσ (θ) 6 eσ (θ). (b) Let σ = 0. Then l0 (θ) 6 c0 (θ) 6 p0 (θ). Remark. Let us apply the theorem for two well-known methods. In the case of σ = 0 we have the finite difference method. Then the method (2) is contractive in the maximum norm for all n if and only if 06q 6
2−θ . 4(1 − θ)2
(14)
The choice θ = 0 corresponds to the forward Euler method (in this case the condition is q 6 0.5), θ = 0.5 corresponds to the Crank–Nicolson method (in this case the condition is q 6 1.5) and the choice θ = 1 corresponds to the backward Euler method (in this case q is optional). We have got the upper bound which was obtained by Kraaijevanger [10] with other tools. The method of that reference does not make possible the examination of the contractivity bound depending on n. Our method as it is shown in Section 3 does carry out this.
Fig. 1.
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
457
The other special choice is σ = 16 which is the finite element method with linear elements. Then the method (2) is contractive in maximum norm for all n if and only if θ > 13 and the condition √ 1 θ + 2 + 9θ 2 − 12θ + 12 (15) 6q 6 6θ 24(1 − θ)2 holds. (If θ = 1 then q > 16 .) The dark area shows the possible parameter choices in Fig. 1. Remark. According to the continuous case, the method (2) is called non-negativity preserving, if the vectors y j (j = 1, 2, . . .) are non-negative for all non-negative initial vectors y 0 . Considering the iteration (3) we can see that the necessary and sufficient condition of the non-negativity preservation is the inequality X > 0. This condition is fulfilled if and only if z > 0 and (X)i,i > 0 (i = 1, 2, . . . , n). It is apparent from the proof of Lemma 6 that the non-negativity preservation is a sufficient condition of the maximum norm contractivity. 3. Numerical examples Let us assume that (σ, θ) are fixed. To define the mesh on which this method is contractive in the maximum norm by fixed n, we have to solve the implicit inequality (11). In our examples we execute it numerically. The necessary and sufficient condition of both the maximum norm contractivity and the weak convergence by fixed n, θ, σ is the following: (a) Assume that σ 6= 0. If θ > 2σ then cσ,n (θ) > q > eσ (θ). Conversely, if θ < 2σ then eσ (θ) > q > pσ,n (θ) and q 6 lσ,n (θ). (b) In the case σ = 0 we have the condition 0 6 q 6 c0,n (θ). Here cσ,n (θ) denotes such a value of q for which the relation (11) turns into the equality. On the base of the proof of Lemma 6 we can verify that by n → ∞ cσ,n (θ) → cσ (θ) and cσ,n (θ) forms a monotonically decreasing sequence in n. Thus we can give for all fixed n a greater bound for q than cσ (θ) defined in [10]. Table 1 The case of σ = 0 n
θ = 0.00
θ = 0.30
θ = 0.45
θ = 0.50
1
1
2.5
10
∞
3
0.5
1.0763
2.9894
990237
5
0.5
0.8961
1.5035
2.0000
7
0.5
0.8721
1.3334
1.6180
9
0.5
0.8682
1.2952
1.5352
11
0.5
0.8675
1.2850
1.5113
13
0.5
0.8674
1.2822
1.5037
15
0.5
0.8674
1.2813
1.5012
25
0.5
0.8673
1.2810
1.5000
c0 (θ )
0.5
0.8673
1.2810
1.5000
458
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
Table 2 σ = 16 , θ = 0.1 n
p1/6,n (0.1)
l1/6,n (0.1)
1
0.3704
0.8333
2
0.2941
0.4167
3
0.2810
0.3156
4
0.2784
0.2743
5
0.2779
0.2533
7
0.2778
0.2331
9
0.2778
0.2240
∀n
c1/6 (0.1) = 0.2778
l1/6 (0.1) = 0.2083
Table 3 σ = 16 , θ = 0.5 n
3
5
7
9
11
15
25
∀n
c1/6,n (0.5)
2.3333
0.9658
0.9034
0.8973
0.8957
0.8954
0.8954
0.8954
Moreover, lσ,n (θ) denotes the bound derived from (4) and pσ,n (θ) is the solution of the equality X ν,ν = 0 for q. It can be verified that pσ,1 (θ) = pσ (θ) and by n → ∞ both pσ,n (θ) → cσ (θ) and lσ,n (θ) → lσ (θ) monotonically decreasing. Because of pσ (θ) > lσ (θ) if θ < 2σ there does not exist q which yields both contractivity and weak convergence for all n. First of all let us consider the case σ = 0. Table 1 shows the computed values c0,n (θ) in the cases θ = 0, 0.3, 0.45, 0.50. We remark that these values satisfy the condition of the weak convergence as well. Now let us consider the case σ = 16 . We must consider two different cases according to θ 6 2σ or θ > 2σ . Let θ be for example 0.1 < 2σ = 23 . Then by fixed n the conditions of the maximum norm contractivity and the weak convergence are satisfied if and only if 53 > q > p1/6,n (0.1) and q 6 l1/6,n (0.1). Table 2 shows that weak convergence and contractivity are fulfilled only in the cases n = 1, 2, 3. Now let θ be 0.5 > 13 . Then by fixed n the conditions of the maximum norm contractivity and weak convergence are satisfied if and only if 13 6 q 6 c1/6,n(0.5). Table 3 shows the values c1/6,n (0.5) computed with numerical methods. Appendix Proof of Lemma 6. Let us denote by X i (i = 1, . . . , n) the sum of the elements in absolute value in the ith row of the matrix X. The values Xi using Lemma 3 can be expressed in the following form: Xi =
n n q q X q X s (Γ ) + (X) = (Γ ) + (G) − , i,j i,i i,j i,i z2 z2 j =1 z2 j =1 z j 6=i
j 6=i
i = 1, . . . , n.
(A.1)
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
459
Our aim is to compute the value of kXk∞ and to give the condition kXk∞ 6 1. We have established in Lemma 6 that there exist only six sign-distributions in the matrix X. We examine these cases. We apply the expression (A.1) and we use Lemmas 3 and 5 in the calculations. Case I (z > 0 and (X)i,i > 0, i = 1, . . . , n). In this case we get the the following estimation: Xi =
n n q X q s s q X (Γ ) + (Γ ) − (Γ )i,j − = i,j i,i 2 2 2 z j =1 z z z j =1 z j 6=i
q 1 s q s q −s z = 2 (1 − γ1,i − γi,n ) − = (1 − γ1,i − γi,n ) − 6 = = 1. z |x| − 2 z z z z z The inequality above shows that the condition z > 0 and (X)i,i > 0 (i = 1, . . . , n) is a sufficient condition for the maximum norm contractivity. Case II (z > 0 and (X)i,i < 0, i = 1, . . . , n). In this case we have Xi =
n n q X q s 2q s q X (Γ ) − (Γ ) + (Γ )i,j − 2 (Γ )i,i + = i,j i,i 2 2 2 z j =1 z z z j =1 z z j 6=i
q 1 f (i) s = 2 1− + , z |x| − 2 sh((n + 1)α) z where 2(|x| − 2) sh (n + 1 − i)α sh(iα), Df = {1, 2, . . . , n}. sh α Extending the function f on the set R we get a differentiable function. With differentiation we can easily obtain that this function has a local minimum point at i = (n + 1)/2 and it has local maximum points at
f (i) = sh (n + 1 − i)α + sh(iα) +
n+1 1 1 n+1 i= ± arch sh α sh α . 2 α 2 2 We verify that the value X i (i = 1, . . . , n) takes its maximum value by the choice i = [(n + 1)/2] = ν. Thus we have to show that f (1) > f (ν). We prove this statement for odd values of n. In this case ν = (n + 1)/2 and the inequality can be rewritten in the following form:
2(|x| − 2) n+1 2(|x| − 2) 2 n + 1 sh α + sh(nα) + sh(nα) sh α > 2 sh α + sh α . sh α 2 sh α 2 Using the identities of hyperbolic functions we get the inequality
n+1 n−1 n+1 n+1 (|x| − 2) sh sh2 α ch α − sh α > α sh α sh(nα) . 2 2 2 sh α 2 This is fulfilled by the choice n = 1. Let us analyze both functions standing on the right and on the left hand sides, respectively. As earlier we assume that they are functions on R. Since
n+1 (|x| − 2) sh (n + 1)α − 2 sh α ch(nα) (A.2) α > 2 sh α therefore the derivative of the left hand side is greater than the derivative of the right hand side. To prove (A.2) it is sufficient to use the well-known identities for the hyperbolic functions.
ch(nα) − ch
460
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
We have proved that in case II, the condition of the maximum norm contractivity can be defined as follows: q f (ν) s kXk∞ = 1− + 6 1. (A.3) z sh((n + 1)α) z We seek the condition which guarantees the contractivity for all n. Let n be an odd natural number. Then the condition (A.3) yields q 1 n+1 1 s |x| − 2 1− + 6 1. th α − n+1 2 z |x| − 2 sh α 2 z ch( 2 α) This is satisfied for all n if and only if q 1 |x| − 2 s 1− + 6 1. 2 z |x| − 2 sh α z Substituting the values of the parameters we get the condition
4q 2 (1 − θ)2 + q 8σ (1 − θ) − 2 + θ + 4σ 2 − σ 6 0 from which we obtain the condition p 8σ (θ − 1) + 2 − θ + (8σ (1 − θ) − 2 + θ)2 − 16(1 − θ)2 σ (4σ − 1) q6 . (A.4) 8(1 − θ)2 If θ = 1 then no condition arises. Case III (z > 0 and there exists k ∈ {1, . . . , ν − 1} such that (X)1,1 , . . . , (X)k,k < 0 and (X)k+1,k+1 , . . . , (X)ν,ν > 0). In this case we consider the following two functions on the set {1, 2, . . . , n}: n q X s g1 (i) = 2 (Γ )i,j − , z j =1 z
1 q f (i) s 1− + . g2 (i) = 2 z |x| − 2 sh((n + 1)α) z
The function g1 (i) results in the sum of the ith row of the matrix X and g2 (i) gives the sum of the ith row, taking the element in the main diagonal with opposite sign. Thus if i 6 k then g1 (i) < g2 (i) = Xi and if k < i 6 ν then X i = g1 (i) > g2 (i). As we know from Lemma 4 and case II these functions take their maximum values at the point i = ν. So g1 (ν) > g1 (i) and g1 (ν) > g2 (i) (i = 1, . . . , n). On the other hand, g1 (ν) 6 1, that is the condition kXk∞ 6 1 is fulfilled in this case as well. Case IV (z < 0 and (X)i,i < 0, i = 1, . . . , n). In this case we can calculate as follows: n n q X s q X s q Xi = 2 (Γ )i,j + − 2 (G)i,i = 2 (Γ )i,j + z j =1 z z z j =1 z j 6=i
q 1 sh((n + 1 − i)α) + sh(iα) s = 2 1− + z |x| − 2 sh((n + 1)α) z q 1 s 4s − 1 2(2(s − z) − 1) 6 2 + = =1+ . z |x| − 2 z 4z + 1 4z + 1 Because of z > − 14 we have to prove the relation 2(s − z) 6 1. Obviously this inequality is fulfilled in the case θ > 12 . In case of θ < 12 we get the upper bound q6
1 − 4σ . 2(1 − 2θ)
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
461
This condition is satisfied automatically if the matrix X is convergent. Case V (z < 0 and (X)i,i > 0, i = 1, . . . , n). In this case the expression X i can be reformulated as follows: Xi =
n n q X q s 2q s q X (Γ ) − (Γ ) − (Γ )i,j − 2 (Γ )i,i − = i,j i,i 2 2 2 z j =1 z z z j =1 z z j 6=i
sh((n + 1 − i)α) sh(iα) sh((n + 1 − i)α) + sh(iα) + 2(|x|−2) q 1 s sh α = 2 1− − . z |x| − 2 sh((n + 1)α) z X i takes its maximum value at i = ν on the set {1, . . . , n}, that is kXk∞ = X ν . We have the estimation
1 1 n+1 q s |x| − 2 kXk∞ = 2 1− − th α + n+1 z |x| − 2 sh α 2 z ch 2 α q 1 |x| − 1 s −q s > 2 1− − = 2 − > 1. z |x| − 2 ch α z z |x| z (The last inequality can be rewritten in the form −1 < z|x|. This is true because z|x| = −2z − 1 > −1.) Thus the method is not contractive in the maximum norm in this case. Case VI (z < 0 and there exists l ∈ {1, . . . , ν − 1} such that (X)1,1 , . . . , (X)l,l > 0 and (X)l+1,l+1 , . . . , (X)ν,ν < 0). In this case our procedure is similar to the case III. Let us consider the following two functions defined on the set {1, 2, . . . , n}: n q X s g3 (i) = 2 (Γ )i,j + , z j =1 z
1 q f (i) s 1− − . g4 (i) = 2 z |x| − 2 sh((n + 1)α) z
The function g3 (i) results in the sum of the ith row of the matrix X in absolute value (without the element in the main diagonal) plus the element in the main diagonal with the opposite sign and g4 (i) gives the sum of the ith row in maximum value (without the element in the main diagonal) plus the ith element of the row. Thus if i 6 l then g3 (i) < g4 (i) = X i and if l < i 6 ν then Xi = g3 (i) > g4 (i). As we know from Lemma 4 and the case V these functions take their maximum values at the point i = ν. So g3 (ν) > g3 (i) and g3 (ν) > g4 (i) (i = 1, . . . , n). On the other hand, g3 (ν) 6 1 (see case IV), that is the condition kXk∞ 6 1 is fulfilled in this case as well. 2 Acknowledgements The author would like to thank István Faragó for the useful advice. References [1] O. Axelsson, T. Steihaug, Some computational aspects in the numerical solution of parabolic equations, J. Comput. Appl. Math. 4 (1978) 129–142. [2] K. Dekker, J.G. Verwer, Stability of Runge–Kutta Methods for Stiff Nonlinear Differential Equations, NorthHolland, Amsterdam, 1984. [3] I. Faragó, T. Pfeil, Preserving concavity in initial-boundary value problems of parabolic type and its numerical solution, Periodica Math. Hung. 30 (1995) 135–139.
462
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
[4] I. Faragó, On the non-negative solution of some system of linear algebraic equations and its application to the numerical solution of partial differential equations, Quaderni del Dipartimento di Matematica, Statistica, Informatica e Applicazioni, Bergamo, 1992. [5] I. Faragó, Non-negativity of the difference schemes, Pure Math. Appl. 6 (1996) 38–50. [6] I. Faragó, Qualitative properties of the numerical solution of linear parabolic problems with non-homogeneous boundary conditions, Comput. Math. Appl. 31 (1996) 143–150. [7] I. Faragó, One-step methods of solving a parabolic problem and their qualitative properties, Publications on Applied Analysis, Eötvös L. University, Department of Applied Analysis, Budapest, 1996. [8] A. Friedmann, Partial Differential Equations of Parabolic Type, Prentice-Hall, Englewood Cliffs, NJ, 1964. [9] R. Horváth, On the positivity of iterative methods, Ann. Univ. Budapest. Sect. Comp., Budapest, 1995, to appear. [10] J.F.B.M. Kraaijevanger, Maximum norm contractivity of discretization schemes for the heat equation, Appl. Numer. Math. 9 (1992) 475–492. [11] T. Pfeil, On the monotonicity in time of the solutions of linear second order homogeneous parabolic equations, Ann. Univ. Sci. Budapest. Eötvös Sect. Math. 36 (1993). [12] P. Rózsa, Linear Algebra and Its Applications, Tankönyvkiadó, Budapest, 1991 (in Hungarian). [13] A.A. Samarskii, Theory of the Difference Schemes, Nauka, Moscow, 1977. [14] J. Smoller, Shock Waves and Reaction–Diffusion Equations, Springer, Berlin, 1982. [15] V. Thomée, Finite difference methods for linear parabolic equations, in: P.G. Ciarlet, J.L. Lions (Eds.), Handbook of Numerical Analysis, North-Holland, Amsterdam, 1990.
Maximum norm contractivity in the numerical solution of the one-dimensional heat equation Róbert Horváth ∗ University of Sopron, Institute of Mathematics, Ady E. u. 5, Sopron, H-9400, Hungary
Abstract In this paper we consider the one-dimensional heat conduction equation (Friedmann, 1964). To the numerical solution of the problem we apply the so-called (σ, θ )-method (Faragó, 1996; Thomée, 1990) which unites a few numerical methods. With the choice σ = 0 we get the finite difference θ -method and the choice σ = 16 results in the finite element method with linear elements. The most important question is the choice of the suitable mesh-parameters. The basic condition arises from the condition of the convergence (Faragó, 1996; Samarskii, 1977; Thomée, 1990). Further conditions can be obtained aiming at preserving some qualitative properties of the continuous problem. Some of them are the following: non-negativity, convexity, concavity, shape preservation and contractivity in some norms (Dekker and Verwer, 1984). Now we shall study the maximum norm contractivity. There are some results in the literature for the parameter choices which guarantee this property (Kraaijevanger, 1992; Samarskii, 1977; Thomée, 1990). However these papers specialize only on the finite difference methods and give sufficient conditions. We determine the necessary and sufficient conditions related to the (σ, θ )-method. We close the paper with numerical examples. 1999 Elsevier Science B.V. and IMACS. All rights reserved. Keywords: Numerical solution; Qualitative properties; Heat equation; Maximum norm contractivity
1. Introduction Let us consider the following parabolic partial differential equation [8], the so-called one-dimensional heat conduction equation ∂u ∂ 2 u = 2, ∂t ∂ξ u(t, 0) = u(t, 1) = 0,
t > 0,
u(0, ξ ) = u0 (ξ ),
ξ ∈ [0, 1],
ξ ∈ (0, 1),
t > 0,
∗ E-mail: [email protected].
0168-9274/99/$20.00 1999 Elsevier Science B.V. and IMACS. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 9 9 ) 0 0 0 0 7 - 0
(1)
452
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
in the domain Ω = (0, ∞) × (0, 1), where u0 is a suitable smooth function for which there exists a unique solution on Ω. The function u : Ω → R is called the solution of the problem (1) if u ∈ C 1,2 (Ω) ∩ C(Ω) and it satisfies the equations in (1). This problem has some characteristic qualitative properties as the physical model [3,11,14]. Applying the well-known maximum principle [14] we can show that the function t 7→ maxx∈[0,1] |u(t, x)| is monotonically decreasing, that is the absolute value of the temperature does not increase in time. This property is called maximum norm contractivity. It can be shown that the function t 7→ u(t, x) tends to 0 for all x ∈ [0, 1] if t → ∞. (If u0 is a concave or a convex function then the convergence is monotone with respect to x.) It is said that the solution composes itself. Moreover, we say that the problem (1) preserves the non-negativity if u0 > 0 implies that u > 0 in Ω. Among the above properties we shall deal with the maximum norm contractivity. For a numerical solution to the problem (1) we define the following mesh for fixed τ ∈ R+ and h = 1/(n + 1) (n ∈ N) [13,15]:
Ωτ,h := (tj , xi ) ∈ Ω | xi = ih (i = 0, . . . , n + 1), tj = j τ (j ∈ N) . Denoting the approximation to the exact solution at the j th time level by y j ∈ Rn and using the so-called (σ, θ)-method [7] we define the following one-step iteration process: y j +1 − y j 1 1 (2) = −θ 2 Qy j +1 − (1 − θ) 2 Qy j , j = 0, 1, . . . . τ h h Here Q = tridiag[−1, 2, −1] ∈ Rn×n and M = I − σ Q ∈ Rn×n are matrices (the matrix I denotes the unit matrix), θ ∈ [0, 1] and σ ∈ [0, 14 ) are arbitrary parameters and the vector y 0 is a suitable approximation of the initial function u0 . This method unites a few numerical methods. If σ = 0 then we get the finite difference method and in the case of σ = 16 the method (2) results in the finite element method with linear elements. First of all we reformulate our problem to a problem of linear algebra. Let us introduce the notations q = τ/ h2 , z = θq − σ . If z = 0 then the method (2) is called explicit. Let us denote the matrix I − qQ by X. If z 6= 0 then we call the method implicit and introduce the notations s = q − z, x = 2 + 1/z, X 1 = Q + (x − 2)I and X 2 = I − sQ. We notice that s > 0 as s = q(1 − θ) + σ and |x| > 2 because z > − 14 . In this case X denotes the matrix (1/z)X−1 1 X 2 . The matrix X is well-defined because the matrix X 1 is invertible for all values |x| > 2 [12]. Thus the iteration (2) can be rewritten in the form M
y j +1 = Xy j ,
j = 0, 1, . . . .
(3)
It is an interesting question to examine the behaviour of the numerical solution defined by the iteration (3). We seek conditions for both the method and the mesh under which the qualitative properties of the continuous problem are preserved in the numerical solution [3–7,9,10]. We require the iteration (3) to be composing, too. This means that the vector-sequence y j (j = 0, 1, . . .) has to tend to the zero vector for all initial vectors y 0 , that is the matrix X must be a so-called convergent matrix. In this case we call the method (2) weak convergent method. We get the following condition for the weak convergence. Lemma 1 [7]. Assume that σ ∈ [0, 14 ) and θ ∈ [0, 0.5). Then the (σ, θ)-method is weak convergent for fixed n if and only if the condition 1 1 − 4σ cos2 (π/(2(n + 1))) 2(1 − 2θ) cos2 (π/(2(n + 1))) is satisfied. Moreover, it is weak convergent for all n if and only if the condition 1 − 4σ q6 2(1 − 2θ) q6
(4)
(5)
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
453
holds. For the parameters σ ∈ [0, 14 ) and θ ∈ [0.5, 1] the method is always weak convergent. Naturally in this paper we consider schemes satisfying (4) and we analyze the properties of the discrete maximum norm contractivity. Definition 2. We call the numerical process (2) contractive in the maximum norm if the inequalities
j +1
y
∞
6 y j ∞
are fulfilled for each j ∈ N. We note that the maximum norm contractivity is not required for neither the stability nor the convergence of the method, it is desirable only from qualitative point of view. Obviously the method (3) is contractive in the maximum norm if and only if the maximum norm of the matrix X is not greater than 1, that is if the condition kXk∞ 6 1 is fulfilled. Our aim is to give conditions for the values of the mesh-parameters by a fixed (σ, θ)-method under which the numerical solution is contractive in maximum norm. 2. The condition of the contractivity In this section we give the condition of the maximum norm contractivity both by fixed n and for all values n. Let us consider the case of z = 0. It follows from the form of the matrix X, that if n = 1 or n = 2, then the necessary and sufficient condition of the contractivity is q 6 1 or q 6 23 , respectively. If n > 3 then the condition is q 6 12 . Now let us pass over the implicit methods. Let G ∈ Rn×n be a symmetric, one-pair matrix defined as follows: (G)i,j = (Γ )i,j if x > 2 and (G)i,j = (−1)i+j −1 (Γ )i,j if x < −2, where sh(iα) sh((n + 1 − j )α) |x| γi,j , if i 6 j , (Γ )i,j = γi,j = , α = arch . (6) γj,i , if i > j , sh α sh((n + 1)α) 2 For our purpose we rewrite the matrix X to an explicit form. The lemma in question can be stated as follows. Lemma 3. Assume that z 6= 0. For the matrix X defined by (1/z)X−1 1 X 2 we have the relation
1 q X= G − sI . z z
(7)
Proof. The matrix G is the inverse of the matrix Q + (x − 2)I [12]. Multiplying the matrix X by the matrix z · (Q + (x − 2)I ) from the left we obtain the matrix I − sQ. 2 To find the maximum norm of the matrix X we have to examine the structure of the matrix and to sum the absolute values of the elements in each row. Using the notation ν = [(n + 1)/2] (the integer part of (n + 1)/2) we can state the following: Lemma 4. For the diagonal elements of the matrix G are valid the following relations: γi,i = γn+1−i,n+1−i ,
i = 1, . . . , ν,
γ1,1 < γ2,2 < · · · < γν,ν .
(8)
454
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
Proof. To prove this lemma it is sufficient to analyze the continuous function ω 7→ sh(ωα) sh((n + 1 − ω)α) on the interval [1, n]. 2 Introducing the notation Gi =
Pn
j =1 |(G)i,j |
=
Pn
j =1 (Γ )i,j
we have
Lemma 5. The sum of the modulus of the elements in the ith row of the matrix G can be expressed as follows: 1 Gi = (9) (1 − γ1,i − γi,n ), i = 1, . . . , n. |x| − 2 Proof. Using the notation bz = 1/(|x| − 2) we have the identity
T
Q + (|x| − 2)I · bz, . . . , bz
T
= (1, . . . , 1)T + bz, 0, . . . , 0, bz .
Multiplying (10) by the matrix Γ we obtain the statement.
(10)
2
Because of Lemma 4 for the sign of z and for the sign distribution of diagonal elements of the matrix X the possible cases are the following: Case I. z > 0 and (X)i,i > 0, i = 1, . . . , n. Case II. z > 0 and (X)i,i < 0, i = 1, . . . , n. Case III. z > 0 and there exists k ∈ {1, . . . , ν −1} such that (X)1,1 , . . . , (X)k,k < 0 and (X)k+1,k+1 , . . . , (X)ν,ν > 0. Case IV. z < 0 and (X)i,i < 0, i = 1, . . . , n. Case V. z < 0 and (X)i,i > 0, i = 1, . . . , n. Case VI. z < 0 and there exists l ∈ {1, . . . , ν − 1} such that (X)1,1 , . . . , (X)l,l > 0 and (X)l+1,l+1 , . . . , (X)ν,ν < 0. For the different cases we examine the condition of the maximum norm contractivity. Our results can be summarized as follows. Lemma 6. Assume that the method (2) is weak convergent. Then in the case V the maximum norm contractivity is not fulfilled. In cases I, III, IV and VI the contractivity is valid without further restrictions for the parameters. Finally, in the case II we get the inequality
sh((n + 1 − ν)α) sh(να) sh((n + 1 − ν)α) + sh(να) + 2(|x|−2) q s sh α 1− + 61 z sh((n + 1)α) z by fixed n and the condition 8σ (θ − 1) + 2 − θ +
p
(8σ (θ − 1) + 2 − θ)2 − 16(1 − θ)2 σ (4σ − 1) 8(1 − θ)2 for all n ∈ N which guarantee the contractivity. In case of θ = 1 no restrictions arise. q6
(11)
(12)
Proof. Because of the simple but long calculations of the proof, the steps are not detailed here. For the detailed proof we refer to Appendix. 2 It is remains to examine the six cases above, by which conditions for the parameters are realized. Comparing the bounds, we can give a necessary and sufficient condition of the maximum norm
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
455
contractivity. We are not able to give the conditions by fixed n because most of our inequalities are implicit. To see some numerical results by fixed n, we refer to Section 3. We can formulate the most important result of our paper as follows: Theorem. Let θ ∈ (0, 1) and σ ∈ [0, 14 ) be fixed. The method (2) is both weak convergent and contractive in maximum norm for all values of n if and only if θ > 2σ and the parameter q satisfies the following condition: σ 8σ (θ − 1) + 2 − θ + 6q 6 θ
p
(8σ (1 − θ) − 2 + θ)2 − 16(1 − θ)2 σ (4σ − 1) . 8(1 − θ)2
(13)
In case of θ = 1 the condition is q > σ and in the case θ = 0 (consequently σ = 0) the condition is 0 6 q 6 12 . Proof. The proof is based on Lemma 6. First let us consider the case z < 0. According to Lemma 6 the weak convergent method (2) is contractive if and only if there exists a negative element in the main diagonal of the matrix X (in cases IV and VI). From this condition we get the inequality (X)ν,ν < 0. This is fulfilled for all n if it is fulfilled for n = 1. It yields the condition q > (1 − 2σ )/(2(1 − θ)). We have to compare this bound with the bounds z < 0 and (5). From the first comparison we get the necessary condition θ < 2σ . On the other hand, in this case 1 − 2σ 1 − 4σ > 2(1 − θ) 2(1 − 2θ) so the conditions 1 − 2σ q> 2(1 − θ)
and
q6
1 − 4σ 2(1 − 2θ)
cannot be fulfilled simultaneously. Thus in the case z < 0 there is no parameter choice σ, θ, q which yields both weak convergence and maximum norm contractivity for all n. Now let us examine the case z > 0. Because of the conditions (5) and z > 0 we get the necessary condition θ > 2σ . In cases I and III the contractivity is fulfilled automatically. These cases are realized if and only if (X)ν,ν > 0. This is fulfilled for all n if it is fulfilled for n = 1. It yields the condition q 6 (1 − 2σ )/(2(1 − θ)). Comparing this condition with the condition θ > 2σ we obtain the inequality σ/θ < q 6 (1 − 2σ )/(2(1 − θ)). (If θ = 1 then q > σ .) If q > (1 − 2σ )/(2(1 − θ)) then there exists a natural number n0 that for all n < n0 the elements in the main diagonal of the matrix X are negative. But in this case we got the following condition for q in Lemma 6: q6
8σ (θ − 1) + 2 − θ +
p
(8σ (1 − θ) − 2 + θ)2 − 16(1 − θ)2 σ (4σ − 1) . 8(1 − θ)2
Thus in the case of z > 0, the conditions of the contractivity are θ > 2σ and σ 8σ (θ − 1) + 2 − θ +
p
(8σ (1 − θ) − 2 + θ)2 − 16(1 − θ)2 σ (4σ − 1) . 8(1 − θ)2
(If θ = 1 then q > σ .) Let us take into consideration the explicit methods as well. In this case we have got the condition q 6 12 . Thus the theorem is proved. 2
456
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
To analyze the conditions by fixed σ we introduce the following functions: lσ (θ) = cσ (θ) =
1 − 4σ 2(1 − 2θ)
θ<
1 2
pσ (θ) =
,
8σ (θ − 1) + 2 − θ +
p
1 − 2σ 2(1 − θ)
(θ < 1),
eσ (θ) =
σ θ
(8σ (1 − θ) − 2 + θ)2 − 16(1 − θ)2 σ (4σ − 1) 8(1 − θ)2
(θ > 0), (θ < 1).
They denote the barrier functions of the convergence. By the other values of θ let us define the values of the functions as ∞. Then the following statements are true: (a) Let σ 6= 0. Then lσ (2σ ) = cσ (2σ ) = pσ (2σ ) = eσ (2σ ) = 12 , if θ > 2σ then lσ (θ) > cσ (θ) > pσ (θ) > eσ (θ) and if θ 6 2σ then lσ (θ) 6 cσ (θ) 6 pσ (θ) 6 eσ (θ). (b) Let σ = 0. Then l0 (θ) 6 c0 (θ) 6 p0 (θ). Remark. Let us apply the theorem for two well-known methods. In the case of σ = 0 we have the finite difference method. Then the method (2) is contractive in the maximum norm for all n if and only if 06q 6
2−θ . 4(1 − θ)2
(14)
The choice θ = 0 corresponds to the forward Euler method (in this case the condition is q 6 0.5), θ = 0.5 corresponds to the Crank–Nicolson method (in this case the condition is q 6 1.5) and the choice θ = 1 corresponds to the backward Euler method (in this case q is optional). We have got the upper bound which was obtained by Kraaijevanger [10] with other tools. The method of that reference does not make possible the examination of the contractivity bound depending on n. Our method as it is shown in Section 3 does carry out this.
Fig. 1.
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
457
The other special choice is σ = 16 which is the finite element method with linear elements. Then the method (2) is contractive in maximum norm for all n if and only if θ > 13 and the condition √ 1 θ + 2 + 9θ 2 − 12θ + 12 (15) 6q 6 6θ 24(1 − θ)2 holds. (If θ = 1 then q > 16 .) The dark area shows the possible parameter choices in Fig. 1. Remark. According to the continuous case, the method (2) is called non-negativity preserving, if the vectors y j (j = 1, 2, . . .) are non-negative for all non-negative initial vectors y 0 . Considering the iteration (3) we can see that the necessary and sufficient condition of the non-negativity preservation is the inequality X > 0. This condition is fulfilled if and only if z > 0 and (X)i,i > 0 (i = 1, 2, . . . , n). It is apparent from the proof of Lemma 6 that the non-negativity preservation is a sufficient condition of the maximum norm contractivity. 3. Numerical examples Let us assume that (σ, θ) are fixed. To define the mesh on which this method is contractive in the maximum norm by fixed n, we have to solve the implicit inequality (11). In our examples we execute it numerically. The necessary and sufficient condition of both the maximum norm contractivity and the weak convergence by fixed n, θ, σ is the following: (a) Assume that σ 6= 0. If θ > 2σ then cσ,n (θ) > q > eσ (θ). Conversely, if θ < 2σ then eσ (θ) > q > pσ,n (θ) and q 6 lσ,n (θ). (b) In the case σ = 0 we have the condition 0 6 q 6 c0,n (θ). Here cσ,n (θ) denotes such a value of q for which the relation (11) turns into the equality. On the base of the proof of Lemma 6 we can verify that by n → ∞ cσ,n (θ) → cσ (θ) and cσ,n (θ) forms a monotonically decreasing sequence in n. Thus we can give for all fixed n a greater bound for q than cσ (θ) defined in [10]. Table 1 The case of σ = 0 n
θ = 0.00
θ = 0.30
θ = 0.45
θ = 0.50
1
1
2.5
10
∞
3
0.5
1.0763
2.9894
990237
5
0.5
0.8961
1.5035
2.0000
7
0.5
0.8721
1.3334
1.6180
9
0.5
0.8682
1.2952
1.5352
11
0.5
0.8675
1.2850
1.5113
13
0.5
0.8674
1.2822
1.5037
15
0.5
0.8674
1.2813
1.5012
25
0.5
0.8673
1.2810
1.5000
c0 (θ )
0.5
0.8673
1.2810
1.5000
458
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
Table 2 σ = 16 , θ = 0.1 n
p1/6,n (0.1)
l1/6,n (0.1)
1
0.3704
0.8333
2
0.2941
0.4167
3
0.2810
0.3156
4
0.2784
0.2743
5
0.2779
0.2533
7
0.2778
0.2331
9
0.2778
0.2240
∀n
c1/6 (0.1) = 0.2778
l1/6 (0.1) = 0.2083
Table 3 σ = 16 , θ = 0.5 n
3
5
7
9
11
15
25
∀n
c1/6,n (0.5)
2.3333
0.9658
0.9034
0.8973
0.8957
0.8954
0.8954
0.8954
Moreover, lσ,n (θ) denotes the bound derived from (4) and pσ,n (θ) is the solution of the equality X ν,ν = 0 for q. It can be verified that pσ,1 (θ) = pσ (θ) and by n → ∞ both pσ,n (θ) → cσ (θ) and lσ,n (θ) → lσ (θ) monotonically decreasing. Because of pσ (θ) > lσ (θ) if θ < 2σ there does not exist q which yields both contractivity and weak convergence for all n. First of all let us consider the case σ = 0. Table 1 shows the computed values c0,n (θ) in the cases θ = 0, 0.3, 0.45, 0.50. We remark that these values satisfy the condition of the weak convergence as well. Now let us consider the case σ = 16 . We must consider two different cases according to θ 6 2σ or θ > 2σ . Let θ be for example 0.1 < 2σ = 23 . Then by fixed n the conditions of the maximum norm contractivity and the weak convergence are satisfied if and only if 53 > q > p1/6,n (0.1) and q 6 l1/6,n (0.1). Table 2 shows that weak convergence and contractivity are fulfilled only in the cases n = 1, 2, 3. Now let θ be 0.5 > 13 . Then by fixed n the conditions of the maximum norm contractivity and weak convergence are satisfied if and only if 13 6 q 6 c1/6,n(0.5). Table 3 shows the values c1/6,n (0.5) computed with numerical methods. Appendix Proof of Lemma 6. Let us denote by X i (i = 1, . . . , n) the sum of the elements in absolute value in the ith row of the matrix X. The values Xi using Lemma 3 can be expressed in the following form: Xi =
n n q q X q X s (Γ ) + (X) = (Γ ) + (G) − , i,j i,i i,j i,i z2 z2 j =1 z2 j =1 z j 6=i
j 6=i
i = 1, . . . , n.
(A.1)
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
459
Our aim is to compute the value of kXk∞ and to give the condition kXk∞ 6 1. We have established in Lemma 6 that there exist only six sign-distributions in the matrix X. We examine these cases. We apply the expression (A.1) and we use Lemmas 3 and 5 in the calculations. Case I (z > 0 and (X)i,i > 0, i = 1, . . . , n). In this case we get the the following estimation: Xi =
n n q X q s s q X (Γ ) + (Γ ) − (Γ )i,j − = i,j i,i 2 2 2 z j =1 z z z j =1 z j 6=i
q 1 s q s q −s z = 2 (1 − γ1,i − γi,n ) − = (1 − γ1,i − γi,n ) − 6 = = 1. z |x| − 2 z z z z z The inequality above shows that the condition z > 0 and (X)i,i > 0 (i = 1, . . . , n) is a sufficient condition for the maximum norm contractivity. Case II (z > 0 and (X)i,i < 0, i = 1, . . . , n). In this case we have Xi =
n n q X q s 2q s q X (Γ ) − (Γ ) + (Γ )i,j − 2 (Γ )i,i + = i,j i,i 2 2 2 z j =1 z z z j =1 z z j 6=i
q 1 f (i) s = 2 1− + , z |x| − 2 sh((n + 1)α) z where 2(|x| − 2) sh (n + 1 − i)α sh(iα), Df = {1, 2, . . . , n}. sh α Extending the function f on the set R we get a differentiable function. With differentiation we can easily obtain that this function has a local minimum point at i = (n + 1)/2 and it has local maximum points at
f (i) = sh (n + 1 − i)α + sh(iα) +
n+1 1 1 n+1 i= ± arch sh α sh α . 2 α 2 2 We verify that the value X i (i = 1, . . . , n) takes its maximum value by the choice i = [(n + 1)/2] = ν. Thus we have to show that f (1) > f (ν). We prove this statement for odd values of n. In this case ν = (n + 1)/2 and the inequality can be rewritten in the following form:
2(|x| − 2) n+1 2(|x| − 2) 2 n + 1 sh α + sh(nα) + sh(nα) sh α > 2 sh α + sh α . sh α 2 sh α 2 Using the identities of hyperbolic functions we get the inequality
n+1 n−1 n+1 n+1 (|x| − 2) sh sh2 α ch α − sh α > α sh α sh(nα) . 2 2 2 sh α 2 This is fulfilled by the choice n = 1. Let us analyze both functions standing on the right and on the left hand sides, respectively. As earlier we assume that they are functions on R. Since
n+1 (|x| − 2) sh (n + 1)α − 2 sh α ch(nα) (A.2) α > 2 sh α therefore the derivative of the left hand side is greater than the derivative of the right hand side. To prove (A.2) it is sufficient to use the well-known identities for the hyperbolic functions.
ch(nα) − ch
460
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
We have proved that in case II, the condition of the maximum norm contractivity can be defined as follows: q f (ν) s kXk∞ = 1− + 6 1. (A.3) z sh((n + 1)α) z We seek the condition which guarantees the contractivity for all n. Let n be an odd natural number. Then the condition (A.3) yields q 1 n+1 1 s |x| − 2 1− + 6 1. th α − n+1 2 z |x| − 2 sh α 2 z ch( 2 α) This is satisfied for all n if and only if q 1 |x| − 2 s 1− + 6 1. 2 z |x| − 2 sh α z Substituting the values of the parameters we get the condition
4q 2 (1 − θ)2 + q 8σ (1 − θ) − 2 + θ + 4σ 2 − σ 6 0 from which we obtain the condition p 8σ (θ − 1) + 2 − θ + (8σ (1 − θ) − 2 + θ)2 − 16(1 − θ)2 σ (4σ − 1) q6 . (A.4) 8(1 − θ)2 If θ = 1 then no condition arises. Case III (z > 0 and there exists k ∈ {1, . . . , ν − 1} such that (X)1,1 , . . . , (X)k,k < 0 and (X)k+1,k+1 , . . . , (X)ν,ν > 0). In this case we consider the following two functions on the set {1, 2, . . . , n}: n q X s g1 (i) = 2 (Γ )i,j − , z j =1 z
1 q f (i) s 1− + . g2 (i) = 2 z |x| − 2 sh((n + 1)α) z
The function g1 (i) results in the sum of the ith row of the matrix X and g2 (i) gives the sum of the ith row, taking the element in the main diagonal with opposite sign. Thus if i 6 k then g1 (i) < g2 (i) = Xi and if k < i 6 ν then X i = g1 (i) > g2 (i). As we know from Lemma 4 and case II these functions take their maximum values at the point i = ν. So g1 (ν) > g1 (i) and g1 (ν) > g2 (i) (i = 1, . . . , n). On the other hand, g1 (ν) 6 1, that is the condition kXk∞ 6 1 is fulfilled in this case as well. Case IV (z < 0 and (X)i,i < 0, i = 1, . . . , n). In this case we can calculate as follows: n n q X s q X s q Xi = 2 (Γ )i,j + − 2 (G)i,i = 2 (Γ )i,j + z j =1 z z z j =1 z j 6=i
q 1 sh((n + 1 − i)α) + sh(iα) s = 2 1− + z |x| − 2 sh((n + 1)α) z q 1 s 4s − 1 2(2(s − z) − 1) 6 2 + = =1+ . z |x| − 2 z 4z + 1 4z + 1 Because of z > − 14 we have to prove the relation 2(s − z) 6 1. Obviously this inequality is fulfilled in the case θ > 12 . In case of θ < 12 we get the upper bound q6
1 − 4σ . 2(1 − 2θ)
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
461
This condition is satisfied automatically if the matrix X is convergent. Case V (z < 0 and (X)i,i > 0, i = 1, . . . , n). In this case the expression X i can be reformulated as follows: Xi =
n n q X q s 2q s q X (Γ ) − (Γ ) − (Γ )i,j − 2 (Γ )i,i − = i,j i,i 2 2 2 z j =1 z z z j =1 z z j 6=i
sh((n + 1 − i)α) sh(iα) sh((n + 1 − i)α) + sh(iα) + 2(|x|−2) q 1 s sh α = 2 1− − . z |x| − 2 sh((n + 1)α) z X i takes its maximum value at i = ν on the set {1, . . . , n}, that is kXk∞ = X ν . We have the estimation
1 1 n+1 q s |x| − 2 kXk∞ = 2 1− − th α + n+1 z |x| − 2 sh α 2 z ch 2 α q 1 |x| − 1 s −q s > 2 1− − = 2 − > 1. z |x| − 2 ch α z z |x| z (The last inequality can be rewritten in the form −1 < z|x|. This is true because z|x| = −2z − 1 > −1.) Thus the method is not contractive in the maximum norm in this case. Case VI (z < 0 and there exists l ∈ {1, . . . , ν − 1} such that (X)1,1 , . . . , (X)l,l > 0 and (X)l+1,l+1 , . . . , (X)ν,ν < 0). In this case our procedure is similar to the case III. Let us consider the following two functions defined on the set {1, 2, . . . , n}: n q X s g3 (i) = 2 (Γ )i,j + , z j =1 z
1 q f (i) s 1− − . g4 (i) = 2 z |x| − 2 sh((n + 1)α) z
The function g3 (i) results in the sum of the ith row of the matrix X in absolute value (without the element in the main diagonal) plus the element in the main diagonal with the opposite sign and g4 (i) gives the sum of the ith row in maximum value (without the element in the main diagonal) plus the ith element of the row. Thus if i 6 l then g3 (i) < g4 (i) = X i and if l < i 6 ν then Xi = g3 (i) > g4 (i). As we know from Lemma 4 and the case V these functions take their maximum values at the point i = ν. So g3 (ν) > g3 (i) and g3 (ν) > g4 (i) (i = 1, . . . , n). On the other hand, g3 (ν) 6 1 (see case IV), that is the condition kXk∞ 6 1 is fulfilled in this case as well. 2 Acknowledgements The author would like to thank István Faragó for the useful advice. References [1] O. Axelsson, T. Steihaug, Some computational aspects in the numerical solution of parabolic equations, J. Comput. Appl. Math. 4 (1978) 129–142. [2] K. Dekker, J.G. Verwer, Stability of Runge–Kutta Methods for Stiff Nonlinear Differential Equations, NorthHolland, Amsterdam, 1984. [3] I. Faragó, T. Pfeil, Preserving concavity in initial-boundary value problems of parabolic type and its numerical solution, Periodica Math. Hung. 30 (1995) 135–139.
462
R. Horváth / Applied Numerical Mathematics 31 (1999) 451–462
[4] I. Faragó, On the non-negative solution of some system of linear algebraic equations and its application to the numerical solution of partial differential equations, Quaderni del Dipartimento di Matematica, Statistica, Informatica e Applicazioni, Bergamo, 1992. [5] I. Faragó, Non-negativity of the difference schemes, Pure Math. Appl. 6 (1996) 38–50. [6] I. Faragó, Qualitative properties of the numerical solution of linear parabolic problems with non-homogeneous boundary conditions, Comput. Math. Appl. 31 (1996) 143–150. [7] I. Faragó, One-step methods of solving a parabolic problem and their qualitative properties, Publications on Applied Analysis, Eötvös L. University, Department of Applied Analysis, Budapest, 1996. [8] A. Friedmann, Partial Differential Equations of Parabolic Type, Prentice-Hall, Englewood Cliffs, NJ, 1964. [9] R. Horváth, On the positivity of iterative methods, Ann. Univ. Budapest. Sect. Comp., Budapest, 1995, to appear. [10] J.F.B.M. Kraaijevanger, Maximum norm contractivity of discretization schemes for the heat equation, Appl. Numer. Math. 9 (1992) 475–492. [11] T. Pfeil, On the monotonicity in time of the solutions of linear second order homogeneous parabolic equations, Ann. Univ. Sci. Budapest. Eötvös Sect. Math. 36 (1993). [12] P. Rózsa, Linear Algebra and Its Applications, Tankönyvkiadó, Budapest, 1991 (in Hungarian). [13] A.A. Samarskii, Theory of the Difference Schemes, Nauka, Moscow, 1977. [14] J. Smoller, Shock Waves and Reaction–Diffusion Equations, Springer, Berlin, 1982. [15] V. Thomée, Finite difference methods for linear parabolic equations, in: P.G. Ciarlet, J.L. Lions (Eds.), Handbook of Numerical Analysis, North-Holland, Amsterdam, 1990.